Friday, 19 August 2016

dem - Zonal Statistics for Polygons with R


I'd like to do zonal statistics / terrain analysis with polygon overlays - I achieved to get min / max like so:


library(raster)
library(rgdal)

alt <- getData('alt', country = "AT")
gadm <- getData('GADM', country = "AT", level = 2)
gadm_sub <- gadm[1:3, ]


plot(alt)
plot(gadm_sub, add=T)

extract(alt, gadm_sub, fun = max, na.rm = T, small = T, df = T)

But for slope / aspect I get this error which I don't grasp how to prevent:


extract(alt, gadm_sub, fun = function(x) terrain(x, opt = "aspect"), na.rm = T, small = T, df = T)
# > Fehler in FUN(X[[1L]], ...) : unbenutzte(s) Argument(e) (na.rm = TRUE)

Answer



Sry - that one was too easy.. For the record I'll post it here:



    library(raster)
library(rgdal)

alt <- getData('alt', country = "AT")
gadm <- getData('GADM', country = "AT", level = 2)
gadm_sub <- gadm[1:3, ]

plot(alt)
plot(gadm_sub, add=T)


asp <- terrain(alt, opt = "aspect", unit = "degrees", df = F)
slo <- terrain(alt, opt = "slope", unit = "degrees", df = F)

extract(slo, gadm_sub, fun = mean, na.rm = T, small = T, df = T)
extract(asp, gadm_sub, fun = mean, na.rm = T, small = T, df = T)

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...