Use sf to cut multipolygons by regional multipolygon

sf_cutregions.R

#################################################-
## Subdivide multipolygon by regional multipolygon with sf ----
## W.K. Petry
#################################################-
## Preliminaries ----
#################################################-
library(sf)
library(tidyverse)
library(cowplot)

theme_set(theme_bw())
#################################################-
## Read in/transform spatial projection ----
#################################################-
# use North Carolina counties example dataset
nc <- st_read(system.file("shape/nc.shp", package="sf")) %>%
  st_transform(crs = 26917) %>%
  mutate(lon = map_dbl(geometry, ~st_centroid(.x)[[1]]),
         lat = map_dbl(geometry, ~st_centroid(.x)[[2]]))

# fetch 115th congressional districts to use as the "regions" for splitting county polygons
# from https://www.census.gov/geo/maps-data/data/cbf/cbf_cds.html
temp <- tempfile()
download.file("http://www2.census.gov/geo/tiger/GENZ2017/shp/cb_2017_us_cd115_500k.zip", temp)
tempd <- tempdir()
unzip(temp, exdir = tempd)
cd <- read_sf(paste0(tempd, "/cb_2017_us_cd115_500k.shp")) %>%
  st_transform(crs = 26917) %>%
  filter(STATEFP == "37")

# plot to check that everything came in correctly
ggplot()+
  geom_sf(data = nc)+
  geom_sf(data = cd, aes(fill = CD115FP), alpha = 0.5)+
  geom_text(data = nc, aes(x = lon, y = lat, label = NAME))

#################################################-
## Split county polygons by census district ----
#################################################-
# Keep an eye on Cumberland and Bladen counties in the south -- they should end up with >1 entries each
nc_cd <- st_intersection(nc, cd) %>%
  mutate(lon = map_dbl(geometry, ~st_centroid(.x)[[1]]),
         lat = map_dbl(geometry, ~st_centroid(.x)[[2]]))

origplot <- ggplot()+
  geom_sf(data = nc)+
  geom_text(data = nc, aes(x = lon, y = lat, label = NAME), size = 2)+
  coord_sf(datum = NA)

cutplot <- ggplot()+
  geom_sf(data = nc_cd, aes(fill = CD115FP))+
  geom_text(data = nc_cd, aes(x = lon, y = lat, label = NAME), size = 2)+
  coord_sf(datum = NA)

plot_grid(origplot, cutplot, ncol = 1)

#################################################-
## Subset by congressional district ----
#################################################-
sub_nc <- tibble(cd = unique(cd$CD115FP)) %>%
  mutate(district_mpoly = map(.x = cd, .f = ~nc_cd %>% filter(CD115FP == .x)))

ggplot(data = sub_nc %>% filter(cd == "07") %>% pull(district_mpoly) %>% .[[1]])+
  geom_sf(aes(fill = CD115FP))+
  geom_text(aes(x = lon, y = lat, label = NAME), size = 2)+
  coord_sf(datum = NA)