Calculate percentage overlap of 2 sets of polygons in R

Question

I have a map of NC counties. I have overlayed two circles, representing the coverage area of two radio transmitters in the state:

### Make transmitters df
rm(list = ls())
library(tidyverse)
library(sf)


transmitters <-
  data.frame(
    ID = 1:2,
    longitude = c(-80.631974, -77.808488),
    latitude = c(35.838583, 35.526252),
    radius = c(50, 100))

### Turn transmitters df into an sf object
tr_wgs <- 
  st_as_sf(transmitters, 
           coords = c("longitude", "latitude"), 
           crs = 4326, 
           dim = "XY")

### Set units (a projection that preserves area in NC)
tr_wgs <- 
  st_transform(tr_wgs, 
               "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD27 +units=km")

### Create a buffer (not sure of the original km measurements, so right now it is just making a 100km buffer)
tr_buff <- 
  st_buffer(tr_wgs, c(transmitters$radius))

### Read North Carolina (and coerce to be same projection as tr_buff)
nc <-
  st_read(system.file("shape/nc.shp", package="sf")) %>%
  st_transform(., "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD27 +units=km")

### Plot
ggplot() +
  geom_sf(data = nc) +
  geom_sf(data = tr_buff)

What I now want to calculate is: what percentage of each of the 100 counties is covered by the circle. I have tried the following:

intersect_pct <- st_intersection(nc, tr_buff)

But this just gives me the intersection areas. How do I get a "proportion covered" variable attached to each county in the nc dataframe?

Answer 1

You're almost there! You can calculate the area of your new intersecting shapes, and then merge this back into your nc object:

# Calculate area and tidy up
intersect_pct <- st_intersection(nc, tr_buff) %>% 
   mutate(intersect_area = st_area(.)) %>%   # create new column with shape area
   dplyr::select(NAME, intersect_area) %>%   # only select columns needed to merge
   st_drop_geometry()  # drop geometry as we don't need it

Now you have a simple dataframe with the county name and size of the intersection with the buffer (by default should be in km2 units). Merge it into nc and create a new column that calculates coverage (I also recalculate county area in case the AREA value is wrong):

# Create a fresh area variable for counties
nc <- mutate(nc, county_area = st_area(nc))

# Merge by county name
nc <- merge(nc, intersect_pct, by = "NAME", all.x = TRUE)

# Calculate coverage
nc <- nc %>% 
   mutate(coverage = as.numeric(intersect_area/county_area))