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

Edzer Pebesma edzer.pebesma at uni-muenster.de
Wed Dec 2 20:06:30 CET 2015


This is a little off-topic for this list, but I usually do something like:

> library(raster)
Loading required package: sp
> overlay # what is overlay? Ah, it is an S4 method:
standardGeneric for "overlay" defined from package "raster"

function (x, y, ...)
standardGeneric("overlay")
<bytecode: 0x41c5a48>
<environment: 0x41d9d48>
Methods may be defined for arguments: x, y
Use  showMethods("overlay")  for currently available ones.
> showMethods("overlay") # which methods are available?
Function: overlay (package raster)
x="Raster", y="missing"
x="Raster", y="Raster"

> getMethod("overlay", c("Raster", "missing")) # a function called
Method Definition:

function (x, y, ...)
... # rest omitted

Not nearly as nice as julia, thanks for the pointer to @which!

On 02/12/15 19:48, Vijay Lulla wrote:
> Thanks Edzer.  Makes sense!  Is there an easier way to inspect what
> specific function will be resolved to with a particular call other
> than "getMethod"?  For e.g., considering our discussion, how would I
> have used getMethod (or some other method) to learn how
> raster::overlay is being called (i.e. actual arguments and what the
> final call would look like)?  Julia has a great macro @which which
> does the resolution and shows you what specific method will be called.
> Check out @which 1 + 2 vs @which 1 + 2 + 3 vs @which 1 + 2 + 3 + 4 vs
> @which 1 + (2 + (3 + 4)) on a julia interpreter.  We need something
> like this for R, especially for spatial objects.
> 
> Regardless, thanks for your explanation.  It's very helpful.
> Cordially,
> Vijay.
> 
> On Wed, Dec 2, 2015 at 12:50 PM, Edzer Pebesma
> <edzer.pebesma at uni-muenster.de> wrote:
>>
>>
>> On 02/12/15 18:28, Loïc Dutrieux wrote:
>>> 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.
>>
>> That is right: the signature of raster::overlay is
>>
>>      overlay(x, y, ..., fun, filename="", recycle=TRUE)
>>
>> and the do.call call with the named list resolves in a call to overlay
>> of the form
>>
>>       overlay(s1 = s1, s2 = s2, s3 = s3, s4 = s4, fun = mean)
>>
>> which has missing x and y. Not naming the arguments passes s1 to x, s2
>> to y and the others to ... allowing overlay to choose the right method.
>>
>>>
>>> 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
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> R-sig-Geo at r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>> --
>> Edzer Pebesma
>> Institute for Geoinformatics (ifgi),  University of Münster,
>> Heisenbergstraße 2, 48149 Münster, Germany; +49 251 83 33081
>> Journal of Statistical Software:   http://www.jstatsoft.org/
>> Computers & Geosciences:   http://elsevier.com/locate/cageo/
>> Spatial Statistics Society http://www.spatialstatistics.info
>>
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

-- 
Edzer Pebesma
Institute for Geoinformatics (ifgi),  University of Münster,
Heisenbergstraße 2, 48149 Münster, Germany; +49 251 83 33081
Journal of Statistical Software:   http://www.jstatsoft.org/
Computers & Geosciences:   http://elsevier.com/locate/cageo/
Spatial Statistics Society http://www.spatialstatistics.info

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 490 bytes
Desc: OpenPGP digital signature
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20151202/e5463db8/attachment.bin>


More information about the R-sig-Geo mailing list