6

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?

1
  • 1
    Please don't put rm(list = ls()) in your scripts. library(tidyverse) is also a bad idea - it pulls in loads of packages that you aren't using which can confuse things. If you need constituent packages from tidyverse then include them individually. In this case I don't think you need it at all since sf brings in the pipe from magrittr.
    – Spacedman
    May 21, 2020 at 11:59

1 Answer 1

12

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))

Your Answer

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

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