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


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.

  • use fasterize to create an index raster with the polygon id, then it's data frame stuff in a table of the cells/values/polygon
    – mdsumner
    Aug 27, 2018 at 2:37
  • @mdsumner I'll give the rasterize function a try (it's running now). Then it should be straightforward to convert everything into simple features. If I get it all to work, I'll post the completed code too, for future reference.
    – Michael A
    Aug 27, 2018 at 19:15
  • 1
    @mdsumner I assume "fasterize" is a slip of the fingers? Sounds like a very useful function...
    – Spacedman
    Aug 28, 2018 at 11:20
  • 1
    @mdsumner My mistake. According to the documentation, fasterize converts an sf object to a raster object, which sounds like the opposite of what I want. Can you give me a short example of how it works in this context?
    – Michael A
    Aug 28, 2018 at 15:54
  • 1
    @dbaston For now I'm trying to summarize the values within each polygon, in a way that isn't limited to just the mean, median, etc. (ideally I'd be able to apply custom functions) but longer-term I'll need to do some nearest-neighbor sort of calculations with the elevation data.
    – Michael A
    Aug 28, 2018 at 18:41

2 Answers 2


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(
 fun=function(x, na.rm) {
   if (na.rm) { 
   } else {

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(
  fun=function(x, w) { weighted.mean(x[x > 400]) })
  • Thanks, I'll look into this. I guess the Simple Features package doesn't support methods like this yet, right?
    – Michael A
    Aug 28, 2018 at 19:02
  • Both raster and extactextractr operate on sf objects from the Simple Features package. raster does so by copying the data from sf to sp objects; exactextractr operates on sf objects directly.
    – dbaston
    Aug 28, 2018 at 19:06

Following from your code I'd use fasterize to very quickly get an index of what is polygon and what is not. You don't need the grouping so even the NA identification would be enough, but if you need elevation grouped by polygon this index is already suitable.

Please note the destructive step that updates the elevatr data set in this code, just so it looks right (masked).

counties_sf$POLYID <- 1:nrow(counties_sf)
polymap <- fasterize(counties_sf, elevation_data, field = "POLYID")
## mask out elevation where there's no polygon
elevation_data[is.na(values(polymap))] <- NA

Grouping, for example:

## to get zonal stats
tibble(value = values(elevation_data), POLYID = values(polymap)) %>% 
dplyr::filter(!is.na(value)) %>% group_by(POLYID) %>% 

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.