# census_workshop_ref.R
#
# Reference code for MI-Support's "Using Census Bureau data products with R"
# workshop - July 24, 2026 - hosted by Ingham County Health Department
#
# Authors: Jack Jacobs (jckjcbs@umich.edu)
#          MI-Support (micom-misupport@umich.edu)

here::i_am("census_workshop_ref.R")


#### Setup ####

library(here)
library(yaml)
library(tidyverse)
library(tidycensus)  # A great library for querying the Census API
library(sf)  # To my knowledge, the primary geospatial R library
library(tigris)  # An API wrapper for Census's geospatial "TIGERweb" service
library(ggmap)  # For loading basemaps (and maybe making density maps)

# For good measure, let's load up the fake pertussis MDSS export again.
mdss <- read_csv(here("data/FAKE_pertussis_mdss.csv"), skip = 1)

# Load in our API key(s)
auth <- read_yaml(here("auth.yml"))


#### Using tidycensus ####

# First, register my API key for this session
census_api_key(auth$census_api_key)

# Our first task: let's get populations for every year that we have in our 
# MDSS query.

# Find the years in our dataset
years <- mdss |> 
  mutate(year = Diagnosis_Date |> mdy() |> year()) |> 
  pull(year) |> 
  unique()

# Wait a minute, we need to know what variables we're looking for.
# Well wait a second, we need to know what survey we're looking in.
# Ingham Co's population is greater than 65,000, so let's use ACS-1.
# Quick refresher: total population by sex
# https://api.census.gov/data/2024/acs/acs1/variables.html
vars <- c("B01001_002","B01001_026")

# Now construct the tidycensus get_acs() query
pop_2015 <- get_acs(
  geography = "county",
  variables = vars,
  year = 2015,
  state = "MI",
  county = "Ingham",
  survey = "acs1"
)

# We can do better than copy-pasting this 10 times, though. Introducing, map()
# https://purrr.tidyverse.org/reference/map.html

# To use map(), we need to define our own function
# Note to facilitator, set this as "acs1" first and let the error come
get_acs5 <- function(
  year, variables = vars,
  state = "MI", county = "Ingham"
) {
  get_acs(
    geography = "county",
    variables = variables,
    year = year,
    state = state,
    county = county,
    survey = "acs5"
  ) |> 
    mutate(year = year)
}

# Now run this function over our vector of years, then stack them all up
pop <- years |> 
  map(get_acs5) |> 
  bind_rows() |> 
  
  # Make an easier-to-read 'sex' variable
  mutate(sex = variable |> recode_values(
    vars[1] ~ "Male",
    vars[2] ~ "Female"
  )) |> 
  
  # Keep only the stuff we care about
  select(year, sex, estimate, moe)


##### Doing it live! #####

# Make a bar graph representing the populations of 65+ over time broken down
#   by the top-level Race populations, and break White down into Hispanic/not
# B01001B_001 - B01001G_031


#### Using tigris ####

# Get all Michigan's counties
mi_counties <- counties("MI")

##### First, using sf #####

# Easiest way to see your geo object
plot(mi_counties)

# Ahh, that looks kind-of bad, huh. Oh well!
# State of Michigan has a good GIS portal where you can get better shapefiles
mi_counties <- here("data/mi_counties.geojson") |> 
  st_read() |> 
  st_transform(4326)

# Let's plot in a better way, too
mi_counties |> 
  ggplot() +
  geom_sf() +
  theme_void()

# (While we're at it, here's a shapefile of Michigan's LHD jurisdictions)
mi_lhd <- here("data/MI_LHD.geojson") |> 
  st_read() |> 
  st_transform(4326)

mi_lhd |> 
  ggplot() +
  geom_sf() +
  theme_void()


##### Back to tigris #####

###### Tracts ######

ingham_tracts <- tracts("MI", "Ingham")


###### ZIP Code Tabulation Areas (ZCTAs) ######

# Grab just Ingham's boundaries
ingham <- mi_counties |> 
  filter(Name == "Ingham")

# Get all of Michigan's ZIP codes that begin with 48 or 49
zip_48_49 <- zctas(starts_with = c("48","49")) |> 
  st_transform(4326)

# Just get Ingham
zip_intersects_ingham <- zip_48_49 |> 
  filter(st_intersects(ingham, zip_48_49, sparse = FALSE)[1, ])
zip_intersection_ingham <- zip_48_49 |> 
  st_intersection(ingham)


#### Adding basemaps to your ggplot maps ####

# Find the coordinate boundaries of our geography
ingham_bbox <- st_bbox(ingham)

# Register my stadiamaps key
register_stadiamaps(auth$stadia_api_key)

# Get a basemap
ingham_basemap <- get_stadiamap(
  bbox = c(left = -84.65, bottom = 42.4128, right = -84.07, top = 42.7773),
  zoom = 11,
  maptype = "alidade_smooth"
)

# Plot this and plot the zip codes on top
ggmap(ingham_basemap) +
  geom_sf(
    data = zip_intersection_ingham,
    fill = NA,
    inherit.aes = FALSE
  ) +
  geom_sf_text(
    aes(label = GEOID20),
    data = zip_intersection_ingham,
    size = 3,
    inherit.aes = FALSE
  )


#### We'll do leaflet live if we get to it! ####