[R-sig-Geo] Take mean of list of raster stacks
Edzer Pebesma
edzer.pebesma at uni-muenster.de
Wed Dec 2 18:50:20 CET 2015
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
-------------- 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/347ec22d/attachment.bin>
More information about the R-sig-Geo
mailing list