puma-mapping-acs1.R

library(tidycensus)
library(tigris)
library(tidyverse)
options(tigris_use_cache = TRUE)

# Grab the data
us_income <- get_acs(
  geography = "puma",
  variables = "B19013_001",
  year = 2023,
  survey = "acs1",
  geometry = TRUE
) 

# Simple ggplot with shifted geometry
ggplot(shift_geometry(us_income)) +
  geom_sf(aes(fill = estimate), color = NA) +
  scale_fill_viridis_c(option = "plasma", name = "Median Income") +
  theme_void() +
  labs(title = "Median Household Income by PUMA")

# Interactive map with mapgl
library(mapgl)

# Format the popup
popup_content <- glue::glue(
  "<strong>{us_income$NAME}</strong><br>",
  "Median Income: {scales::dollar_format()(us_income$estimate)}"
)

us_income$popup <- popup_content

# Build the interactive map
mapboxgl(
  style = mapbox_style("light"), 
  center = c(-98.5795, 39.8283), 
  zoom = 3
) %>%
  add_fill_layer(
    id = "puma_income",
    source = us_income,
    fill_color = interpolate(
      column = "estimate",
      values = c(25000, 50000, 100000, 150000, 200000),
      stops = viridisLite::plasma(5),
      na_color = "lightgrey"
    ),
    fill_opacity = 0.7,
    popup = "popup", 
    hover_options = list(
      fill_color = "cyan",
      fill_opacity = 1
    )
  ) %>%
  add_legend(
    "Median HH Income by PUMA",
    values = c("$25k", "$50k", "$100k", "$150k", "$200k"),
    colors = viridisLite::plasma(5)
  )