[R-sig-Geo] Creating a maximum raster layer from RasterStack

[Agustin Lobo]

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