# ----------------------------------------------------------
# Step 1: Install & load required packages
# ----------------------------------------------------------
if (!require("tidyverse"))
install.packages("tidyverse")
if (!require("gt"))
install.packages("gt")
if (!require("leaflet"))
install.packages("leaflet")
if (!require("leafpop"))
install.packages("leafpop")
if (!require("sf"))
install.packages("sf")
if (!require("RColorBrewer"))
install.packages("RColorBrewer")
if (!require("classInt"))
install.packages("classInt")
if (!require("scales"))
install.packages("scales")
if (!require("htmlwidgets"))
install.packages("htmlwidgets")
if (!require("tidycensus"))
install.packages("tidycensus")
library(tidyverse)
library(gt)
library(sf)
library(leaflet)
library(leafpop)
library(RColorBrewer)
library(classInt)
library(scales)
library(htmlwidgets)
library(tidycensus)
# ----------------------------------------------------------
# Step 2: Area ZIP Codes
# ----------------------------------------------------------
ZIPList <- c(
"37127",
"37128",
"37129",
"37130",
"37132",
"37085",
"37118",
"37149",
"37037",
"37153",
"37167",
"37086"
)
# ----------------------------------------------------------
# Step 3: Download HUD SAFMR Excel file
# ----------------------------------------------------------
download.file(
"https://www.huduser.gov/portal/datasets/fmr/fmr2026/fy2026_safmrs.xlsx",
"rent.xlsx",
mode = "wb"
)
# ----------------------------------------------------------
# Step 4: Read Excel data
# ----------------------------------------------------------
FMR_Area <- readxl::read_xlsx(path = "rent.xlsx",
.name_repair = "universal")
# ----------------------------------------------------------
# Step 5: Filter FMR data for "ZIPList" ZIP codes,
# select and rename columns,
# and remove duplicates
# ----------------------------------------------------------
FMR_Area <- FMR_Area %>%
transmute(
ZIP = ZIP.Code,
Studio = SAFMR.0BR,
BR1 = SAFMR.1BR,
BR2 = SAFMR.2BR,
BR3 = SAFMR.3BR,
BR4 = SAFMR.4BR
) %>%
filter(ZIP %in% ZIPList) %>%
distinct()
# ----------------------------------------------------------
# Step 6: Download and unzip the ZCTA shapefile
# ----------------------------------------------------------
download.file(
"https://www2.census.gov/geo/tiger/GENZ2020/shp/cb_2020_us_zcta520_500k.zip",
"ZCTAs2020.zip",
mode = "wb"
)
unzip("ZCTAs2020.zip")
# ----------------------------------------------------------
# Step 7: Load ZCTA shapefile into R
# ----------------------------------------------------------
ZCTAMap <- read_sf("cb_2020_us_zcta520_500k.shp")
# ----------------------------------------------------------
# Step 8: Prepare ZIP column for joining
# ----------------------------------------------------------
FMR_Area$ZIP <- as.character(FMR_Area$ZIP)
# ----------------------------------------------------------
# Step 9: Join the FMR data to the ZCTA polygons
# ----------------------------------------------------------
FMR_Area_Map <- left_join(FMR_Area, ZCTAMap, by = c("ZIP" = "ZCTA5CE20"))
# ----------------------------------------------------------
# Step 10: Drop unneeded Census columns
# ----------------------------------------------------------
FMR_Area_Map <- FMR_Area_Map %>%
select(-c(AFFGEOID20, GEOID20, NAME20, LSAD20, ALAND20, AWATER20))
# ----------------------------------------------------------
# Step 11: Convert to sf object and reproject to WGS84 (EPSG:4326),
# then filter to valid geometries (fix for hover/popup)
# ----------------------------------------------------------
FMR_Area_Map <- st_as_sf(FMR_Area_Map)
if (!is.na(sf::st_crs(FMR_Area_Map))) {
FMR_Area_Map <- st_transform(FMR_Area_Map, 4326)
}
# >>> Fix: keep only rows with valid geometry so popups/labels match features <<<
FMR_Area_Map <- FMR_Area_Map %>%
dplyr::filter(!sf::st_is_empty(geometry) & !is.na(sf::st_geometry_type(geometry)))
# (Optional) See which ZIPs didn't have polygons
# missing_zips <- setdiff(ZIPList, ZCTAMap$ZCTA5CE20)
# if (length(missing_zips)) {
# message("ZIPs without ZCTA polygons (not mapped): ", paste(missing_zips, collapse = ", "))
# }
# ----------------------------------------------------------
# Step 12: Fetch ACS data
# ----------------------------------------------------------
# census_api_key("PasteYourAPIKeyBetweenTheseQuoteMarks") # <- Add your API key
Census_Data <- get_acs(
geography = "zcta",
variables = c("DP04_0047", "DP04_0045"),
year = 2024,
survey = "acs5",
output = "wide",
geometry = FALSE
)
Census_Data <- Census_Data %>%
transmute(
ZIP = GEOID,
Rentals = DP04_0047E,
Rentals_MOE = DP04_0047M,
Households = DP04_0045E,
Households_MOE = DP04_0045M
)
# ----------------------------------------------------------
# Step 13: Filter ACS data for "ZIPList" ZIP codes
# ----------------------------------------------------------
Census_Data <- Census_Data %>%
filter(ZIP %in% ZIPList)
# Left-join ACS counts & MOEs into the map data by ZIP
FMR_Area_Map <- FMR_Area_Map %>%
left_join(Census_Data, by = "ZIP")
# ==========================================================
# Step 14: EDIT MAP CONTROL POINTS AS NEEDED
# ==========================================================
ShadeBy <- "BR2" # Which numeric column to shade by: "Studio","BR1","BR2","BR3","BR4"
PaletteName <- "Blues" # Any sequential/diverging RColorBrewer palette (e.g., "Blues","OrRd","PuBuGn","GnBu")
legend_classes <- 7 # Number of legend classes/bins (typical: 4–7)
# Customizable popup labels:
# Keys MUST match columns we will pass (ZIP and the *_fmt columns created below).
# Values are the human-friendly headers shown in the popup. Reorder to change popup order.
popup_labels <- c(
ZIP = "ZIP",
Studio_fmt = "Studio",
BR1_fmt = "1-Bed",
BR2_fmt = "2-Bed",
BR3_fmt = "3-Bed",
BR4_fmt = "4-Bed",
# --- NEW FIELDS (comma-formatted strings defined below) ---
Rentals_fmt = "Renter-occupied units",
Rentals_MOE_fmt = "Renter MOE (±)",
Households_fmt = "Occupied housing units",
Households_MOE_fmt = "Occupied MOE (±)"
)
# ==========================================================
# Friendly labels for legend title (optional)
friendly_names <- c(
Studio = "Studio",
BR1 = "1-Bed",
BR2 = "2-Bed",
BR3 = "3-Bed",
BR4 = "4-Bed"
)
# Legend title shows ONLY the selected ShadeBy field (friendly name if available)
legend_title <- if (ShadeBy %in% names(friendly_names)) {
friendly_names[[ShadeBy]]
} else {
ShadeBy
}
# ----------------------------------------------------------
# Step 15: Helper: Build a Brewer palette safely for a requested size
# - Uses up to the palette's native max colors
# - Interpolates if you ask for more than the palette provides
# ----------------------------------------------------------
build_brewer_colors <- function(name, k) {
info <- RColorBrewer::brewer.pal.info
if (!name %in% rownames(info)) {
stop(sprintf("Palette '%s' not found in RColorBrewer.", name))
}
max_n <- info[name, "maxcolors"]
base <- RColorBrewer::brewer.pal(min(max_n, max(3, k)), name)
if (k <= length(base)) {
base[seq_len(k)]
} else {
grDevices::colorRampPalette(base)(k)
}
}
# ----------------------------------------------------------
# Step 16: Jenks breaks + robust fallback (prevents non-unique breaks)
# ----------------------------------------------------------
vals <- FMR_Area_Map[[ShadeBy]]
vals <- vals[!is.na(vals)]
# 1) Try Jenks
ci <- classInt::classIntervals(vals, n = legend_classes, style = "jenks")
breaks <- sort(unique(ci$brks))
# 2) If Jenks couldn't produce enough unique breaks, fall back to quantiles, then pretty
if (length(breaks) < 3) {
qbreaks <- quantile(
vals,
probs = seq(0, 1, length.out = legend_classes + 1),
na.rm = TRUE,
type = 7
)
qbreaks <- sort(unique(as.numeric(qbreaks)))
if (length(qbreaks) >= 3) {
breaks <- qbreaks
} else {
pbreaks <- pretty(range(vals, na.rm = TRUE), n = legend_classes)
pbreaks <- sort(unique(as.numeric(pbreaks)))
if (length(pbreaks) >= 3) {
breaks <- pbreaks
} else {
# As a last resort, ensure at least two unique break points around a single value
rng <- range(vals, na.rm = TRUE)
if (rng[1] == rng[2]) {
b0 <- rng[1]
eps <- if (abs(b0) < 1) 1e-9 else abs(b0) * 1e-9
breaks <- c(b0 - eps, b0 + eps)
} else {
breaks <- rng
}
}
}
}
# 3) Final guard: ensure strictly increasing, unique breaks (and at least two)
breaks <- sort(unique(breaks))
if (length(breaks) < 2) {
b0 <- vals[1]
eps <- if (abs(b0) < 1) 1e-9 else abs(b0) * 1e-9
breaks <- c(b0 - eps, b0 + eps)
}
# Palette length must match # of bins (breaks - 1)
n_bins <- max(1, length(breaks) - 1)
pal_colors <- build_brewer_colors(PaletteName, n_bins)
pal_bin <- colorBin(
palette = pal_colors,
domain = FMR_Area_Map[[ShadeBy]],
bins = breaks,
na.color = "#cccccc",
right = FALSE
)
# ----------------------------------------------------------
# Step 17: Precompute FillColor (avoid hard-coding a field name in leaflet)
# ----------------------------------------------------------
FMR_Area_Map$FillColor <- pal_bin(FMR_Area_Map[[ShadeBy]])
# Build a clean hover label with comma formatting (no dollar signs)
FMR_Area_Map$HoverLabel <- sprintf(
"ZIP %s: %s = %s",
FMR_Area_Map$ZIP,
legend_title,
ifelse(is.na(FMR_Area_Map[[ShadeBy]]), "NA", scales::comma(FMR_Area_Map[[ShadeBy]]))
)
# ----------------------------------------------------------
# Step 18: Popup table with comma formatting (no dollar signs)
# Create a formatted copy for display; keep numerics unchanged in FMR_Area_Map
# ----------------------------------------------------------
popup_data <- FMR_Area_Map %>%
mutate(
Studio_fmt = ifelse(is.na(Studio), NA, scales::comma(Studio)),
BR1_fmt = ifelse(is.na(BR1), NA, scales::comma(BR1)),
BR2_fmt = ifelse(is.na(BR2), NA, scales::comma(BR2)),
BR3_fmt = ifelse(is.na(BR3), NA, scales::comma(BR3)),
BR4_fmt = ifelse(is.na(BR4), NA, scales::comma(BR4)),
# --- NEW: formatted ACS fields for the popup ---
Rentals_fmt = ifelse(is.na(Rentals), NA, scales::comma(Rentals)),
Rentals_MOE_fmt = ifelse(is.na(Rentals_MOE), NA, scales::comma(Rentals_MOE)),
Households_fmt = ifelse(is.na(Households), NA, scales::comma(Households)),
Households_MOE_fmt = ifelse(is.na(Households_MOE), NA, scales::comma(Households_MOE))
)
# Build the popup data with user-defined labels as actual column names.
# IMPORTANT: We pass THIS object to popupTable and use its own colnames in zcol.
popup_labels <- popup_labels # keep as defined above
popup_keys <- intersect(names(popup_labels), names(popup_data)) # columns available to show
popup_display <- popup_data %>%
sf::st_drop_geometry() %>%
select(all_of(popup_keys))
colnames(popup_display) <- unname(popup_labels[popup_keys])
# ----------------------------------------------------------
# Step 19: Build the Leaflet interactive map with base layer options
# - Default: CartoDB.Positron (your original)
# - Options: Esri.WorldStreetMap and Esri.WorldImagery
# ----------------------------------------------------------
Rent_Category_Map <- leaflet(FMR_Area_Map, options = leafletOptions(preferCanvas = TRUE)) %>%
# Default/base layer (original view)
addProviderTiles(providers$CartoDB.Positron, group = "Streets (CartoDB Positron)") %>%
# Additional selectable base layers
addProviderTiles(providers$Esri.WorldStreetMap, group = "Streets (Esri World Street Map)") %>%
addProviderTiles(providers$Esri.WorldImagery, group = "Satellite (Esri World Imagery)") %>%
# Your data layer (overlay group)
addPolygons(
fillColor = ~ FillColor,
color = "#444444",
weight = 1,
opacity = 1,
fillOpacity = 0.7,
label = ~ HoverLabel,
labelOptions = labelOptions(
style = list("font-weight" = "bold"),
textsize = "12px",
direction = "auto"
),
popup = leafpop::popupTable(
popup_display,
feature.id = FALSE,
row.numbers = FALSE,
zcol = colnames(popup_display)
),
highlight = highlightOptions(
weight = 2,
color = "#000000",
fillOpacity = 0.8,
bringToFront = TRUE
),
group = "FMR by ZIP"
) %>%
# Legend
addLegend(
position = "bottomright",
pal = pal_bin,
values = FMR_Area_Map[[ShadeBy]],
title = legend_title,
opacity = 0.7,
labFormat = labelFormat(
big.mark = ",",
digits = 0,
between = " – "
)
) %>%
# Layer control to switch basemaps ONLY (overlay toggle removed)
addLayersControl(
baseGroups = c(
"Streets (CartoDB Positron)",
"Streets (Esri World Street Map)",
"Satellite (Esri World Imagery)"
),
options = layersControlOptions(collapsed = FALSE)
)
# Note: No overlayGroups parameter, so "FMR by ZIP" can't be toggled off.
# ----------------------------------------------------------
# Step 20: Display the map
# ----------------------------------------------------------
Rent_Category_Map
# ----------------------------------------------------------
# Step 21: (Optional) Save the map as an HTML file
# ----------------------------------------------------------
outfile <- paste0("FMR_",
ShadeBy,
"_",
PaletteName,
"_",
legend_classes,
"Classes.html")
htmlwidgets::saveWidget(widget = Rent_Category_Map,
file = outfile,
selfcontained = TRUE)
message("Saved map to: ", normalizePath(outfile))
# ----------------------------------------------------------
# Step 22: Downshifting map file to a data file &
# showing contents as a table
# ----------------------------------------------------------
Data_From_Map <- st_drop_geometry(FMR_Area_Map) %>%
select(-c(FillColor,HoverLabel))
Data_From_Map_Table <- gt(Data_From_Map) %>%
tab_header(title = "Data") %>%
cols_align(align = "left")
Data_From_Map_Table