Creating a maximum raster layer from RasterStack

library(raster)

r1 <- raster(matrix((1:9000000),3000,3000))
r2 <- raster(matrix(sample(1:9000000,9000000),3000,3000))
r3 <- raster(matrix(sample(1:9000000,9000000),3000,3000))
s <- stack(r1,r2,r3)

system.time(mx <- stackApply(s,indices=rep(1, nlayers(s)),max))

system.time(mx3 <- max(s, na.rm=TRUE))

s <- stack(r1,r2,r3,r1,r2,r3,r1,r2,r3,r1,r2,r3,r1,r2,r3,r1,r2,r3,r1,r2,r3,r1,r2,r3)
system.time(mx <- stackApply(s,indices=rep(1, nlayers(s)),max))

system.time(mx3 <- max(s))

Robert,
I keep on because this is an important point for the analysis of time
series of image:

Actually I get stackApply() faster than max() if the nb of layers increases:

Example with 3 larger layers:
r1 <- raster(matrix((1:9000000),3000,3000))
r2 <- raster(matrix(sample(1:9000000,9000000),3000,3000))
r3 <- raster(matrix(sample(1:9000000,9000000),3000,3000))
s <- stack(r1,r2,r3)
system.time(mx <- stackApply(s,indices=c(1,1,1),max))
   user  system elapsed
  2.052   0.508   2.566
system.time(mx3 <- max(s, na.rm=TRUE))
   user  system elapsed
  1.892   0.400   2.300

But with 24 layers:

s <- stack(r1,r2,r3,r1,r2,r3,r1,r2,r3,r1,r2,r3,r1,r2,r3,r1,r2,r3,r1,r2,r3,r1,r2,r3)
system.time(mx <- stackApply(s,indices=c(1,1,1),max))
   user  system elapsed
 17.896   0.744  18.694
system.time(mx3 <- max(s, na.rm=TRUE))
   user  system elapsed
 37.948  17.436  55.519


Agus

On Tue, Oct 15, 2013 at 7:06 PM, Robert J. Hijmans <r.hijmans at gmail.com> wrote:

Agus, yes, that was a mistake, thank you for pointing that out. And
yes, you can use stackApply (just as you could always use tapply
instead of apply), but there are simpler, more direct, and faster,
approaches such as

which.max(s)
max(s)

calc(s, max)

Robert


On Tue, Oct 15, 2013 at 1:06 AM, Agustin Lobo <alobolistas at gmail.com> wrote:

Eddie:

Please note "You need one element in "indices" for each input layer"
and do read the help page.

Robert: yo made a typo right? it is not "indices=1:nlayers(s)", but
"indices=rep(1,nlayers(s))" ?

i.e:
r1 <- raster(matrix((1:9),3,3))
r2 <- raster(matrix(sample(1:9,9),3,3))
r3 <- raster(matrix(sample(1:9,9),3,3))
s <- stack(r1,r2,r3)
as.array(s)
mx <- stackApply(s,indices=c(1,1,1),max)
plot(round(mx))
as.matrix(mx)

I understand Eddie wants one single layer with the maximum value for
each pixel across layers
and the layer with the actual maximum:
miwhichmax <- function(x,na.rm=TRUE) which.max(x)
mxd <- stackApply(s,indices=c(1,1,1),miwhichmax)
as.matrix(mxd)

Regarding performance, you know better of course.

Agus


On Mon, Oct 14, 2013 at 7:06 PM, Robert J. Hijmans <r.hijmans at gmail.com> wrote:

Eddie,

stackApply is not the appropriate function in this case (stackApply is
equivalent to tapply or aggregate for a matrix; see the manual),
although you could make it do what you want trough
"indices=1:nlayers(s)", but that is rather inefficient. Instead, to
get the max value of cells across layers, you can do:

s <- stack(system.file("external/rlogo.grd", package="raster"))
m <- max(s, na.rm=TRUE)

or

m2 <-  calc(s, max, na.rm=TRUE)

Robert



On Mon, Oct 14, 2013 at 9:13 AM, Eddie Smith <eddieatr at gmail.com> wrote:

Thanks again Agustin.

The code works fine for 5 raster layers.

library(raster)
img <- list.files(pattern='\\.img$')
s <- stack(img)
mx <- stackApply(s,indices=c(1,1,1,1,1),max)

class       : RasterLayer
dimensions  : 1032, 1656, 1708992  (nrow, ncol, ncell)
resolution  : 0.04166667, 0.04166667  (x, y)
extent      : 96.5, 165.5, -16.5, 26.5  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
+towgs84=0,0,0
data source : in memory
names       : layer
values      : 12.48, 34.35  (min, max)


