I have an existing vector shapefile with all the counties in the United States (called counties.shp) I also have a dataset that contains many species-county presence data points (species.csv) which looks something like this:

 GEOID  county    state         Species 
1: 21007 Ballard Kentucky     plant1
2: 21007 Ballard Kentucky     plant2
3: 21007 Ballard Kentucky     plant3
4: 21007 Ballard Kentucky     plant4 
5: 21017 Bourbon Kentucky     plant5

I'd like to create a new polygon for each species that includes all of the counties that they are paired with in the species.csv file. Each row is a unique species-county combination, so both counties and species are repeated throughout that file, yet I'd like to get a single polygon file for each species. There are several thousand species, so I'd like to have a way to do this all at once and compile the resulting polygons in a relatively neat way so I can easily use them for zonal statistics.

I'm imagining this is a somewhat simple task, but I'm completely new to GIS in R.

  • what do you mean for "new polygon", a multipolygon of the counties? a convex hull? a new file?
    – Elio Diaz
    Apr 16, 2021 at 2:05
  • Do you know the basics of spatial data in R, using the sf package? eg how to read shapefiles, plot them, inspect them etc?
    – Spacedman
    Apr 16, 2021 at 7:51
  • The first two rows you listed don't look a "unique species-county" - is row 2 meant to say "plant2"?
    – Spacedman
    Apr 16, 2021 at 7:52
  • @ElioDiaz a multipolygon is what I'm after.
    – pfadenhw
    Apr 16, 2021 at 12:12
  • @Spacedman I've been using rgdal instead of sf but I've been able to read shapefiles and rasters and plot them using that. And yes, that was a typo, thanks for the catch I've corrected the original.
    – pfadenhw
    Apr 16, 2021 at 12:12

1 Answer 1


This approach filters the county polygons and them puts them together using lapply to iterate for each plant


nc = st_read(system.file("shape/nc.shp", package="sf"))

plants = c("p1", "p2", "p3", "p4")
# sample the counties twice to get more than one plant by county
plants_by_county = data.frame(county = sample(nc$NAME, 200, replace = T), plant =  sample(plants, 200, replace = T))

# data frame looks like this
          county plant
1         Yancey    p4
2          Davie    p2
3       Buncombe    p2
4       Alamance    p1

plant_multipolygons = lapply(plants, function(x) {
  nc %>%
    filter(NAME %in% plants_by_county$county[plants_by_county$plant == x]) %>% # filter counties
    select(geometry) %>%
    summarise() %>% # comment out summarise to get single polygons
    mutate(plant = x) })

plant_multipolygons = do.call(rbind, plant_multipolygons) # bind them into a data.frame

tm_shape(plant_multipolygons) + tm_polygons() + tm_facets("plant")

enter image description here

  • 1
    you must transform your spatialPolygons object to sf with st_as_sf(your_polygon_object)
    – Elio Diaz
    Apr 16, 2021 at 16:44
  • 1
    pass dplyr::select() to instantiate properly
    – Elio Diaz
    Apr 16, 2021 at 16:55
  • 1
    you have 2 options: stick with sf which calculates st_area from a EPSG:4326 CRS and returns ellipsoidal area in m²; or use a projected CRS such as UTM or lambert
    – Elio Diaz
    Apr 16, 2021 at 17:37
  • 1
    I think sticking with this process should work as long as I can figure out how to make sure my rasters are in the same CRS. Thank you for all your help!
    – pfadenhw
    Apr 16, 2021 at 17:41
  • 1
    I just updated R (now version 4.1.0) and all of a sudden this code is running VERY SLOWLY - it used to take about 3 minutes to run, now it takes upwards of 2 hours... Any idea why this happened?
    – pfadenhw
    Jun 30, 2021 at 14:50

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.