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

[Dominik Schneider]

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.  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. 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]]