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

[Vijay Lulla]

Nice use of do.call Loic (Sorry, don't know how to do the accents).  I
wasn't aware that you can send in a list to args!  So, thanks for it.

Now here's my quesiton.  If I change the statement

l <- list(s1,s2,s3,s4,fun=mean)
to
l <- list(s1=s1,s2=s2,s3=s3,s4=s4,fun=mean) # similar to Thiago's models.list

I get an error with missing "overlay" function for signature
"missing","missing".  I tried to find some ways around it but have had
no luck.  This is due to R's S3/S4 object system and how function
arguments are passed, I think.  I don't understand the S4 system that
well.  Basically, how would I use a list with names in do.call with
raster::overlay?

On Wed, Dec 2, 2015 at 11:26 AM, Loïc Dutrieux <loic.dutrieux at wur.nl> wrote:
> 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
>>
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo