[R-sig-Geo] Creating a maximum raster layer from RasterStack
Robert J. Hijmans
r.hijmans at gmail.com
Tue Oct 15 19:06:34 CEST 2013
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
More information about the R-sig-Geo
mailing list