Take mean of list of raster stacks

Hi Thiago,

Otherwise overlay wrapped in do.call seems to work.

See example below:

Cheers,
Lo?c

library(raster)

# Create multiple raterStacks
s1 <- stack(system.file("external/rlogo.grd", package="raster"))
s2 <- s1 * 2
s3 <- s1 * 3
s4 <- s1 * 4

# Create a unnamed list (not sure why it didn't work with a named list)
l <- list(s1, s2, s3, s4, fun = mean)

# do.call call
do.call(what = raster::overlay, args = l)




On 12/02/2015 03:39 PM, lyndon.estes at gmail.com wrote:

Hi Thiago,




Done in haste, but I think this might do it (it?s on an 8X8 problem
though):





stlist <- lapply(1:8, function(x) {

  rl <- stack(lapply(1:8, function(y) {

    r <- raster(nrow = 10, ncol = 10)

    r[] <- sample(1:100, size = ncell(r), replace = TRUE)

    r

  }))

  rl

})

names(stlist) <- paste0("s", 1:8)




stack(lapply(1:8, function(x) {

   calc(stack(lapply(1:8, function(y) stlist[[y]][[x]])), mean)

}))




Hope this helps.




Cheers, Lyndon




?
Sent from Mailbox

On Wed, Dec 2, 2015 at 9:23 AM, Thiago V. dos Santos
<thi_veloso at yahoo.com.br> wrote:

Hi all,
I have a list with five raster stacks, each of them containing 9 layers:

models.list

$CanESM2
class : RasterBrick
dimensions : 23, 19, 437, 9 (nrow, ncol, ncell, nlayers)
resolution : 0.5, 0.5 (x, y)
extent : -57.5, -48, -34, -22.5 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names : layer.1, layer.2, layer.3, layer.4, layer.5, layer.6, layer.7,
layer.8, layer.9
min values : 137.51260, 103.75805, 85.07232, 114.59114, 88.59638,
82.38541, 98.64818, 91.78697, 74.57888
max values : 526.1966, 490.5268, 537.6004, 536.0648, 526.3977, 509.5332,
557.7880, 503.1330, 531.5689
$`GFDL-ESM2M`
class : RasterBrick
dimensions : 23, 19, 437, 9 (nrow, ncol, ncell, nlayers)
resolution : 0.5, 0.5 (x, y)
extent : -57.5, -48, -34, -22.5 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names : layer.1, layer.2, layer.3, layer.4, layer.5, layer.6, layer.7,
layer.8, layer.9
min values : 99.87192, 84.49617, 82.04732, 91.23503, 82.46968, 78.61677,
91.31480, 84.72990, 77.58408
max values : 549.9278, 550.9575, 555.1746, 542.2581, 526.9369, 543.8348,
532.9768, 524.7191, 525.7651
$inmcm4
class : RasterBrick
dimensions : 23, 19, 437, 9 (nrow, ncol, ncell, nlayers)
resolution : 0.5, 0.5 (x, y)
extent : -57.5, -48, -34, -22.5 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names : layer.1, layer.2, layer.3, layer.4, layer.5, layer.6, layer.7,
layer.8, layer.9
min values : 153.1610, 180.0696, 165.8414, 155.4981, 210.9747, 131.2129,
205.0893, 149.3376, 164.3868
max values : 548.4998, 521.2526, 532.5670, 551.9284, 561.8148, 523.1451,
534.9090, 561.0131, 551.4501
$`MRI-CGCM3`
class : RasterBrick
dimensions : 23, 19, 437, 9 (nrow, ncol, ncell, nlayers)
resolution : 0.5, 0.5 (x, y)
extent : -57.5, -48, -34, -22.5 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names : layer.1, layer.2, layer.3, layer.4, layer.5, layer.6, layer.7,
layer.8, layer.9
min values : 206.9614, 205.4357, 173.1827, 139.5373, 169.0720, 172.5434,
195.4526, 160.2298, 182.6004
max values : 687.7671, 686.6686, 689.2235, 687.3547, 645.3307, 671.9138,
669.0936, 665.2333, 669.0399
$`NorESM1-M`
class : RasterBrick
dimensions : 23, 19, 437, 9 (nrow, ncol, ncell, nlayers)
resolution : 0.5, 0.5 (x, y)
extent : -57.5, -48, -34, -22.5 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names : layer.1, layer.2, layer.3, layer.4, layer.5, layer.6, layer.7,
layer.8, layer.9
min values : 211.6625, 185.8265, 187.7064, 187.3369, 186.3985, 149.3203,
156.6462, 153.4485, 116.1606
max values : 605.5658, 603.2598, 569.0408, 599.4353, 589.8222, 601.7283,
617.0612, 603.3071, 645.2594
What I need to do is to come out with a single stack, also with 9 layers,
that is composed by the mean of the correspondent layers of all elements in
the list.
For example, the first layer of the resulting stack would be the average
of the first layer of the five elements of the list.
In terms of code, it would be something like this:
1st layer of result stack <- mean (1st layer of 1st element, 1st layer of
2nd element, 1st layer of 3rd element, 1st layer of 4th element, 1st layer
of 5th element)
2nd layer of result stack <- mean (2nd layer of 1st element, 2nd layer of
2nd element, 2nd layer of 3rd element, 2nd layer of 4th element, 2nd layer
of 5th element)
3rd layer of result stack <- mean (3rd layer of 1st element, 3rd layer of
2nd element, 3rd layer of 3rd element, 3rd layer of 4th element, 3rd layer
of 5th element)
...
8th layer of result stack <- mean (8th layer of 1st element, 8th layer of
2nd element, 8th layer of 3rd element, 8th layer of 4th element, 8th layer
of 5th element)
9th layer of result stack <- mean (9th layer of 1st element, 9th layer of
2nd element, 9th layer of 3rd element, 9th layer of 4th element, 9th layer
of 5th element)
Any hints on how I can accomplish that?
Greetings,
-- Thiago V. dos Santos
PhD student
Land and Atmospheric Science
University of Minnesota

        [[alternative HTML version deleted]]