us_states_zip-codes

us_states_zip-codes

library(tidyverse)
library(tigris)
library(sf)
library(janitor)
library(tmap)

tmap_mode('view')

sf_use_s2(FALSE)

set_crs <- 4326

df <- tigris::zctas(cb = TRUE) %>% 
  st_as_sf() %>% 
  st_transform(set_crs)

states <- tigris::states(cb = TRUE) %>% 
  st_as_sf() %>% 
  st_transform(set_crs)


st_crs(states) == st_crs(df)
st_crs(states)


# Intersect ZIP code polygons and US state polygons
intersect_pct <- df %>% 
  st_intersection(., states)

# Includes duplicated ZIP-state matches
# Calculate area of overlap
intersect_pct <- intersect_pct %>% 
  mutate(int_area = st_area(.))

intersect_pct <- intersect_pct %>% 
  arrange(ZCTA5CE20)

# Check the duplicates - they are hopefully along state borders
t <- intersect_pct %>% 
  group_by(ZCTA5CE20) %>% 
  mutate(t = n()) %>% 
  filter(t > 1) %>% 
  ungroup()

ggplot() +
  geom_sf(data = t)

rm(t)

df %>% 
  filter(ZCTA5CE20 == '37380') %>% 
  tm_shape(.) +
  tm_polygons(alpha = .6)


# Clear out duplicates by selecting (within ZIP code polygons) the state match with the largest state overlap
zip_states <- intersect_pct %>% 
  mutate(int_area = unclass(int_area)) %>% 
  group_by(ZCTA5CE20) %>% 
  slice_max(int_area) %>% 
  ungroup() %>% 
  select(-int_area)


# Get their coordinates
zip_states_cents <- zip_states %>% 
  st_centroid() %>% 
  st_coordinates() %>% 
  as_tibble() %>% 
  rename(longitude_cent = X,
         latitude_cent = Y)

zip_states
zip_states_cents

zip_states <- zip_states %>% 
  bind_cols(., zip_states_cents)

rm(zip_states_cents)

# Drop geometry column
zip_states <- zip_states %>% 
  st_drop_geometry()