[R-sig-Geo] raster::zonal with more than 1 zonal layer
Dominik.Schneider at colorado.edu
Fri Mar 4 17:36:47 CET 2016
On Mar 3, 2016 11:15 AM, "Loïc Dutrieux" <loic.dutrieux at wur.nl> wrote:
> On 03/03/2016 05:19 PM, Dominik Schneider wrote:
>> That looks great, thanks!
>> to close the question:
>> b=as.data.frame(zonal(dat, zoneLayer, 'mean') )
>> 1 question,
>> in this line: u <- data.frame(unique(s))
>> does unique(s) give all unique combinations in the 3rd dimension of the
>> stack? that’s a neat trick.
> I'm not 100% sure, we can check with a simplified case.
> r1 <- r2 <- raster()
> n <- ncell(r1)
> r1 <- c(rep(1, n/3), rep(2, n/3), rep(3, n/3))
> r2 <- c(rep(4, n/2), rep(5, n/2))
> s <- stack(r1, r2)
> So yes, it looks like unique combinations in the 3rd dimension is what
> unique() returns in the case of a multilayer raster object.
>> On Wed, Mar 02, 2016 at 3:24 AM "Loïc Dutrieux" <">"Loïc Dutrieux"
>> <mailto:>> wrote:
>> On 03/02/2016 01:48 AM, Dominik Schneider wrote:
>> > I'd like to summarise a raster using elevation and watershed. I was
>> > originally using extract() with a shape file and then each
>> elevation band
>> > within each polygon but it's very slow. zonal() is much faster
>> and I can
>> > rasterize my polygons to use it.
>> Is it really much faster? I sounds very similar to what extract() does
>> under the hood.
>> > But how do I robustly combine multiple
>> > rasterized shapefiles?
>> > e.g.
>> > dat=raster(matrix(runif(64),nrow=8))
>> > z1=raster(matrix( sample(1:4, 64, replace=T),nrow=8))#represents
>> > bands
>> > z2=raster(matrix(sample(1401:1408,64,replace=T),nrow=8))#represents
>> > watersheds
>> > zonal(dat,z1*z2,'mean')
>> > this works well if you are certain that each combination of the
>> values in
>> > z1 and z2 are unique and each combination is present. otherwise
>> it gets
>> > messy. are there any suggestions for this use case?
>> > ideally one could do: zonal(dat,stack(z1,z2),'mean') and all the
>> > bookkeeping would be taken care of.
>> What about:
>> s <- stack(z1, z2)
>> u <- data.frame(unique(s))
>> u$newVal <- seq(nrow(u))
>> zoneLayer <- subs(s, u, by = 1:2, which = 3)
>> # Let's validate
>> a <- zonal(dat,z1*z2,'mean')
>> b <- zonal(dat, zoneLayer, 'mean')
>> all.equal(b[order(b[,2]),2], a[order(a[,2]),2]) # TRUE
>> > my other thought is to extract all the
>> > values into data frames and use dplyr but I was wondering if
>> there was a
>> > raster way to do this.
>> > Thanks
>> > Dominik
>> > [[alternative HTML version deleted]]
>> > _______________________________________________
>> > R-sig-Geo mailing list
>> > R-sig-Geo at r-project.org <mailto: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 <mailto:R-sig-Geo at r-project.org>
[[alternative HTML version deleted]]
More information about the R-sig-Geo