city_halls.csv

tract_distance.R

library(tigris)
library(sf)
library(tidyverse)
options(tigris_use_cache = TRUE)
options(tigris_class = "sf")

halls <- read_csv("city_halls.csv") %>%
  st_as_sf(coords = c("X", "Y"), crs = 4269) %>%
  st_transform("+proj=eqdc +lat_0=39 +lon_0=-96 +lat_1=33 +lat_2=45 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs")

us <- unique(fips_codes$state)[1:51]

# Get all tracts in the US
us_tracts <- map(us, function(x) {
  tracts(x, cb = TRUE)
}) %>%
  rbind_tigris()

# Get CBSAs, then identify those in the halls dataset
cbsa <- core_based_statistical_areas(cb = TRUE) %>%
  filter(GEOID %in% halls$metro_id) %>%
  select(metro_id = GEOID)

# Run an inner spatial join to allocate CBSA GEOIDs to tracts
metro_tracts <- us_tracts %>%
  st_join(cbsa, join = st_within, left = FALSE) %>%
  st_transform("+proj=eqdc +lat_0=39 +lon_0=-96 +lat_1=33 +lat_2=45 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs")

# Now, iterate through each metro and get the distances - then combine them back
tract_distance <- map(cbsa$metro_id, function(x) {
  t1 <- filter(metro_tracts, metro_id == x)
  hall <- filter(halls, metro_id == x)
  if (nrow(hall) == 1) {
    t2 <- t1 %>%
      mutate(dist = as.numeric(
        st_distance(
          st_centroid(.), hall
        )
      ))
  } else if (nrow(hall) > 1) {
    # Need the closest distance here
    dist <- st_distance(st_centroid(t1), hall)
    mindist <- apply(dist, 1, min)
    t2 <- mutate(t1, dist = mindist)
  }
  t2
}) %>%
  reduce(rbind)