Graphics and maps

The interactive map.

The plot graphic.

Code

This is the code I used.

############################################################# STEP 0: Save your API key

# The script will not work unless you first paste your
# API key into the API key setup code, shown above, and
# run the code. The code needs to run only once per computer.

############################################################

############################################################# STEP 1: Load Packages
############################################################

if (!require("tidycensus")) install.packages("tidycensus")
if (!require("tidyverse")) install.packages("tidyverse")
if (!require("kableExtra")) install.packages("kableExtra")
if (!require("sf")) install.packages("sf")
if (!require("leaflet")) install.packages("leaflet")
if (!require("htmlwidgets")) install.packages("htmlwidgets")
if (!require("plotly")) install.packages("plotly")

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

############################################################
# STEP 2: Verify Census API Key
############################################################

if (!nzchar(Sys.getenv("CENSUS_API_KEY"))) {
  stop(
    "No Census API key was found.\n",
    "Restart RStudio and try again."
  )
}

census_key <- Sys.getenv("CENSUS_API_KEY")

############################################################
# STEP 3: Configure Your Analysis
############################################################

# ACS year

acs_year <- 2024

# State abbreviation

state_abbrev <- "TN"

# ACS variable
#
# Examples:
#
# DP05_0001  = Total population
# DP03_0062  = Median Household Income
# DP03_0128P = Poverty Rate
# DP02_0068P = Bachelor's Degree or Higher Rate
# DP05_0018  = Median Age
# DP04_0046P = Homeownership Rate
# DP04_0134  = Median Gross Rent
# DP04_0142P = Rent Burden Rate
# DP05_0090P = Hispanic or Latino Population Percentage
# DP02_0094P = Foreign-Born Percentage
# DP03_0009P = Unemployment Rate
# DP03_0025  = Average Commute Time

my_variable <- "DP03_0062"

# Geography level
#
# Recommended ACS geography options:
#
# "county"
# "county subdivision"
# "tract"

geography_level <- "county"

# County filter
#
# Use NULL to include all counties in the state.
#
# Examples:
#
# selected_counties <- NULL
# selected_counties <- "Rutherford"
# selected_counties <- c(
#   "Davidson",
#   "Rutherford",
#   "Williamson",
#   "Cheatham",
#   "Robertson",
#   "Sumner",
#   "Wilson"
# )

selected_counties <- c(
  "Davidson",
  "Rutherford",
  "Williamson",
  "Cheatham",
  "Robertson",
  "Sumner",
  "Wilson"
)

############################################################
# STEP 4: Explore ACS Codebooks (Optional)
############################################################

ProfileTables <- load_variables(
  acs_year,
  "acs5/profile"
)

# DetailedTables <- load_variables(
#   acs_year,
#   "acs5"
# )

# SubjectTables <- load_variables(
#   acs_year,
#   "acs5/subject"
# )

############################################################
# STEP 5: Download ACS Data
############################################################

if (is.null(selected_counties)) {
  
  mydata <- get_acs(
    geography = geography_level,
    state = state_abbrev,
    year = acs_year,
    survey = "acs5",
    variables = my_variable,
    geometry = TRUE
  )
  
} else {
  
  mydata <- get_acs(
    geography = geography_level,
    state = state_abbrev,
    county = selected_counties,
    year = acs_year,
    survey = "acs5",
    variables = my_variable,
    geometry = TRUE
  )
  
}

############################################################
# STEP 6: Keep Needed Variables and Geometry
############################################################

mydata <- mydata |>
  select(
    Area = NAME,
    Variable = variable,
    Estimate = estimate,
    Margin_of_Error = moe,
    geometry
  )

############################################################
# STEP 7: Transform Coordinates for Leaflet
############################################################

mydata <- st_transform(
  mydata,
  crs = 4326
)

############################################################
# STEP 8: Download County Boundaries (Optional Overlay)
############################################################

