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

Wed Dec 2 18:28:53 CET 2015

Hi Vijay,

I was also wondering about that. The reason I think is that the first 
raster* object passed to overlay must be named x= (or unnamed)

So that
l <- list(x = s1, s2 = s2, s3 = s3, s4 = s4, fun = mean)

works but,

l <- list(s1 = s1, s2 = s2, s3 = s3, s4 = s4, fun = mean)

doesn't.

Cheers,
Loïc

On 12/02/2015 06:17 PM, Vijay Lulla wrote:
> 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