Just wondering if I have hundreds of raster layer in a folder.

I tried mx <- stackApply(s,indices=c(1:5),max) but it gave me different
result

class       : RasterBrick
dimensions  : 1032, 1656, 1708992, 5  (nrow, ncol, ncell, nlayers)
resolution  : 0.04166667, 0.04166667  (x, y)
extent      : 96.5, 165.5, -16.5, 26.5  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
+towgs84=0,0,0
data source : in memory
names       : layer.1, layer.2, layer.3, layer.4, layer.5
min values  :   16.75,   16.26,   15.00,   12.48,   15.00
max values  :   34.35,   34.02,   32.88,   33.48,   34.05

Anyone?

Thanks in advance.







On Tue, Oct 8, 2013 at 9:17 PM, Agustin Lobo <alobolistas at gmail.com> wrote:

Eddie,
(I'm answering through the list, perhaps somebody else is interested)

You need one element in "indices" for each input layer, so if you have
5 layers you should use
mx <- stackApply(s,indices=c(1,1,1,1,1),max)
The actual value of each index indicates to which output layer that
input layer is contributing.
In your case, there is only one output layer, thus all values must be
1. But see the
help page for an example with different index values.

Agus

On Tue, Oct 8, 2013 at 9:41 PM, Agustin Lobo <Agustin.Lobo at ictja.csic.es>
wrote:

---------- Forwarded message ----------
From: Eddie Smith <eddieatr at gmail.com>
Date: Tue, Oct 8, 2013 at 7:40 PM
Subject: Re: [R-sig-Geo] Creating a maximum raster layer from RasterStack
To: Agustin.Lobo at ictja.csic.es


Dear Agus,

Thank you very much for your reply. Your script below is working
however, would you mind if I ask a basic question on what is the
purpose of [indices=c(1,1)] ? I tried this code with 5 raster layers
and the error message is

mx <- stackApply(s,indices=c(1,1),max)

Warning message:
In stackApply(s, indices = c(1, 1), max) :
  number of items to replace is not a multiple of replacement length

I am sorry if my question is too basic. I am new to R.

Thank you.

library(raster)
r1 <- raster(matrix(rnorm(100,5,10),10,10))
r2 <- raster(matrix(rnorm(100,2,30),10,10))
s <- stack(r1,r2)
mx <- stackApply(s,indices=c(1,1),max)
miwhichmax <- function(x,na.rm=TRUE) which.max(x)
mxd <- stackApply(s,indices=c(1,1),miwhichmax)


On Mon, Oct 7, 2013 at 7:20 PM, Agustin Lobo <alobolistas at gmail.com>

wrote:

oops, copied the wrong line,
it must be
mx <- stackApply(s,indices=c(1,1),max)

Agus


On Mon, Oct 7, 2013 at 8:18 PM, Agustin Lobo <alobolistas at gmail.com>

wrote:

What about this:
r1 <- raster(matrix(rnorm(100,5,10),10,10))
r2 <- raster(matrix(rnorm(100,2,30),10,10))
s <- stack(r1,r2)
mx <- stackApply(s,max)

which.max() fails because raster requires the na.rm argument (I

think, Robert?)

mxd <- stackApply(s,indices=c(1,1),which.max)
Error in FUN(newX[, i], ...) : unused argument (na.rm = TRUE)

so make your own function icluding the na.rm argument:

miwhichmax <- function(x,na.rm=TRUE) which.max(x)
mxd <- stackApply(s,indices=c(1,1),miwhichmax)

Agus


On Thu, Oct 3, 2013 at 6:55 PM, Eddie Smith <eddieatr at gmail.com>

wrote:

Dear list,

I have a RasterLayer that contains 365 bands with a dimension of

1032,

1656, 1708992  (nrow, ncol, ncell).

library(raster)
img <- list.files(pattern='\\.img$')
stack <- stack(img)
data2002 <- writeRaster(stack, filename="2002.tif", format="GTiff",
overwrite=TRUE)

1. How could I produce a raster layer that contains a maximum value

from

365 bands for each pixel?
I tried maxValue() but that only give me a maximum value for each

band.

2. If I manage to do step 1, is it possible to know from which band

is

actually the maximum value originate?

        [[alternative HTML version deleted]]

--
--
Dr. Agustin Lobo
Institut de Ciencies de la Terra "Jaume Almera" (CSIC)
Lluis Sole Sabaris s/n
08028 Barcelona
Spain
Tel. 34 934095410
Fax. 34 934110012
e-mail Agustin.Lobo at ictja.csic.es

        [[alternative HTML version deleted]]