[R-sig-Geo] Take mean of list of raster stacks
Vijay Lulla
vijaylulla at gmail.com
Wed Dec 2 19:48:56 CET 2015
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
More information about the R-sig-Geo
mailing list