Figure 1: Estimated proportions of Nashville-area renter households with GRAPIs of 35 percent or more, 2020-2024

Code:

# ----------------------------------------------------------
# Step 1: Install required packages (if missing)
# ----------------------------------------------------------

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("htmlwidgets"))
  install.packages("htmlwidgets")
if (!require("plotly"))
  install.packages("plotly")   # For the interactive dot plot

# ----------------------------------------------------------
# Step 2: Load libraries
# ----------------------------------------------------------

library(tidyverse)
library(tidycensus)
library(sf)
library(leaflet)
library(htmlwidgets)
library(plotly)

# ----------------------------------------------------------
# Step 3: Transmit Census API key
# ----------------------------------------------------------

census_api_key("PasteYourAPIKeyBetweenTheseQuoteMarks")

# ----------------------------------------------------------
# Step 4: Fetch ACS codebooks (for variable lookup if needed)
# ----------------------------------------------------------

# DetailedTables <- load_variables(2024, "acs5", cache = TRUE)
# SubjectTables  <- load_variables(2024, "acs5/subject", cache = TRUE)
ProfileTables  <- load_variables(2024, "acs5/profile", cache = TRUE)

# DP04_0142P: 
# Percent!!GROSS RENT AS A PERCENTAGE OF HOUSEHOLD INCOME (GRAPI)
# !!Occupied units paying rent (excluding units where GRAPI cannot
# be computed)!!35.0 percent or more

# ----------------------------------------------------------
# Step 5: Specify target variable(s)
# ----------------------------------------------------------

VariableList <- c(Estimate_ = "DP04_0142P")

# ----------------------------------------------------------
# Step 6: Fetch ACS data (PUMAs, Tennessee)
# ----------------------------------------------------------

mydata <- get_acs(
  geography = "public use microdata area",
  state = "TN",
  variables = VariableList,
  year = 2024,
  survey = "acs5",
  output = "wide",
  geometry = TRUE
)

# ----------------------------------------------------------
# Step 7: Rename the "NAME" field as "Area"
# ----------------------------------------------------------

mydata <- mydata %>% 
  rename(Area = NAME)

# ----------------------------------------------------------
# Step 8: Filter flexibly, based on the "Area" field 
# ----------------------------------------------------------

search_terms <- c("Rutherford",
                  "Davidson",
                  "Williamson",
                  "Murfreesboro")

# Note: use c("ALL") to capture all areas

if (length(search_terms) == 1 && toupper(search_terms) == "ALL") {
  filtereddata <- mydata
} else {
  or_pattern <- paste(
    stringr::str_replace_all(search_terms, "([\\W])", "\\\\\\1"), 
    collapse = "|"
  )
  
  filtereddata <- mydata %>%
    dplyr::filter(
      stringr::str_detect(
        Area %>% tidyr::replace_na(""),
        stringr::regex(or_pattern, ignore_case = TRUE)
      )
    )
}

# ----------------------------------------------------------
# Step 9: Prepare data for mapping (rename, as sf, CRS)
# ----------------------------------------------------------

mapdata <- filtereddata %>%
  rename(
    Estimate = Estimate_E,
    Range = Estimate_M
  ) %>%
  st_as_sf()

# Ensure CRS is WGS84 for Leaflet
mapdata <- st_transform(mapdata, 4326)

# ----------------------------------------------------------
# Step 10: Build color palette with quantile-based breaks
# ----------------------------------------------------------

qs <- quantile(mapdata$Estimate, probs = seq(0, 1, length.out = 6), na.rm = TRUE)

pal <- colorBin(
  palette = "Blues", # Can specify other palettes here
  domain = mapdata$Estimate,
  bins = qs,
  pretty = FALSE
)

# ----------------------------------------------------------
# Step 11: Build the plotly dot plot with error bars
# ----------------------------------------------------------

plotdf <- filtereddata %>%
  sf::st_drop_geometry() %>% 
  dplyr::mutate(
    point_color = pal(Estimate_E),
    y_ordered   = reorder(Area, Estimate_E),            
    hover_text  = paste0("Area: ", Area)
  )

mygraph <- plotly::plot_ly(
  data = plotdf,
  x = ~Estimate_E,
  y = ~as.character(y_ordered),   
  type = "scatter",
  mode = "markers",
  showlegend = FALSE,             
  marker = list(
    color = ~point_color,
    size  = 8,
    line  = list(color = "rgba(120,120,120,0.9)", width = 0.5)
  ),
  error_x = list(
    type       = "data",
    array      = ~Estimate_M,      
    arrayminus = ~Estimate_M,      
    color      = "rgba(0,0,0,0.65)",
    thickness  = 1
  ),
  text = ~hover_text,
  hovertemplate = "%{text}<br>%{x:,}<extra></extra>"
) %>%
  plotly::layout(
    title = list(
      text = "Estimates by area<br><sup>County subdivisions. Brackets show error margins.</sup>"
    ),
    xaxis = list(
      title      = "ACS estimate",
      automargin = TRUE,
      tickformat = ",.0f"   
    ),
    yaxis = list(
      title         = "",
      automargin    = TRUE,
      categoryorder = "array",
      # preserve the reorder() order on the y-axis
      categoryarray = levels(plotdf$y_ordered)
    ),
    margin = list(l = 180, r = 20, b = 60, t = 60, pad = 2)
  )

mygraph

# ----------------------------------------------------------
# Step 12: Create popup content for the map
# ----------------------------------------------------------

mapdata$popup <- paste0(
  "<strong>", mapdata$Area, "</strong><br/>",
  "<hr>",
  "Estimate: ", format(mapdata$Estimate, big.mark = ","), "<br/>",
  "Plus/Minus: ", format(mapdata$Range, big.mark = ",")
)

# ----------------------------------------------------------
# Step 13: Build the Leaflet map
# ----------------------------------------------------------

DivisionMap <- leaflet(mapdata) %>%
  # Choose one basemap:
  addProviderTiles(providers$CartoDB.Positron) %>%
  # addProviderTiles(providers$Esri.WorldStreetMap, group = "Streets (Esri World Street Map)") %>%
  # addProviderTiles(providers$Esri.WorldImagery,   group = "Satellite (Esri World Imagery)") %>%
  addPolygons(
    fillColor   = ~pal(Estimate),
    fillOpacity = 0.5, 
    color       = "black",
    weight      = 1,
    popup       = ~popup
  ) %>%
  addLegend(
    pal    = pal,
    values = ~Estimate,
    title  = "Estimate",
    labFormat = labelFormat(big.mark = ",")
  )

DivisionMap

# ----------------------------------------------------------
# Step 14: Export graph as a standalone HTML file
# ----------------------------------------------------------
# This creates a fully self-contained HTML file for the dot plot.

saveWidget(
  widget = as_widget(mygraph),
  file = "ACSGraph.html",
  selfcontained = TRUE
)

# ----------------------------------------------------------
# Step 15: Export map as a standalone HTML file
# ----------------------------------------------------------
# This creates a fully self-contained HTML you can open or share.
# Adjust the path/filename as you like.

saveWidget(
  widget = DivisionMap,
  file = "ACSMap.html",
  selfcontained = TRUE
)