This R Markdown file details the coding documentation for the ECO Action Mapping Project conducted by Nick An from October, 2022 - March, 2023. More details on the actual data itself can be found in a separate data documentation file. The four variables mapped from this dataset include heat and flood factors (from Risk Factor) and low income and ozone percentiles (from EJ index).

Loading packages that are required to create the maps.

pacman::p_load(tidyverse,         # general data wrangling
               tidycensus,        # importing Census attribute data into R
               sf,                # Spatial data classes
               tmap,              # Mapping/cartography
               leaflet,           # creating maps over open street maps 
               viridis,           # color maps 
               leaflet.extras     #extra leaflet feature to add search bar for map
)

Reading the csv files for this project

#EJ Screen (from EJ Index), Flooding data (from Risk Factor), and Heat data (from Risk Factor)

EJSCREEN <- read.csv ("EJSCREEN.csv")
flood <- read.csv ("flood.csv")
heat <- read.csv ('heat.csv')

Acquiring the geography data for Georgia

#Georgia's tract and block data were obtained from the American Community Survey (ACS). 
#This was done so that data from EJScreen and Risk Factor could be combined into a dataset with spatial info.
#This ACS dataset is a combination of tract and block levels, so any additional data using this code and dataset must have this combo. 
#The GEOID variable was recoded as numeric so it would match the other datasets when joining. 

georgia <- get_acs(
  geography = "tract",
  variables = "B19013_001",
  year = 2020,
  state = "GA",
  geometry = TRUE)
## Getting data from the 2016-2020 5-year ACS
## Downloading feature geometry from the Census website.  To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.
georgia2 <- georgia %>% mutate (GEOID = as.numeric(GEOID))

Coding flooding and heat variables

#The "fips" variable was renamed as GEOID for both flood and heat variables so that it matches the GEOID name for the final dataset when merging. 
#For both heat and flood variables, "major_heat" and "major_flood" variables were combined by combining the       "count_floodfactor" and "count_heatfactor" variables that were 5 or higher. 
#This was done by using information from the Risk Factor website where the authors assigned properties that had a score of "5" or higher to be at "major" risk or higher (severe and extreme). 
#Only flood and heat factors that had scores of 5 or higher were considered to only focus on properties that were at a greater degree of risk. 

flood2 <- rename (flood, "GEOID" = "fips")
flood3 <- mutate (flood2, major_flood = count_floodfactor5 +
                    count_floodfactor6  + 
                    count_floodfactor7  + 
                    count_floodfactor8  + 
                    count_floodfactor9 + 
                    count_floodfactor10)
flood4 <- mutate (flood3, p_major = ntile(major_flood, 5))

heat2 <- rename (heat, "GEOID" = "fips")
heat3 <- mutate (heat2, major_heat = count_heatfactor5 +
                    count_heatfactor6  + 
                    count_heatfactor7  + 
                    count_heatfactor8  + 
                    count_heatfactor9 + 
                    count_heatfactor10)
heat4 <- mutate (heat3, p_majorheat = ntile(major_heat, 5))

Merging all datasets into one datasets

#EJScreen was recoded so that the "ID" variable matches the "GEOID" variable of all the other datasets.
#All the datasets were joined using the left_join function by the "GEOID" variable. 
#The final combined dataset was projected to 4326 (WGS84) so that it matches the coordinate system used by the             leaflet package 
#The final dataset (totalmap.aea2) includes only metro atlanta counties.

EJSCREEN2 <- rename (EJSCREEN, "GEOID" = "ID")

EJmap <- georgia2 %>% left_join (EJSCREEN2, by = 'GEOID')
EJfloodmap <- EJmap %>% left_join (flood4, by = "GEOID")
totalmap <- EJfloodmap %>% left_join (heat4, by = "GEOID")

totalmap.aea <- st_transform(totalmap, 4326)
totalmap.aea2 <- filter (totalmap.aea, CNTY_NAME == "DeKalb" | CNTY_NAME == "Fulton" | CNTY_NAME == "Cherokee" | CNTY_NAME == "Clayton" | CNTY_NAME =="Cobb" |CNTY_NAME == "Douglas" | CNTY_NAME =="Fayette" |CNTY_NAME == "Forsyth" | CNTY_NAME =="Gwinnett" | CNTY_NAME =="Rockdale" |CNTY_NAME == "Henry")

Proctor Creek and Intrenchment Creek Watershed Boundaries

#Proctor Creek Watershed Boundaries. 
#The boundary was created by looking at the watershed boundary image from EPA's waterway website. 
#Points were mapped on google maps to create a rough sketch of the same boundary, and latitude and longitude            information was obtained from these points and inputted into R. 
#These points were stored into the "s1" object. 

 s1 <- "[(-84.39024, 33.75442),
