I recently investigated a subset of Social Vulnerability Index (SVI) data from CDC/ATSDR for Wake County, NC. I’m working on building my mapping skills in R, so this project felt like a great opportunity to play around with the tigris and sf packages.
Load necessary packages to use functions throughout the document.
library(tidyverse) # data wrangling
library(janitor) # clean up
library(labelled) # variable labels
library(scales) # scaling functions for visualization
library(openxlsx) # xslx reading
library(expss) # labels
library(sf) # mapping
library(tidycensus) # for using Census data
library(tigris) # download TIGER/Line shapefiles from Census
library(ggmap) # spatial viz
library(ggthemes) # ggplot themes
library(GGally) # ggplot2 extension
library(grateful) # package citations
library(knitr)After importing the file, I filtered the data to Wake County, NC for the purposes of this analysis. Using the tigris package, I created a shapefile tibble with Wake County’s census tracts. I also incorporated supplemental data by importing a file with the health service zones in Wake County, which allow me to further spatially aggregate the data.
## SVI file
wake_tbl <- read.xlsx('Data/svi-data.xlsx') %>%
as_tibble() %>% # Convert to tibble
clean_names() %>% # Clean variable names
filter(county == "Wake County") # Filter to Wake County only
## Wake census tracts
wake_tracts <- tracts("NC", "Wake")
## Health service zones
wake_service_zones <- read.xlsx('Data/health_service_zones.xlsx') %>%
as_tibble() %>% # Convert to tibble
clean_names() %>% # Clean variable names
mutate(across(geoid:tract, ~as.character(.x))) # Convert variables to characterSteps to clean the data included: (1) converting the indicators provided to numeric, (2) creating overall variables for Wake County’s total population and housing units, (3) cleaning the location variable to only include the census tract, (4) creating a population density variable, (5) converting indicators that are counts to percents for more consistent comparison across tracts, (6) calculating a count of overcrowded housing units since it was the only indicator provided as a percent, and (7) multiplying percentages by 100 and round to 1 decimal for easier interpretation in figures to be developed.
wake_tbl_clean <- wake_tbl %>%
mutate(
# Convert variables to numeric
across(area_in_sq_mi:people_25_w_o_high_school_diploma, ~as.numeric(.x)),
# Create variables for total county population and total housing units
total_pop = sum(population),
total_housing_units = sum(housing_units),
# Clean location variable
location = str_extract(location, "^[^;]+"),
# Create variable for population density
pop_density = population/area_in_sq_mi,
# Update indicators to percents for comparison
## Denominator = Population
perc_below_150_percent_poverty = ppl_below_150_percent_poverty/population,
perc_bipoc_residents = bipoc_residents/population,
perc_w_o_high_school_diploma = people_25_w_o_high_school_diploma/population,
## Denominator = Households
perc_hh_no_vehicle = households_with_no_vehicle/households,
# Calculate count of overcrowded housing units
overcrowded_housing_units = percent_of_overcrowded_housing_units*housing_units,
# Multiply percentages by 100 and round to 1 decimal
across(starts_with("perc_"), ~round((.x*100), digits = 1))
)I merged the cleaned dataset and the wake_tracts shapefile, joining by the “GEOID” field in the shapefile and the “fips_code” field in the cleaned data. Then, I merged the resulting file with the wake_service_zones shapefile to incorporate service zone data.
wake_sf <- left_join(
wake_tracts,
wake_tbl_clean,
by = c("GEOID" = "fips_code")) %>%
drop_na(starts_with("perc"))
wake_service_zones_sf <- left_join(
wake_sf,
wake_service_zones,
by = c("GEOID" = "geoid")) %>%
drop_na(starts_with("perc")) %>%
mutate(zone_name = ifelse(zone_name=="East Central", "East\nCentral", zone_name))The most socially vulnerable areas are census tracts in the top 25% for rates of no high school diploma (age 25+), individuals below 150% FPL, overcrowded housing, households without a vehicle, and BIPOC residents. The least socially vulnerable areas are in the bottom 25% for these indicators.
# Create layered map with service zones, tracts, and most/least socially vulnerable areas
## Create base map
base_map <- ggplot(data = wake_service_zones_sf) +
geom_sf(aes(geometry = geometry), fill = "white", color = "gray", size = 0.2) +
theme_map()
## Add the most and least vulnerable tracts layer
map_top_tracts <- base_map +
geom_sf(data = highest_quartile_tracts_geo, aes(fill = "Most Vulnerable"), color = "gray", alpha = .75) +
geom_sf(data = lowest_quartile_tracts_geo, aes(fill = "Least Vulnerable"), color = "gray", alpha = .75) +
scale_fill_manual(
name = "Tract-Level Social Vulnerability",
values = c("Least Vulnerable" = "#2c7fb8", "Most Vulnerable" = "#edf8b1"),
labels = c("Least Vulnerable Areas", "Most Vulnerable Areas")
)
## Add the service zones layer
zones <- wake_service_zones_sf %>%
group_by(zone_name) %>%
summarize(geometry = st_union(geometry))
map_with_zones <- map_top_tracts +
geom_sf(data = zones, fill = NA, color = "black", size = 20)
## Add the service zone label layer
zone_centroids <- zones %>%
st_centroid()
map_with_labels <- map_with_zones +
geom_sf_text(data = zone_centroids, aes(label = zone_name), size = 3, color = "black", fontface = "bold") +
labs(
title = "Spatial Disparities in Social Vulnerability Factors in Wake County",
caption = caption) +
theme(
legend.position = "right",
legend.direction = "vertical",
legend.title = element_text(size = 10, face = "bold"),
legend.text = element_text(size = 10),
plot.title = element_text(size = 12, face = "bold"),
plot.caption = element_text(size = 10, face = "italic", hjust=0, vjust=1))
map_with_labels
Wake County is divided into 8 service zones. This map shows that the most socially vulnerable communities are clustered together, concentrated in the East Central service zone and nearby tracts on the borders of the North Central, South Central, and West Central service zones. On the other hand, the least socially vulnerable communities are more distinctly spread across the county in the Northern, North Central, West, and South Central service zones. We can use this information to help inform which tracts and service zones county-level investments should go to.
Established by Congress in 2017 through the Tax Cuts & Jobs Act, Opportunity Zones are a nationwide incentive program to encourage long-term investment and economic growth in communities that have experienced historical disinvestment. The Opportunity Zones Program currently defines eligible census tracts as areas where the poverty rate is 20% or higher and/or family income is less than 80% of the area’s median income. In Raleigh, there are 10 census tracts designated as Opportunity Zones.
The shaded areas in this map are the most socially vulnerable areas, or census tracts in the top 25% for rates of no high school diploma (age 25+), individuals below 150% FPL, overcrowded housing, households without a vehicle, and BIPOC residents. Among these tracts, green areas are designated Opportunity Zones while gray areas are not.
highest_quartile_tracts_geo_oz <- highest_quartile_tracts_geo %>%
mutate(oz = ifelse(
NAME %in% c(
'540.18','527.04','540.08','527.06','524.07','509', '545', '530.09', '528.03', '528.07'
), "1", "0"))
# Create layered map with service zones, tracts, and most/least socially vulnerable areas
## Create base map
base_map <- wake_service_zones_sf %>%
filter(zone_name %in% c('East\nCentral', 'North Central', 'West Central', 'South Central')) %>%
ggplot() +
geom_sf(aes(geometry = geometry), fill = "white", color = "gray", size = 0.2) +
theme_map()
## Add the most vulnerable tracts layer with Opportunity Zones highlighted
map_oz <- base_map +
geom_sf(data = highest_quartile_tracts_geo_oz, aes(fill = oz), color = "black", alpha = .75) +
scale_fill_manual(
name = "Opportunity Zone Status",
values = c("1" = "#2ca25f", "0" = "gray"),
labels = c("Not Opportunity Zone", "Opportunity Zone")
)
## Add the service zones layer
zones <- wake_service_zones_sf %>%
filter(zone_name %in% c('East\nCentral', 'North Central', 'West Central', 'South Central', 'East Central')) %>%
group_by(zone_name) %>%
summarize(geometry = st_union(geometry))
map_with_zones <- map_oz +
geom_sf(data = zones, fill = NA, color = "black", size = 20)
## Add the service zone label layer
zone_centroids <- zones %>%
st_centroid()
map_with_labels <- map_with_zones +
geom_sf_text(data = zone_centroids, aes(label = zone_name), size = 3, color = "black", fontface = "bold") +
labs(
title = "Opportunity Zones Among Top Socially Vulnerable Areas in Wake County",
caption = "CDC/ATSDR Social Vulnerability Index, City of Raleigh Opportunity Zones") +
theme(
legend.position = "right",
legend.direction = "vertical",
legend.title = element_text(size = 10, face = "bold"),
legend.text = element_text(size = 10),
plot.title = element_text(size = 12, face = "bold"),
plot.caption = element_text(size = 10, face = "italic", hjust=0, vjust=1))
map_with_labels
Of the 10 tracts identified as the most socially vulnerable areas, 4 are designated Opportunity Zones by the city of Raleigh. This overlap suggests the initiative is targeting areas with high social vulnerability, but could do better. Expanding Opportunity Zone tract eligibility beyond poverty and income measures could lead to improved social vulnerability among communities that have been historically marginalized. Policymakers should consider additional investment incentive programs for communities with lower rates of educational attainment and higher rates of severe housing problems, for example.
# Add packages used to bibliography
cite_packages(output = "table", out.dir = ".") %>% kable()| Package | Version | Citation |
|---|---|---|
| base | 4.4.2 | R Core Team (2024) |
| expss | 0.11.6 | Demin (2023) |
| GGally | 2.2.1 | Schloerke et al. (2024) |
| ggmap | 4.0.0 | Kahle and Wickham (2013) |
| ggthemes | 5.1.0 | Arnold (2024) |
| janitor | 2.2.0 | Firke (2023) |
| knitr | 1.48 | Xie (2014); Xie (2015); Xie (2024) |
| labelled | 2.13.0 | Larmarange (2024) |
| openxlsx | 4.2.5.2 | Schauberger and Walker (2023) |
| rmarkdown | 2.27 | Xie, Allaire, and Grolemund (2018); Xie, Dervieux, and Riederer (2020); Allaire et al. (2024) |
| scales | 1.3.0 | Wickham, Pedersen, and Seidel (2023) |
| sf | 1.0.19 | Pebesma (2018); Pebesma and Bivand (2023) |
| tidycensus | 1.6.7 | Walker and Herman (2024) |
| tidyverse | 2.0.0 | Wickham et al. (2019) |
| tigris | 2.1 | Walker (2024) |
#knitr::kable(pkgs)