[R-sig-Geo] Creating a maximum raster layer from RasterStack
Robert J. Hijmans
r.hijmans at gmail.com
Thu Oct 17 06:09:46 CEST 2013
Hi Agus,
Yes, you are right, in that case stackApply is faster (but note that
you had a mistake in the indices argument that further increased the
difference ---- but is another reason why I would not use stackAppy
unless I really needed it). I did not think it would be possible that
a more generic function could be faster. This helped me to find an
inefficiency that occurs when using 'max', that does not occur in
'stackApply'. I have now fixed that (raster version 2.1-59) such that
I get:
> library(raster)
Loading required package: sp
> 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))
user system elapsed
1.56 0.28 1.84
> system.time(mx3 <- max(s, na.rm=TRUE))
user system elapsed
1.25 0.29 1.54
>
> 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))
user system elapsed
13.22 2.22 16.48
> system.time(mx3 <- max(s))
user system elapsed
12.15 2.03 14.29
That is, speeds are still about the same. That is disappointing.
Thanks for pointing this out.
Robert
On Tue, Oct 15, 2013 at 11:36 PM, Agustin Lobo <alobolistas at gmail.com> wrote:
> 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]]
>>>>>> >> >>
>>>>>> >> >> _______________________________________________
>>>>>> >> >> R-sig-Geo mailing list
>>>>>> >> >> R-sig-Geo at r-project.org
>>>>>> >> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>> >
>>>>>> >
>>>>>> >
>>>>>> >
>>>>>> > --
>>>>>> > --
>>>>>> > 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]]
>>>>>
>>>>> _______________________________________________
>>>>> 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
More information about the R-sig-Geo
mailing list