I have a polygon sf object of 125 unique polygons (representing different areas scattered throughout a US state). None of the polygons share a border. I also have a 30m hillshaded digital elevation model raster image of the state, which covers the area of the polygons (and then some). I want to use the terrain function from the raster package to calculate the unique slope and aspect for each polygon. My end product would be having the polygon sf object with two new columns, slope and aspect, so each polygon would have its slope and aspect.
I have used the terrain function on the original, state-wide raster image and a masked raster image of just the polygon areas, which returns the min and max slope and aspects for the whole raster image.
Beyond that, I am pretty lost on how to get the slope and aspect calculated for each unique polygon's area, and then get that data into the polygon sf object. If anyone has any advice I'd be happy to hear it.
Thank you!
Example data
library(raster)
elevation <- getData('alt', country='CHE')
x <- terrain(elevation, opt=c('slope', 'aspect'), unit='degrees')
Solution
sw <- getData("GADM", country="CHE", level=1)
e <- extract(x, sw, fun=mean, na.rm=TRUE)
x <- data.frame(sw$NAME_1, e)
head(x)
# sw.NAME_1 slope aspect
#1 Aargau 3.065123 170.2449
#2 Appenzell Ausserrhoden 4.905536 187.5591
#3 Appenzell Innerrhoden 7.380191 170.1497
#4 Basel-Landschaft 3.742552 178.2771
#5 Basel-Stadt 1.002250 125.0466
#6 Bern 9.025566 195.5117
You can assign the values to the SpatialPolygons like this
swe <- cbind(sw, e)