This was a one-off assignment given to me by my work at the ARC. They wanted to check air quality relative to data center locations. Atlanta is tricky when it comes to air quality monitoring because there are only a few (3 I believe) in the metro area.

Libraries Used

library(tidycensus)
library(tidyverse)
library(tmap)
library(sf)
library(dplyr)
library(stringr)
library(openair)
library(RAQSAPI)
library(tigris)
options(tigris_use_cache = TRUE)

Gathering Data

quarterly_data <- aqs_quarterlysummary_by_county(
  parameter = c("88101"),
  bdate = as.Date("20250101", format = "%Y%m%d"),
  edate = as.Date("20251231", format = "%Y%m%d"),
  stateFIPS = "13",
  countycode = c('135','247','121','067','057',
                 '117','063','089','113','151','097'))

ga_counties <- counties(state = "13", cb = TRUE, class = "sf")

aq_clean <- quarterly_data %>%
  mutate(GEOID = paste0(state_code, county_code)) %>%
  group_by(GEOID) %>%
  summarize(avg_pm25 = mean(arithmetic_mean, na.rm = TRUE))

map_data <- ga_counties %>%
  left_join(aq_clean, by = "GEOID")

Building Initial Map

ggplot(map_data) +
  geom_sf(aes(fill = avg_pm25), color = "white", size = 0.2) +
  scale_fill_viridis_c(
    option = "magma", 
    name = "PM2.5 (µg/m³)",
    na.value = "grey90") +
  labs(
    title = "Georgia Air Quality Choropleth",
    subtitle = "Annual Mean PM2.5 Concentration by County",
    caption = "Source: EPA AQS Data Mart"
  ) +
  theme_void()

Examining Data Center Location and EJ Scores with Air Quality

EJ2024 <- st_read("C:/Users/xavier/OneDrive - Atlanta Regional Commission/Desktop/ARP/Data/ATLEJ/ATLEJ2024.shp")
## Reading layer `ATLEJ2024' from data source 
##   `C:\Users\xavier\OneDrive - Atlanta Regional Commission\Desktop\ARP\Data\ATLEJ\ATLEJ2024.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2796 features and 5 fields (with 5 geometries empty)
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -85.60516 ymin: 30.35785 xmax: -80.84038 ymax: 35.00124
## Geodetic CRS:  NAD83
DC <- st_read('C:/Users/xavier/OneDrive - Atlanta Regional Commission/Desktop/ARP/Data/AllDC/all data centers.shp')
## Reading layer `all data centers' from data source 
##   `C:\Users\xavier\OneDrive - Atlanta Regional Commission\Desktop\ARP\Data\AllDC\all data centers.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 77 features and 17 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -84.9617 ymin: 33.38011 xmax: -83.78302 ymax: 34.31076
## Geodetic CRS:  WGS 84
DC <- DC %>% 
  st_transform(st_crs(EJ2024))
map_data <- map_data %>%
  st_transform(st_crs(EJ2024))

tmap_mode('view')

# Layer 1: EJ2024 Polygons#
map1 <- tm_shape(EJ2024) +
  tm_polygons(
    fill = "ej_score", 
    fill.scale = tm_scale_intervals(values = "greys"),
    col = "black", 
    lwd = 0.2,
    fill_alpha = 1) +
  
  # Layer 2: Air Quality Data
  tm_shape(map_data) +
  tm_fill(
    fill = "avg_pm25", 
    fill.scale = tm_scale_intervals(values = "Reds"),
    fill_alpha = 0.5) +
  
  # Layer 3: Data Centers
  tm_shape(DC) +
  tm_dots(
    fill = "status",
    size = 0.4  )


map1