framingham_median_age.R

library(tidycensus)
library(tigris)
library(sf)
library(dplyr)
library(stringr)

options(tigris_use_cache = TRUE)

# For entire city:
median_age <- get_acs(geography = "place", variables = "B01002_001",
                                  state = "MA", year = 2023, survey = "acs5") |>
 filter(stringr::str_detect(NAME, "Framingham"))


# --- Get Framingham boundary (place geometry) ---
framingham_boundary <- tigris::places(state = "MA", year = 2023, class = "sf") |>
  filter(NAME == "Framingham")

# --- Get Middlesex County census block groups with median age ---
middlesex_block_group_data <- get_acs(
  geography  = "block group",
  variables  = "B01002_001",
  state      = "MA",
  county     = "Middlesex",
  year       = 2023,  # Last available year for 5-year American Community Survey
  survey     = "acs5", # 5-year American Community Survey
  geometry   = TRUE,
  cache_table = TRUE
)

# ChatGPT helped write this part of the code
# --- Keep only block groups that fall in (or touch) Framingham ---
# Option A (centroid-in-polygon; safer for edge cases):
framingham_block_group_data <- middlesex_block_group_data |>
  mutate(centroid = st_centroid(geometry)) |>
  st_as_sf() |>
  filter(lengths(st_within(centroid, framingham_boundary)) > 0) |>
  select(-centroid)

framingham_block_group_data_for_mapping <- framingham_block_group_data |>
  select(GEOID, NAME, median_age = estimate, moe, geometry)

# Basic map

# Map:
library(tmap)

# Set tmap to plot mode (instead of interactive view)
tmap_mode("plot")

tm_shape(framingham_block_group_data_for_mapping) +
  tm_fill(
    col = "median_age",
    # palette = "YlOrRd",
    palette = "purple",
    breaks = seq(25, 55, by = 5),
    title = "Median Age"
  ) +
  tm_borders(col = "gray40", lwd = 0.5) +
  tm_shape(framingham_boundary) +
  tm_borders(col = "black", lwd = 1.2) +
  tm_layout(
    title = "Median Age of Residents by Block Group, Framingham MA (ACS 2019–2023)",
    frame = FALSE,
    legend.outside = TRUE
  )

# Map with text of median age on each polygon, with black text on light polygons and white text on dark polygons. ChatGPT helped with that part

library(dplyr)

age_breaks <- seq(20, 70, by = 5)

framingham_block_group_data_for_mapping <- framingham_block_group_data_for_mapping %>%
  mutate(
    age_bin   = cut(median_age, breaks = age_breaks, include.lowest = TRUE, right = FALSE),
    bin_id    = as.integer(age_bin),
    mid_bin   = ceiling((length(age_breaks) - 1) / 2),
    text_col  = ifelse(bin_id > mid_bin, "white", "black"),
    text_col  = ifelse(is.na(text_col), "black", text_col),  # NA guard
    age_label = sprintf("%.1f", median_age)
  )

tmap_mode("view")

tm_shape(framingham_block_group_data_for_mapping) +
  tm_polygons(
    fill = "median_age",
    fill.scale  = tm_scale(values = "Purples", breaks = age_breaks),
    fill.legend = tm_legend(title = "Median Age"),
    border.col = "gray40", border.lwd = 0.5
  ) +
  # ← Keep tm_text attached to the block-group shape
  tm_text(
    text = "age_label",
    col  = "text_col",
    size = 0.8,
    options = opt_tm_text(point.label = TRUE)  # v4 replacement for auto.placement
  ) +
  # Now switch to your boundary shape
  tm_shape(framingham_boundary) +
  tm_borders(col = "black", lwd = 1.2) +
  tm_title("Median Age of Residents by Block Group, Framingham MA (ACS 2019–2023)") +
  tm_layout(frame = FALSE, legend.outside = TRUE)