Here is the R code for mapping and graphing household natural gas heating for the Nashville area:

# ============================================================
#  AMERICAN COMMUNITY SURVEY DATA: HOUSEHOLDS USING GAS HEAT
#  Tennessee Counties and County Subdivisions (2018–2022 ACS)
# ============================================================


# ------------------------------------------------------------
# 1. 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("mapview")) install.packages("mapview")

library(tidyverse)
library(tidycensus)
library(sf)
library(mapview)


# ------------------------------------------------------------
# 2. Transmit Census API key
# ------------------------------------------------------------

# Uncomment and paste your API key between the quote marks below:
# census_api_key("PasteYourAPIKeyBetweenTheseQuoteMarks")


# ------------------------------------------------------------
# 3. Fetch ACS codebooks (for variable lookup)
# ------------------------------------------------------------

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


# ------------------------------------------------------------
# 4. Verify target ACS variables
# ------------------------------------------------------------

ChosenVars <- filter(ProfileTables,
                     name == "DP04_0063P" |  # % households using gas heat
                       name == "DP02_0001")    # Total households

print(ChosenVars$name)
print(ChosenVars$label)
print(ChosenVars$concept)


# ------------------------------------------------------------
# 5. Specify target variables
# ------------------------------------------------------------

VariableList <- c(
  PctGasHeat_ = "DP04_0063P",
  Households_ = "DP02_0001"
)


# ============================================================
#  COUNTY-LEVEL ANALYSIS
# ============================================================


# ------------------------------------------------------------
# 6. Fetch county-level ACS data
# ------------------------------------------------------------

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


# ------------------------------------------------------------
# 7. Reformat data (split "NAME" field)
# ------------------------------------------------------------

mydata <- separate_wider_delim(
  mydata,
  NAME,
  delim = " County, ",
  names = c("County", "State")
)


# ------------------------------------------------------------
# 8. Filter for Nashville-area counties
# ------------------------------------------------------------

filtereddata <- mydata %>%
  filter(County %in% c("Davidson", "Rutherford", "Williamson",
                       "Cheatham", "Robertson", "Sumner", "Wilson"))


# ------------------------------------------------------------
# 9. Plot county-level data (dot plot)
# ------------------------------------------------------------

DotPlot <- ggplot(filtereddata, aes(x = PctGasHeat_E, y = reorder(County, PctGasHeat_E))) +
  geom_errorbarh(aes(xmin = PctGasHeat_E - PctGasHeat_M,
                     xmax = PctGasHeat_E + PctGasHeat_M)) +
  geom_point(size = 3, color = "darkblue") +
  theme_minimal(base_size = 12.5) +
  labs(
    title = "Pct. Households with Gas Heat",
    subtitle = "Nashville-area counties (2018–2022 ACS). Brackets show error margins.",
    x = "ACS Estimate (%)",
    y = ""
  )

DotPlot


# ============================================================
#  COUNTY SUBDIVISION (DISTRICT) ANALYSIS
# ============================================================


# ------------------------------------------------------------
# 10. Fetch county subdivision data
# ------------------------------------------------------------

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


# ------------------------------------------------------------
# 11. Reformat data (split "NAME" field)
# ------------------------------------------------------------

mydata <- separate_wider_delim(
  mydata,
  NAME,
  delim = ", ",
  names = c("District", "County", "State")
)


# ------------------------------------------------------------
# 12. Filter for Nashville-area counties
# ------------------------------------------------------------

filtereddata <- mydata %>%
  filter(County %in% c("Davidson County", "Rutherford County",
                       "Williamson County", "Cheatham County",
                       "Robertson County", "Sumner County",
                       "Wilson County"))


# ------------------------------------------------------------
# 13. Prepare spatial data for mapping
# ------------------------------------------------------------

mapdata <- filtereddata %>%
  rename(
    PctGasHeat = PctGasHeat_E,
    Households = Households_E
  )

mapdata <- st_as_sf(mapdata)
mapviewOptions(basemaps.color.shuffle = FALSE)


# ============================================================
#  INTERACTIVE MAPPING
# ============================================================


# ------------------------------------------------------------
# 14. Install and load additional packages for mapping
# ------------------------------------------------------------

if (!require("leaflet")) install.packages("leaflet")
if (!require("leaflet.extras2")) install.packages("leaflet.extras2")

library(leaflet)
library(leaflet.extras2)
library(leafpop)


# ------------------------------------------------------------
# 15. Define custom color palette
# ------------------------------------------------------------

mypalette <- colorRampPalette(c("lightgray", "darkorange"))


# ------------------------------------------------------------
# 16. Generate interactive map
# ------------------------------------------------------------

BetterMap <- mapview(
  mapdata,
  zcol = "PctGasHeat",
  col.regions = mypalette,
  at = seq(0, 100, 20),
  layer.name = "Pct. Households Using Gas Heat",
  popup = popupTable(
    mapdata,
    feature.id = FALSE,
    row.numbers = FALSE,
    zcol = c("County", "State", "PctGasHeat", "Households")
  )
)

BetterMap


# ============================================================
#  DATA EXPORT
# ============================================================


# ------------------------------------------------------------
# 17. Export attribute data to CSV (no geometry)
# ------------------------------------------------------------

CSVdata <- st_drop_geometry(mapdata)
write.csv(CSVdata, "mydata.csv", row.names = FALSE)