Percent with a bachelor’s degree or higher
Rutherford County, 2019-2023. Source: American Community Survey

Code:

# ============================================================
# 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("gt")) install.packages("gt")
if (!require("gtExtras")) install.packages("gtExtras")
if (!require("plotly")) install.packages("plotly")

library(tidyverse)
library(tidycensus)
library(sf)
library(leaflet)
library(leaflet.extras2)
library(gt)
library(gtExtras)
library(plotly)

# ============================================================
# 1. CENSUS API KEY
# ============================================================

census_api_key("fe54531c204de591e6bd059d9a06b520289573e1")

# ============================================================
# 2. LOAD ACS CODEBOOKS
# ============================================================

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

# ============================================================
# 3. DEFINE VARIABLES OF INTEREST
# ============================================================

VariableList =
  c(
    Count_   = "DP02_0059",
    Percent_ = "DP02_0068P"
  )

# ============================================================
# 4. FETCH COUNTY SUBDIVISION DATA (TENNESSEE)
# ============================================================

mydata <- get_acs(
  geography = "county subdivision",
  state = "TN",
  variables = VariableList,
  year = 2023,
  survey = "acs5",
  output = "wide",
  geometry = TRUE
)

# ============================================================
# 5. CLEAN AND REFORMAT GEOGRAPHIC NAMES
# ============================================================

mydata <- mydata %>%
  separate_wider_delim(
    NAME,
    delim  = ", ",
    names  = c("Division", "County", "State")
  ) %>%
  mutate(County = str_remove(County, " County"))

# ============================================================
# 6. FILTER FOR A SINGLE COUNTY
# ============================================================

filtereddata <- mydata %>%
  filter(County == "Rutherford")

# ============================================================
# 7. CREATE INTERACTIVE MAP
# ============================================================

mapdata <- filtereddata %>%
  rename(
    Percent   = Percent_E,
    PercentEM = Percent_M,
    Count     = Count_E,
    CountEM   = Count_M
  ) %>%
  st_as_sf() %>%
  st_transform(4326)

# Viridis color palette
pal <- colorNumeric(
  palette = "viridis",
  domain = mapdata$Percent
)

# Build leaflet map with three base maps
DivisionMap <- leaflet(mapdata) %>%
  # Default base map (Positron Light)
  addProviderTiles("CartoDB.Positron", group = "Positron (Light)") %>%
  addProviderTiles("Esri.WorldStreetMap", group = "ESRI Streets") %>%
  addProviderTiles("Esri.WorldImagery", group = "ESRI Satellite") %>%
  addPolygons(
    fillColor = ~pal(Percent),
    fillOpacity = 0.7,
    color = "#333333",
    weight = 1,
    group = "Data Layer",
    popup = ~paste0(
      "<b>", Division, "</b><br>",
      "County: ", County, "<br>",
      "Percent: ", Percent, "%<br>",
      "Percent MOE: ±", PercentEM, "<br>",
      "Count: ", Count, "<br>",
      "Count MOE: ±", CountEM
    )
  ) %>%
  addLegend(
    pal = pal,
    values = mapdata$Percent,
    position = "topright",
    title = "Percent"
  ) %>%
  addLayersControl(
    baseGroups = c(
      "Positron (Light)",
      "ESRI Streets",
      "ESRI Satellite"
    ),
    overlayGroups = c("Data Layer", "County Borders"),
    options = layersControlOptions(
      collapsed = TRUE,
      position = "bottomleft"
    )
  )

DivisionMap

# ============================================================
# 8. INTERACTIVE PLOTLY GRAPH OF ESTIMATES WITH ERROR BARS
# ============================================================

mygraph <- plot_ly(
  data = filtereddata,
  x = ~Percent_E,
  y = ~reorder(Division, Percent_E),
  type = 'scatter',
  mode = 'markers',
  error_x = list(
    type = "data",
    array = ~Percent_M,
    visible = TRUE
  ),
  marker = list(color = "#099d91", size = 10)
) %>%
  layout(
    title = "Estimates by Area (Interactive)",
    xaxis = list(title = "2019–2023 ACS Estimate"),
    yaxis = list(title = "", automargin = TRUE)
  )

mygraph