tract-ac-estimates.R

library(tigris)
library(mapgl)
library(freestiler)
library(tidyverse)
options(tigris_use_cache = TRUE)

ac_estimates <- read_csv(
  "https://www2.census.gov/programs-surveys/demo/datasets/lace/2023/LACE_23_Tract.csv"
) |>
  mutate(GEOID = str_sub(GEO_ID, start = -11))

us_tracts <- tracts(cb = TRUE, year = 2023)

us_tracts_ac <- us_tracts |>
  select(GEOID) |>
  inner_join(
    ac_estimates,
    by = "GEOID"
  ) |>
  filter(NO_AC_PE >= 0, AC_PE >= 0)

no_ac_class <- interpolate_palette(
  data = us_tracts_ac,
  column = "NO_AC_PE",
  method = "equal",
  n = 5,
  palette = viridisLite::inferno
)

# Visualize them all - GeoJSON source with tolerance 0 (a bit slower to load)
maplibre(
  style = openfreemap_style("positron"),
  center = c(-98.5, 39.5),
  zoom = 4
) |>
  add_source(
    id = "tracts-source",
    data = us_tracts_ac,
    tolerance = 0
  ) |>
  add_fill_layer(
    id = "no-ac",
    source = "tracts-source",
    fill_color = no_ac_class$expression,
    fill_opacity = 0.7,
    popup = concat(
      "<strong>",
      get_column("NAME"),
      "</strong><br>",
      "Occupied housing units: ",
      number_format("HSE_OCC_E", maximum_fraction_digits = 0),
      "<br>With AC: ",
      number_format("AC_PE", maximum_fraction_digits = 1),
      "%",
      " (±",
      number_format("AC_PM", maximum_fraction_digits = 1),
      ")",
      "<br>Without AC: ",
      number_format("NO_AC_PE", maximum_fraction_digits = 1),
      "%",
      " (±",
      number_format("NO_AC_PM", maximum_fraction_digits = 1),
      ")"
    ),
    hover_options = list(fill_color = "cyan", fill_opacity = 1)
  ) |>
  add_legend(
    "Occupied housing units without AC (%)",
    type = "continuous",
    layer_id = "no-ac",
    interactive = TRUE,
    values = no_ac_class$breaks,
    colors = no_ac_class$colors,
    position = "bottom-left",
    width = 300
  )

# Faster: tile first, and visualize fast!
freestile(
  us_tracts_ac,
  "us_tracts_ac.pmtiles",
  layer_name = "tracts",
  min_zoom = 2,
  max_zoom = 12
)

# Serve the tiles
serve_tiles("us_tracts_ac.pmtiles", port = 8888)

# Visualize with PMTiles:
maplibre(
  style = openfreemap_style("positron"),
  center = c(-98.5, 39.5),
  zoom = 4
) |>
  add_pmtiles_source(
    id = "tracts-source",
    url = "http://localhost:8888/us_tracts_ac.pmtiles",
    tolerance = 0
  ) |>
  add_fill_layer(
    id = "no-ac",
    source = "tracts-source",
    source_layer = "tracts",
    fill_color = no_ac_class$expression,
    fill_opacity = 0.7,
    popup = concat(
      "<strong>",
      get_column("NAME"),
      "</strong><br>",
      "Occupied housing units: ",
      number_format("HSE_OCC_E", maximum_fraction_digits = 0),
      "<br>With AC: ",
      number_format("AC_PE", maximum_fraction_digits = 1),
      "%",
      " (±",
      number_format("AC_PM", maximum_fraction_digits = 1),
      ")",
      "<br>Without AC: ",
      number_format("NO_AC_PE", maximum_fraction_digits = 1),
      "%",
      " (±",
      number_format("NO_AC_PM", maximum_fraction_digits = 1),
      ")"
    ),
    hover_options = list(fill_color = "cyan", fill_opacity = 1)
  ) |>
  add_legend(
    "Occupied housing units without AC (%)",
    type = "continuous",
    layer_id = "no-ac",
    interactive = TRUE,
    values = no_ac_class$breaks,
    colors = no_ac_class$colors,
    position = "bottom-left",
    width = 300
  )