[R-sig-Geo] Creating a maximum raster layer from RasterStack
Agustin Lobo
alobolistas at gmail.com
Wed Oct 16 08:36:28 CEST 2013
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