[R-sig-Geo] Take mean of list of raster stacks

[Loïc Dutrieux]

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
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>