Take mean of list of raster stacks

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

         [[alternative HTML version deleted]]

--
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