An application of R and Leaflet Mapping

Using data from Census American Community Survey (ACS) 5-year estimates, I created this interactive map to show the percentage of housing units of a given county that are burdened by housing costs >30% of their monthly income.

Here's the code:

library(tidycensus)
library(tidyverse)
library(viridis)
library(usmap)
library(RColorBrewer)
library(leaflet)
library(htmltools)
library(htmlwidgets)
library(sf)
library(tigris)
library(stringr)
#To find the variabels you're looking for, they are available here
#https://api.census.gov/data/2019/acs/acs5/profile/groups/DP04.html
mycounties <- st_read("https://catalog.data.gov/dataset/tiger-line-shapefile-2019-nation-u-s-current-county-and-equivalent-national-shapefile")
#https://catalog.data.gov/dataset/tiger-line-shapefile-2019-nation-u-s-current-county-and-equivalent-national-shapefile
#https://www.census.gov/cgi-bin/geo/shapefiles/index.php?year=2021&layergroup=Counties+%28and+equivalent%29
key <- "My Key"
data <- get_acs(geography = "county",
                variables = c("DP04_0111PE",
                              "DP04_0112PE",
                              "DP04_0113PE",
                              "DP04_0114PE",
                              "DP04_0115PE"),
                key = key,
                year = 2019,
                survey = "acs5")

# Now that we've got the raw, lets tidy it up and get the data we want ---------

widedata <- data%>%select(c(GEOID,variable,estimate)) %>% pivot_wider(names_from = variable, values_from = estimate)

colnames(widedata) <- c("GEOID","DP04_0111PE","DP04_0112PE","DP04_0113PE","DP04_0114PE","DP04_0115PE" )

#Creating our main column of interest by adding together two other columns
widedata$highburden <- widedata$DP04_0114PE+widedata$DP04_0115PE

#Just  grabbing Colorado data. 
widedata <- widedata %>% filter(substr(GEOID,1,2) == "08")

#Joining the data on a left inner join.
inter_data <- left_join(mycounties,widedata,by='GEOID')

#These are gonna let the map have hoverable labels
labels <- sprintf("%s
%g% percent burdened by housing ", 
                  inter_data$NAMELSAD,inter_data$highburden) %>% lapply(HTML)

#The colors we want the states to be depending on level
pal <- colorBin(palette = "OrRd", 9, domain = widedata$highburden)

#Creating the map!
map_interactive <- inter_data %>% 
                  st_transform(crs = "+init=epsg:4326") %>%
                  leaflet() %>%
                  addProviderTiles(provider = "CartoDB.Positron") %>%
                  addPolylines(color = "black", opacity = 1, weight = 1.5)%>%
                  addPolygons(label = labels,
                              stroke = FALSE,
                              smoothFactor = 0.3,
                              opacity = 1,
                              fillOpacity = 0.7,
                              fillColor = ~ pal(highburden),
                              highlightOptions = highlightOptions(weight = 5,
                                                                  fillOpacity = 1,
                                                                  color = "black",
                                                                  bringToFront = TRUE))%>%
                  addLegend("bottomright",
                            pal = pal,
                            values = ~ highburden,
                            title = "Percent of Units",
                            opacity = 0.7)