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

[Dominik Schneider]

Thanks  Loïc
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') )
>> 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
>>
>>
>>

	[[alternative HTML version deleted]]