(-84.39363, 33.75842),
(-84.39751, 33.76374),
(-84.39963, 33.77027),
(-84.40604, 33.77417),
(-84.40886, 33.77565),
(-84.41195, 33.78103),
(-84.42538, 33.78688),
(-84.42999, 33.78513),
(-84.43351, 33.78624),
(-84.43514, 33.78806),
(-84.44116, 33.79804),
(-84.44678, 33.80588),
(-84.44524, 33.81305),
(-84.44704, 33.82093),
(-84.45496, 33.82642),
(-84.45976, 33.82856),
(-84.46311, 33.83369),
(-84.47247, 33.83996),
(-84.48663, 33.84346),
(-84.49281, 33.84866),
(-84.49641, 33.84987),
(-84.49839, 33.84831),
(-84.49942, 33.84374),
(-84.49744, 33.84096),
(-84.49135, 33.8339),
(-84.49332, 33.8302),
(-84.49873, 33.81465),
(-84.51587, 33.80358),
(-84.53389, 33.78503),
(-84.52393, 33.78218),
(-84.50145, 33.78302),
(-84.48222, 33.78289),
(-84.47158, 33.77376),
(-84.46059, 33.75207),
(-84.43673, 33.74824),
(-84.42471, 33.73882),
(-84.408, 33.73382),
(-84.40641, 33.74248)]"

#unwanted characters stripped out 

c1 <- chartr('()[]','    ', s1)

#cleaned string converted into two column matrix so the points can easily be mapped and connected with Leaflet 

m1 <- matrix(as.numeric(strsplit(c1,",")[[1]]),ncol=2,byrow=TRUE)


#Intrenchment Creek Watershed Boundaries. Same process as Proctor Creek Watershed Boundaries. 

s2 <- "[(-84.38986, 33.75492),
(-84.40641, 33.74248),
(-84.408, 33.73382),
(-84.39073, 33.72085),
(-84.37322, 33.71384),
(-84.35949, 33.7083),
(-84.35425, 33.70394),
(-84.34988, 33.69637),
(-84.35267, 33.68047),
(-84.35757, 33.67125),
(-84.35817, 33.66904),
(-84.35508, 33.66447),
(-84.33714, 33.65782),
(-84.31833, 33.65546),
(-84.31618, 33.6536),
(-84.31335, 33.65403),
(-84.30245, 33.66017),
(-84.29009, 33.66623),
(-84.28073, 33.66638),
(-84.27696, 33.67323),
(-84.25858, 33.67466),
(-84.25858, 33.67845),
(-84.26501, 33.68131),
(-84.27188, 33.69123),
(-84.26862, 33.70116),
(-84.28019, 33.72837),
(-84.27882, 33.73358),
(-84.28723, 33.73901),
(-84.29375, 33.74707),
(-84.29298, 33.75178),
(-84.29813, 33.75449),
(-84.30173, 33.7629),
(-84.3134, 33.76561),
(-84.32199, 33.76053),
(-84.33632, 33.76276),
(-84.3498, 33.75946),
(-84.35997, 33.75439),
(-84.38186, 33.74983),
(-84.38959, 33.74875)]"

#unwanted characters stripped out 

c2 <- chartr('()[]','    ',s2)

#cleaned string converted into two column matrix so the points can easily be mapped and connected with Leaflet 

m2 <- matrix(as.numeric(strsplit(c2,",")[[1]]),ncol=2,byrow=TRUE)

Finding quantiles for heat and flooding

#Quantile function was done to identify the quantiles for floods and heat that were "major." 
#For flooding, the quantiles were 0, 46, 116, 258, and 18,649. 
#For heat, the quantiles were 0, 0, 668, 1571, and 62775. This was done using the national datasets for flood and heat risk, so that the same cutoffs for data at the national level can be applied to Georgia. 
#However, since there are some outliers using these numbers, the maximum cutoff was found using the max function for the Georgia datasets. 

quantile(flood3$major_flood)
##    0%   25%   50%   75%  100% 
##     0    46   116   258 18649
quantile(heat3$major_heat)
##    0%   25%   50%   75%  100% 
##     0     0   668  1571 62775
max(totalmap.aea$major_flood, na.rm = TRUE)
## [1] 3365
max(totalmap.aea$major_heat, na.rm = TRUE)
## [1] 6414

Assigning colors for variables

#The "colorBin" and "colorQuantile" functions were used to assign chloropeth colors for each of the 4 variables               of this map: flood factor, heat factor, low income percentiles, and ozone percentiles. 
#"colorBin" was used for flood and heat factors using the quantiles assigned from the previous code. 
#Major floods were designated with a blue scale and major heat was designated with a red scale. 
#colorQuantile" function was used for the already existing quantiles from the EJ index for low income percentiles and ozone percentiles.
#Low income percentiles was categorized using a gray scale for 7 categories and ozone percentiles was categorized using a purple scale for 5 categories.

  pal <- colorBin ("Blues",
                     domain = totalmap.aea2$major_flood, bin = c(0, 46, 116, 258, 3365))
  pal2 <- colorQuantile ("Greys",
                      domain = totalmap.aea2$P_LWINCPCT,
                      n=7)
  pal3 <- colorBin ("Reds",
                        domain = totalmap.aea2$major_heat, bin= c(0, 668, 1571, 6414))
  pal4 <- colorQuantile ("Purples",
                         domain = totalmap.aea2$P_OZONE_D2, n=5)