if (is.null(selected_counties)) {
  
  county_boundaries <- get_acs(
    geography = "county",
    state = state_abbrev,
    year = acs_year,
    survey = "acs5",
    variables = my_variable,
    geometry = TRUE
  )
  
} else {
  
  county_boundaries <- get_acs(
    geography = "county",
    state = state_abbrev,
    county = selected_counties,
    year = acs_year,
    survey = "acs5",
    variables = my_variable,
    geometry = TRUE
  )
  
}

county_boundaries <- st_transform(
  county_boundaries,
  4326
)

############################################################
# STEP 9: Sort Results from Highest to Lowest
############################################################

mydata <- mydata |>
  arrange(desc(Estimate))

############################################################
# STEP 10: Display Results as a Table
############################################################

ResultsTable <- mydata |>
  st_drop_geometry() |>
  kbl(
    col.names = c(
      "Area",
      "Variable",
      "Estimate",
      "Margin of Error"
    ),
    format.args = list(big.mark = ",")
  ) |>
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = FALSE
  )

############################################################
# STEP 11: Create Pop-Up Content
############################################################

mydata$popup <- paste0(
  mydata$Area,
  "\n\nEstimate: ",
  format(mydata$Estimate, big.mark = ","),
  "\nPlus/Minus: ",
  format(mydata$Margin_of_Error, big.mark = ",")
)

############################################################
# STEP 12: Build Color Palette
############################################################

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

pal <- colorBin(
  palette = "Blues",
  domain = mydata$Estimate,
  bins = qs,
  pretty = FALSE
)

############################################################
# STEP 13: Create Choropleth Map
############################################################

AreaMap <- leaflet(mydata) |>
  addProviderTiles(
    providers$Esri.WorldTopoMap
  ) |>
  addPolygons(
    fillColor = ~pal(Estimate),
    fillOpacity = 0.5,
    color = "black",
    weight = 1,
    popup = ~popup
  )

# Optional county boundary overlay for
# sub-county geographies

if (geography_level != "county") {
  
  AreaMap <- AreaMap |>
    addPolygons(
      data = county_boundaries,
      fill = FALSE,
      color = "black",
      weight = 2,
      opacity = 1,
      group = "County Boundaries"
    )
  
}

AreaMap <- AreaMap |>
  addLegend(
    pal = pal,
    values = ~Estimate,
    title = "Estimate",
    labFormat = labelFormat(
      big.mark = ","
    )
  )

############################################################
# STEP 14: Create Plotly Dot Plot
############################################################

plotdf <- mydata |>
  st_drop_geometry() |>
  mutate(
    point_color = pal(Estimate),
    
    y_ordered = reorder(
      Area,
      Estimate
    ),
    
    hover_text = paste0(
      Area,
      "\nEstimate: ",
      format(Estimate, big.mark = ","),
      "\nPlus/Minus: ",
      format(Margin_of_Error, big.mark = ",")
    )
  )

DotPlot <- plot_ly(
  data = plotdf,
  x = ~Estimate,
  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 = ~Margin_of_Error,
    arrayminus = ~Margin_of_Error,
    color = "rgba(0,0,0,0.65)",
    thickness = 1
  ),
  text = ~hover_text,
  hovertemplate = "%{text}"
) |>
  layout(
    title = list(
      text = paste0(
        "ACS Estimates by Geography",
        "\nError bars show ACS margins of error."
      )
    ),
    xaxis = list(
      title = "Estimate",
      tickformat = ",.0f",
      automargin = TRUE
    ),
    yaxis = list(
      title = "",
      automargin = TRUE,
      categoryorder = "array",
      categoryarray = levels(plotdf$y_ordered)
    ),
    margin = list(
      l = 200,
      r = 20,
      b = 60,
      t = 60,
      pad = 2
    )
  )

############################################################
# STEP 15: Export Dot Plot (Optional)
############################################################

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

############################################################
# STEP 16: Export Map (Optional)
############################################################

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

############################################################
# STEP 17: Show output
############################################################

AreaMap # Name of the map
DotPlot # Name of the dot plot with error bars