The R code shown below will generate each of the maps shown above.
# Installing and loading 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")
if (!require("leaflet")) install.packages("leaflet")
if (!require("leaflet.extras2")) install.packages("leaflet.extras2")
library(tidyverse)
library(tidycensus)
library(sf)
library(mapview)
library(leaflet)
library(leafpop)
# Transmitting API key obtained from:
# https://api.census.gov/data/key_signup.html
census_api_key("PasteYourAPIKeyBetweenTheseQuoteMarks")
# Fetching 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)
# Specifying target variables
VariableList =
c(Estimate_ = "DP02_0114P",
Households_ = "DP02_0001",
Population_ = "DP02_0124")
# Some potentially interesting ACS variables
# DP02_0124 - Total population
# DP05_0037P - Pct. white alone
# DP03_0074P - Pct. HH with SNAP benefits in past 12 months
# DP02_0001 - Total households
# DP03_0072P - Pct. with cash public assistance
# DP03_0128P - Pct. all people below poverty level
# DP03_0099P - Pct. no health insurance coverage
# DP04_0142P - GRAPI 35% or more
# DP04_0115P - Sel. Mo. Owner Costs 35% or more
# DP03_0092 - Median earnings for workers
# DP03_0093 - Median earnings, full-time year-round male workers
# DP03_0094 - Median earnings, full-time year-round female workers
# DP03_0062 - Median household income
# DP02_0068P - Pct. 25 or older with a bachelor's degree or higher
# DP02_0114P - Pct. 5 and older non-English spoken at home
# DP05_0076P - Pct. Hispanic, of any race
# ---------- County subdivision analysis ----------
# Fetching data
mydata <- get_acs(
geography = "county subdivision",
state = "TN",
variables = VariableList,
year = 2023,
survey = "acs5",
output = "wide",
geometry = TRUE)
# Reformatting data
mydata <-
separate_wider_delim(mydata,
NAME,
delim = ", ",
names = c("Division", "County", "State"))
# Filtering data
filtereddata <- mydata %>%
filter(County %in% c("Rutherford County"))
# Plotting data
mygraph <- ggplot(filtereddata, aes(x = Estimate_E, y = reorder(Division, Estimate_E))) +
geom_errorbarh(aes(xmin = Estimate_E - Estimate_M, xmax = Estimate_E + Estimate_M)) +
geom_point(size = 3, color = "#099d91") +
theme_minimal(base_size = 12.5) +
theme(axis.text.y = element_text(size = 8)) +
labs(title = "Estimates by area",
subtitle = "County subdivisions. Brackets show error margins.",
x = "2019-2023 ACS estimate",
y = "")
mygraph
# Mapping data
mapdata <- filtereddata %>%
rename(Estimate = Estimate_E,
Range = Estimate_M,
Population = Population_E,
Households = Households_E)
mapdata <- st_as_sf(mapdata)
mapviewOptions(basemaps.color.shuffle = FALSE)
DivisionMap <- mapview(mapdata,
zcol = "Estimate",
layer.name = "Estimate",
popup = popupTable(
mapdata,
feature.id = FALSE,
row.numbers = FALSE,
zcol = c(
"State",
"County",
"Division",
"Estimate",
"Range",
"Population",
"Households")))
DivisionMap
# Filtering for Nashville-area counties
filtereddata2 <- mydata %>%
filter(County %in% c("Rutherford County",
"Davidson County",
"Williamson County"))
# Mapping multi-county data
mapdata2 <- filtereddata2 %>%
rename(Estimate = Estimate_E,
Range = Estimate_M,
Population = Population_E,
Households = Households_E)
mapdata2 <- st_as_sf(mapdata2)
# Fetching county borders
Counties <- get_acs(
geography = "county",
state = c("TN"),
county = c("Rutherford", "Davidson", "Williamson"),
variables = VariableList,
year = 2023,
survey = "acs5",
output = "wide",
geometry = TRUE)
mapviewOptions(basemaps.color.shuffle = FALSE)
DivisionMap2 <- mapview(
Counties,
zcol = NULL,
legend = FALSE,
col.regions = "transparent",
alpha.regions = 0,
lwd = 1.5,
color = "black"
) + mapview(
mapdata2,
zcol = "Estimate",
layer.name = "Estimate",
popup = popupTable(
mapdata2,
feature.id = FALSE,
row.numbers = FALSE,
zcol = c(
"State",
"County",
"Division",
"Estimate",
"Range",
"Population",
"Households"
)
)
)
DivisionMap2
# Mapping the Whole state
filtereddata3 <- mydata
# Mapping multi-county data
mapdata3 <- filtereddata3 %>%
rename(Estimate = Estimate_E,
Range = Estimate_M,
Population = Population_E,
Households = Households_E)
mapdata3 <- st_as_sf(mapdata3)
mapviewOptions(basemaps.color.shuffle = FALSE)
DivisionMap3 <- mapview(mapdata3,
zcol = "Estimate",
layer.name = "Estimate",
popup = popupTable(
mapdata3,
feature.id = FALSE,
row.numbers = FALSE,
zcol = c(
"State",
"County",
"Division",
"Estimate",
"Range",
"Population",
"Households")))
DivisionMap3