[R-sig-Geo] Creating a maximum raster layer from RasterStack
Agustin Lobo
alobolistas at gmail.com
Tue Oct 15 10:06:17 CEST 2013
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