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