Tennessee’s 2023 and 2026 U.S. House districts, with tract-level
population density.
Code:
# ============================================================
# Step 0. INSTALL AND LOAD REQUIRED PACKAGES
# ============================================================
if (!require("tidyverse")) install.packages("tidyverse")
if (!require("tidycensus")) install.packages("tidycensus")
if (!require("sf")) install.packages("sf")
if (!require("leaflet")) install.packages("leaflet")
if (!require("leaflet.extras2")) install.packages("leaflet.extras2")
if (!require("kableExtra")) install.packages("kableExtra")
if (!require("htmlwidgets")) install.packages("htmlwidgets")
library(tidyverse)
library(tidycensus)
library(sf)
library(leaflet)
library(leaflet.extras2)
library(kableExtra)
library(htmlwidgets)
# ============================================================
# Step 1. FETCH 2020 CENSUS TRACT DATA (POPULATION + GEOMETRY)
# ============================================================
tract_data <- get_decennial(
geography = "tract",
state = "TN",
variables = "P2_001N",
year = 2020,
geometry = TRUE
) %>%
rename(TotalPop = value) %>%
st_transform(5070) %>% # equal-area CRS
mutate(
land_sq_m2 = as.numeric(st_area(geometry)),
land_sq_mi = land_sq_m2 / 2.58999e6,
pop_density = TotalPop / land_sq_mi
) %>%
st_transform(4326)
# ============================================================
# Step 2. SIMPLIFY TRACT GEOMETRY (RPUBS FIX)
# ============================================================
sf::sf_use_s2(FALSE)
tract_data_simplified <- st_simplify(
tract_data,
dTolerance = 0.0002,
preserveTopology = TRUE
)
sf::sf_use_s2(TRUE)
# ============================================================
# Step 3. COLOR PALETTE + POPUP
# ============================================================
pal <- colorNumeric(
palette = "viridis",
domain = tract_data_simplified$pop_density,
na.color = "transparent"
)
popup_text <- ~paste0(
"<b>", NAME, "</b><br>",
"Population: ", scales::comma(TotalPop), "<br>",
"<b>Density:</b> ",
scales::comma(round(pop_density)), " people / sq mi"
)
# ============================================================
# Step 4. LOAD & STANDARDIZE DISTRICT GEOMETRY
# ============================================================
# --- 2023 Enacted Districts ---
cd_2023 <- get_acs(
geography = "congressional district",
state = "TN",
variables = "B01001_001",
year = 2023,
survey = "acs5",
geometry = TRUE
) %>%
st_transform(4326) %>%
mutate(
district_num = as.integer(stringr::str_extract(NAME, "\\d+")),
cd_name = paste0("District ", district_num),
label = as.character(district_num)
) %>%
select(cd_name, label, district_num, geometry)
# --- 2026 Proposed Districts ---
NewDistricts <- st_read(
"NewCongressional26.shp",
quiet = TRUE
) %>%
st_transform(4326) %>%
st_make_valid() %>%
mutate(
district_num = as.integer(DISTRICT),
cd_name = paste0("District ", district_num),
label = as.character(district_num)
) %>%
select(cd_name, label, district_num, geometry)
# ============================================================
# Step 5. TRACT → DISTRICT MATCHING (OPTION A)
# ============================================================
sf::sf_use_s2(FALSE)
tract_centroids <- tract_data_simplified %>%
st_point_on_surface()
summarise_districts <- function(tracts, districts) {
tracts %>%
st_join(districts, join = st_within) %>%
st_drop_geometry() %>%
group_by(cd_name, district_num) %>%
summarise(
Mean_Density = mean(pop_density, na.rm = TRUE),
Median_Density = median(pop_density, na.rm = TRUE),
.groups = "drop"
) %>%
arrange(district_num)
}
district_2023 <- summarise_districts(tract_centroids, cd_2023)
district_2026 <- summarise_districts(tract_centroids, NewDistricts)
sf::sf_use_s2(TRUE)
# ============================================================
# Step 6. DISTRICT LABEL POINTS
# ============================================================
cd_2023_labels <- cd_2023 %>% st_point_on_surface()
cd_2026_labels <- NewDistricts %>% st_point_on_surface()
# ============================================================
# Step 7. MAP 1: TRACT DENSITY + 2023 DISTRICTS
# ============================================================
Map_2023 <- leaflet(tract_data_simplified) %>%
addProviderTiles("CartoDB.Positron", group = "Positron (Light)") %>%
addProviderTiles("Esri.WorldStreetMap", group = "Street Map") %>%
addProviderTiles("Esri.WorldImagery", group = "Satellite") %>%
addPolygons(
fillColor = ~pal(pop_density),
fillOpacity = 0.75,
color = "#444444",
weight = 0.3,
popup = popup_text,
group = "Census Tracts"
) %>%
addPolylines(
data = cd_2023,
color = "black",
weight = 2,
group = "Congressional Districts"
) %>%
addLabelOnlyMarkers(
data = cd_2023_labels,
label = ~label,
labelOptions = labelOptions(
noHide = TRUE,
direction = "center",
textsize = "12px",
style = list("font-weight" = "bold")
),
group = "District Labels"
) %>%
addLegend(
pal = pal,
values = ~pop_density,
title = "Population density<br>(people per sq mi)",
position = "topright"
) %>%
addLayersControl(
baseGroups = c("Positron (Light)", "Street Map", "Satellite"),
overlayGroups = c("Census Tracts", "Congressional Districts", "District Labels"),
options = layersControlOptions(position = "bottomleft", collapsed = TRUE)
)
Map_2023
# ============================================================
# Step 8. TABLE 1: 2023 DISTRICTS (SORTED BY DISTRICT NUMBER)
# ============================================================
Table_2023 <- district_2023 %>%
mutate(
`Mean Density` = scales::comma(round(Mean_Density)),
`Median Density` = scales::comma(round(Median_Density))
) %>%
select(
District = cd_name,
`Mean Density`,
`Median Density`
) %>%
kbl(
format = "html",
caption = "Tract‑Level Urbanization by Congressional District (2023 Enacted Plan)"
) %>%
kable_styling(
full_width = FALSE,
bootstrap_options = c("striped", "hover", "condensed")
)
Table_2023
# ============================================================
# Step 9. MAP 2: TRACT DENSITY + 2026 DISTRICTS
# ============================================================
Map_2026 <- leaflet(tract_data_simplified) %>%
addProviderTiles("CartoDB.Positron", group = "Positron (Light)") %>%
addProviderTiles("Esri.WorldStreetMap", group = "Street Map") %>%
addProviderTiles("Esri.WorldImagery", group = "Satellite") %>%
addPolygons(
fillColor = ~pal(pop_density),
fillOpacity = 0.75,
color = "#444444",
weight = 0.3,
popup = popup_text,
group = "Census Tracts"
) %>%
addPolylines(
data = NewDistricts,
color = "black",
weight = 2,
group = "Congressional Districts"
) %>%
addLabelOnlyMarkers(
data = cd_2026_labels,
label = ~label,
labelOptions = labelOptions(
noHide = TRUE,
direction = "center",
textsize = "12px",
style = list("font-weight" = "bold")
),
group = "District Labels"
) %>%
addLegend(
pal = pal,
values = ~pop_density,
title = "Population density<br>(people per sq mi)",
position = "topright"
) %>%
addLayersControl(
baseGroups = c("Positron (Light)", "Street Map", "Satellite"),
overlayGroups = c("Census Tracts", "Congressional Districts", "District Labels"),
options = layersControlOptions(position = "bottomleft", collapsed = TRUE)
)
Map_2026
# ============================================================
# Step 10. TABLE 2: 2026 DISTRICTS (SORTED BY DISTRICT NUMBER)
# ============================================================
Table_2026 <- district_2026 %>%
mutate(
`Mean Density` = scales::comma(round(Mean_Density)),
`Median Density` = scales::comma(round(Median_Density))
) %>%
select(
District = cd_name,
`Mean Density`,
`Median Density`
) %>%
kbl(
format = "html",
caption = "Tract‑Level Urbanization by Congressional District (Proposed 2026 Plan)"
) %>%
kable_styling(
full_width = FALSE,
bootstrap_options = c("striped", "hover", "condensed")
)
Table_2026
# ============================================================
# Step 11. EXPORT MAPS AND TABLES
# ============================================================
saveWidget(Map_2023, "Urbanization_2023_Map.html", selfcontained = TRUE)
saveWidget(Map_2026, "Urbanization_2026_Map.html", selfcontained = TRUE)
save_kable(Table_2023, "Urbanization_2023_Table.html")
save_kable(Table_2026, "Urbanization_2026_Table.html")