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

[Loïc Dutrieux]

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
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>