Adding popup features for census tracts

#These popup features were added using the paste0() function so that users could click on each tract and block groups to find the GEOID and values for each variable. 
  mypopup <- paste0("GEOID: ", totalmap.aea2$GEOID, "<br>", "Flood Factor: ", totalmap.aea2$major_flood)
  mypopup3 <- paste0("GEOID: ", totalmap.aea2$GEOID, "<br>", "Heat Factor: ", totalmap.aea2$major_heat)
  mypopup2 <- paste0("GEOID: ", totalmap.aea2$GEOID, "<br>", "Low Income Percentile: ", totalmap.aea2$P_LWINCPCT)
  mypopup4 <- paste0("GEOID: ", totalmap.aea2$GEOID, "<br>", "Ozone Percentile: ", totalmap.aea2$P_OZONE_D2)

Mapping the dataset

#The leaflet function was used to map the combined dataset for the 4 climate change variables. 
#The dataset was inputed into the leaflet () function. 
#Open street Maps ("OSM") was added using the addTiles function so that the maps would be overlayed with an easy to navigate street map.
#Polygons were added for each of the 4 variables, with the fillColor function being used to create the chloropeth colors on the map based on the previous code where I assigned the colors. 
#In addition, the watershed boundaries that I created from the earlier code was inputted, also using the                 addPolygons function.
#Legends were added for each of the variables and formatted in a way (ex: bottom right, bottom left) so that they would not overlap with each other using the addLegend function.
#The function addLayersControl was used so that the users could select if they want each of the 4 variables to                appear over the open street map or not. 
#hideGroup() function was done for each of the 4 variables so that when a user first opens the map, all                       the variables are hidden so that they are not overwhelmed by seeing all of them. 
#Finally, an addSearchOSM() function was added so that users could type in specific addresses that the map can                take the user to. 
  leaflet(totalmap.aea2) %>%
    addTiles(group = "OSM") %>%
    addPolygons (group = "flood",
                 stroke = FALSE, 
                 fillColor = ~pal(major_flood),
                 fillOpacity = 0.4,
                 popup = mypopup)   %>%
    addPolygons (group = "low income",
                 stroke = FALSE, 
                 fillColor = ~pal2(P_LWINCPCT),
                 fillOpacity = 0.4,
                 popup = mypopup2) %>%
    addPolygons (group = "heat",
                 stroke = FALSE, 
                 fillColor = ~pal3(major_heat),
                 fillOpacity = 0.4,
                 popup = mypopup3) %>%
    addPolygons (group = "ozone",
                 stroke = FALSE, 
                 fillColor = ~pal4(P_OZONE_D2),
                 fillOpacity = 0.4,
                 popup = mypopup4) %>%
    addPolygons(data=m1, weight=5,
                fillColor = "Transparent",
                label = "Proctor Creek Watershed") %>%
    addPolygons(data=m2, weight=5,
                fillColor = "Transparent",
                label = "Intrenchment Creek Watershed") %>%
    addLegend(pal = pal, 
              values = ~major_flood, 
              opacity = 0.4, 
              title = "Flood Factor",
              position = "bottomright",
              group = "flood") %>%
    addLegend(pal = pal2, 
              values = ~P_LWINCPCT, 
              opacity = 0.4, 
              title = "Low Income Percentiles",
              position = "bottomright",
              group = "low income") %>%
    addLegend(pal = pal3, 
              values = ~major_heat, 
              opacity = 0.4, 
              title = "Heat Factor",
              position = "bottomleft",
              group = "heat") %>%
    addLegend(pal = pal4, 
              values = ~P_OZONE_D2, 
              opacity = 0.4, 
              title = "Ozone Percentiles",
              position = "bottomleft",
              group = "ozone") %>%
    addLayersControl(baseGroups = c("OSM"), 
                     overlayGroups = c("flood", "low income", "heat", "ozone")) %>%
    hideGroup ("flood") %>%
    hideGroup ("low income") %>%
    hideGroup ("heat") %>%
    hideGroup ("ozone") %>%
    setView(lng = -84.408868, lat = 33.760163,zoom = 11.5) %>%
    addSearchOSM() 
Leaflet | © OpenStreetMap contributors, CC-BY-SA

Other notes: The final leaflet map allows the user to zoom in and zoom out, move around the map, search for locations of interest, and select for which of the 4 variables they are interested in, making the map user friendly. The entire EJ index dataset is included in my combined dataset, so other variables can be examined to be incorporated into the map. In addition, other datasets can also be added into the existing dataset (as long as they have the same distribution of tract and block groups) to add more variables into the map (ex: PM2.5).