This code graphs modeled EPA air quality data, specifically particulate matter 2.5 microns in width (PM2.5), by Virginia census tract for the years 2016 and 2019. The modeled data is relatively course, therefor rather than assessing daily means of PM2.5, we will count the number of days each census tract exceeds the threshold of 15 μg/m3.

library(leaflet)
library(tidyquant)
library(acs)
library(tigris)
library(dplyr)

1. Import and format data

# Read in EPA data
AQ <- read.csv("DSData.csv", header=TRUE)
AQ$Date = as.Date(AQ$Date , format = "%b-%d-%Y")
head(AQ)
##         Date CensTractFips Latitude Longitude Prediction SEpred  Poll
## 1 2016-01-01   51590001200 36.54985 -79.39896    29.9335 5.2464 ozone
## 2 2016-01-01   51590000900 36.55974 -79.45373    29.9841 5.1088 ozone
## 3 2016-01-01   51590001000 36.56341 -79.41633    30.0356 5.1441 ozone
## 4 2016-01-01   51590001100 36.56701 -79.39010    29.8816 5.0633 ozone
## 5 2016-01-01   51590980100 36.56769 -79.34830    29.7589 5.0105 ozone
## 6 2016-01-01   51117930700 36.56988 -78.38239    27.5174 4.8608 ozone
# Split by pollutant
AQ.poll <-split(AQ,AQ$Poll)

# Split PM2.5 by year
pm <- AQ.poll$pm25
pm.year <- split(pm, format(as.Date(pm$Date), "%Y"))
pm.16 <-pm.year$`2016`
pm.19 <-pm.year$`2019`

# Count days over thresholds
# 2016
total.pm.16 <- pm.16 %>%
  group_by(CensTractFips) %>%
  summarize(
    n_gt15 = sum(Prediction > 15))

# 2019
total.pm.19 <- pm.19 %>%
  group_by(CensTractFips) %>%
  summarize(
    n_gt15 = sum(Prediction > 15))

total.pm.16
## # A tibble: 1,898 × 2
##    CensTractFips n_gt15
##            <dbl>  <int>
##  1   51001090100      4
##  2   51001090200      5
##  3   51001090300      4
##  4   51001090400      5
##  5   51001090500      5
##  6   51001090600      5
##  7   51001090700      5
##  8   51001090800      4
##  9   51001980100      4
## 10   51001980200      4
## # … with 1,888 more rows
total.pm.19
## # A tibble: 1,898 × 2
##    CensTractFips n_gt15
##            <dbl>  <int>
##  1   51001090100      5
##  2   51001090200      5
##  3   51001090300      6
##  4   51001090400      5
##  5   51001090500      4
##  6   51001090600      4
##  7   51001090700      4
##  8   51001090800      3
##  9   51001980100      4
## 10   51001980200      5
## # … with 1,888 more rows

2. Generate maps

# Load tracts from tigris
va.tracts <- tracts(state = 'VA', cb=TRUE, year=2015)

# 2016
# Rename to GEOID, save as character, then join AQ data with map data
names(total.pm.16)[names(total.pm.16) == 'CensTractFips'] <- 'GEOID'
total.pm.16$GEOID = as.character(total.pm.16$GEOID)
joined <- left_join(va.tracts,total.pm.16)
## Joining, by = "GEOID"
# Remove IDs with no land
joined <- joined[joined$ALAND>0,]

# Create palette, must use na.color to prevent it from showing on legend
pal <- colorNumeric(
  palette = "YlGnBu",
  domain = joined$n_gt15,
  na.color = NA)

# Pop up info for when you click on the graph
pop1 <- paste0("GEOID: ", joined$GEOID, "<br>",
                 "PM2.5 Days over 15 μg/m3 ", round(joined$n_gt15,2))

# Make leaflet map
pm25.16.map <- leaflet() %>%
  addProviderTiles("Esri.OceanBasemap") %>%
  addPolygons(data = joined, 
              fillColor = ~pal(n_gt15), 
              color = "#b2aeae", # you need to use hex colors
              fillOpacity = 0.7, 
              weight = 1, 
              smoothFactor = 0.2,
              popup = pop1) %>%
  leaflet::addLegend(pal = pal, 
                     values = joined$n_gt15, 
                     position = "bottomright", 
                     title = "Days over <br> 15 μg/m3",
                     labFormat = labelFormat(suffix = " "))

#2019
# Rename to GEOID, save as character, then join AQ data with map data
names(total.pm.19)[names(total.pm.19) == 'CensTractFips'] <- 'GEOID'
total.pm.19$GEOID = as.character(total.pm.19$GEOID)
joined2 <- left_join(va.tracts,total.pm.19)
## Joining, by = "GEOID"
# Remove tracts with no land if present
joined2 <- joined2[joined2$ALAND>0,]

# Create color palette, must use na.color to prevent it from showing on legend
pal2 <- colorNumeric(
  palette = "YlGnBu",
  domain = joined2$n_gt15,
  na.color = NA)

# Pop up info for when you click on the graph
pop2 <- paste0("GEOID: ", joined2$GEOID, "<br>",
                 "PM2.5 Days over 15 μg/m3 ", round(joined2$n_gt15,2))

# Make leaflet map
pm25.19.map <- leaflet() %>%
  addProviderTiles("Esri.OceanBasemap") %>%
  addPolygons(data = joined2, 
              fillColor = ~pal2(n_gt15), 
              color = "#b2aeae", # you need to use hex colors
              fillOpacity = 0.7, 
              weight = 1, 
              smoothFactor = 0.2,
              popup = pop2) %>%
  leaflet::addLegend(pal = pal2, 
                     values = joined2$n_gt15, 
                     position = "bottomright", 
                     title = "Days over <br> 15 μg/m3",
                     labFormat = labelFormat(suffix = " "))

3. Visualize results

pm25.16.map
pm25.19.map