R Programming – Plotting and Analyzing Extracted Elevation Data in R

elevationrrasterzonal statistics

I have a shapefile of county boundaries for New York state and elevation data downloaded through the elevatr package.

library(elevatr)
library(raster)
library(sf)
library(USAboundaries)

counties <- us_counties(map_date = "1930-01-01", resolution = "high", states = c("NY"))

counties_sf <- as(counties, "Spatial")
elevation_data <-get_elev_raster(counties_sf, z=9, src = "aws")
map <- extract(elevation_data, counties)
plot(elevation_data, axes=TRUE)
plot(st_geometry(counties), add=TRUE)

This produces an image like this:

New York elevation data

That's all well and good, but how do I restrict the elevation data to just the area of the county boundaries, e.g. for computing zonal statistics or for creating a map of only the state?

The extract function from the raster package returns a sequential list of lists of elevation points, but as far as I can tell, there isn't a way to link those points back to the unique IDs of the actual counties that come from the shapefile, even using the extent function.

Ideally I'd like to work with everything using the sf (Simple Features) package, as in this question, because that's what I'm most familiar with. It's a little confusing for me to keep track of which packages return raster objects, which return sf objects, which return spatial polygon data frames, etc.

Best Answer

You can summarize values within each polygon by passing a custom function to the extract function in the raster package, or the exact_extract function in the exactextractr package (much faster, handles pixels partially covered by polygons.)

The raster::extract function expects a summarizing function with the signature function(x, na.rm), e.g.:

counties$second_lowest_point <- extract(
 elevation_data,
 counties, 
 fun=function(x, na.rm) {
   if (na.rm) { 
     sort(na.omit(x))[2]
   } else {
     sort(x)[2]
   }
 })

The exactextractr::exact_extract function expects a summarizing function with the signature function(x, w), where w is the fraction of the pixel that is covered by the polygon. Here we're taking the area-weighted mean of elevations > 400m: 

counties$mean_elevation_over_400 <- exact_extract(
  elevation_data,
  counties,
  fun=function(x, w) { weighted.mean(x[x > 400]) })
Related Question