3D aggregation of RasterBrick or RasterStack

Question

I have some data in a 3D grid identified by simple i,j,k locations (no real-world spatial information). These data are in a RasterStack right now.

b <- stack(system.file("external/rlogo.grd", package="raster"))
# add more layers
b <- stack(b,b)
# dimensions
dim(b)
[1]  77 101   6

yields 77 rows, 101 columns, 6 layers.

# upscale by 2
up <- aggregate(b,fact=2)
dim(up)
[1] 39 51  6

yields 39 rows, 51 columns, 6 layers.

Hoped-for behavior: 3 layers.

I'm looking for a method to aggregate across layers in addition to the present behavior, which is to aggregate within each layer. I'm open to other data structures, but would prefer an existing upscaling/resampling/aggregation algorithm to one I write myself.

Potentially related are http://quantitative-advice.gg.mq.edu.au/t/fast-way-to-grid-and-sum-coordinates/110/5 or the spacetime package, which assumes the layers are temporal rather than spatial, adding more complexity.

Matt · Answer 1 · 2017-04-27T21:22:51.937

I did not read the documentation correctly. To aggregate across layers:

For example, fact=2 will result in a new Raster* object with 2*2=4 times fewer cells. If two numbers are supplied, e.g., fact=c(2,3), the first will be used for aggregating in the horizontal direction, and the second for aggregating in the vertical direction, and the returned object will have 2*3=6 times fewer cells. Likewise, fact=c(2,3,4) aggregates cells in groups of 2 (rows) by 3 (columns) and 4 (layers).

It may be necessary to play with expand=TRUE vs expand=FALSE to get it to work, but this seems inconsistent (I have reported it as a bug).

score 0 · Answer 2 · answered Apr 27 '17 at 18:24

Supouse you define agg.fact variable to denote the value 2:

agg.fact  <- 2
up <- aggregate(b, fact = agg.fact)
dim(up)
[1] 39 51  6

Now we generate a table which indicates which layers will be aggregate with anothers using agg.fact:

positions <- matrix(1:nlayers(b), nrow = nlayers(b)/agg.fact, byrow = TRUE)

     [,1] [,2]
[1,]    1    2
[2,]    3    4
[3,]    5    6

And apply a function(in this case mean but could be max``,sum` or another ...) to each pair of layers

up2 <- stack(apply(positions, 1, function(x){
  mean(b[[x[1]]], b[[x[2]]])
}))

dim(up2)
[1]  77 101   3

Or if want to aggregate in 3 dimensions (choose if want aggregate 1-2d and then 3d or viceverza):

up3 <- stack(apply(positions, 1, function(x){
  aggregate(mean(b[[x[1]]], b[[x[2]]]), fact = agg.fact) #first 3d
  #mean(aggregate(b[[x[1]]], fact = agg.fact), aggregate(b[[x[2]]]), fact = agg.fact) ##first 1d-2d
}))
dim(up3)
[1] 39 51  3

3D aggregation of RasterBrick or RasterStack

2 Answers2