[R-sig-Geo] raster::zonal with more than 1 zonal layer

[Loïc Dutrieux]

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') )
> full_join(u,b,by=c('newVal'='zone’))
>
> 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)
plot(s)

unique(s)

So yes, it looks like unique combinations in the 3rd dimension is what 
unique() returns in the case of a multilayer raster object.

Cheers,
Loïc

>
> ds
>
>
>
> 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
>     elevation
>      > 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
>
>     Cheers,
>     Loïc
>
>      > 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>
>     https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>