The figure and map below show 2022 five-year American Community Survey estimates of the proportion of housholds with gas furnaces in the Nashville area.

Percentages of households with gas furnaces, by county subdivision

The code for the data extraction, graphic and map:

# Installing and loading required packages

if (!require("tidyverse")) install.packages("tidyverse")
if (!require("tidycensus")) install.packages("tidyverse")
if (!require("sf")) install.packages("tidyverse")
if (!require("mapview")) install.packages("tidyverse")

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

# Transmitting API key

#census_api_key("PasteYourAPIKeyBetweenTheseQuoteMarks")

# Fetching ACS codebooks

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

# Double checking target variables

ChosenVars <- filter(ProfileTables,name == "DP04_0063P"|
                       name == "DP02_0001")
print(ChosenVars$name)
print(ChosenVars$label)
print(ChosenVars$concept)

# Specifying target variables

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

# Fetching data

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

# Reformatting data

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

# Filtering data

filtereddata <- mydata %>% 
  filter(County == "Davidson"|
           County == "Rutherford"|
           County == "Williamson"|
           County == "Cheatham"|
           County == "Robertson"|
           County == "Sumner"|
           County == "Wilson")

# Plotting data

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. Brackets show error margins.", 
       x = "2018-2022 ACS estimate", 
       y = "")

DotPlot

# Mapping data

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

mapdata <- st_as_sf(mapdata)

mapviewOptions(basemaps.color.shuffle = FALSE)

# Installing and activating some additional packages

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

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

# Defining a custome color palette for the map

mypalette = colorRampPalette(c('lightgray', 'darkorange'))

# Generating the map

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

BetterMap

# Exporting data in .csv format

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