Required packages

Install the following packages using the CRAN

library(geofacet)
library(ggplot2)
library(sf)
library(ggplot2)
library(tidyr)
library(dplyr)
library(ggthemes)
library(ggrepel)
library(stringr)
library(rio)
library(tidycensus)
library(readxl)
library(tidyverse)
library(tmap)
library(tidycensus)
library(tigris)
library(maps)
library(leaflet)
library(biscale)
library(cowplot)
library(ggridges)
library(ggpubr)
library(ggspatial)

Data Download

Both the FEMA National Risk Index https://hazards.fema.gov/nri/ ranks and The CDC Social Vulnerability Index (SVI) is freely available from these websites. All tract level FEMA data can be downloaded from here https://www.atsdr.cdc.gov/placeandhealth/svi/index.html https://hazards.fema.gov/nri/Content/StaticDocuments/DataDownload//NRI_Shapefile_CensusTracts/NRI_Shapefile_CensusTracts.zip and the CDC data https://www.atsdr.cdc.gov/placeandhealth/svi/data_documentation_download.html

### Once you have downloaded the data, simply alter the directory to location of where the data is stored on your computer
risk <- st_read("NRI_Shapefile_CensusTracts.shp")
## Reading layer `NRI_Shapefile_CensusTracts' from data source `NRI_Shapefile_CensusTracts.shp' using driver `ESRI Shapefile'
## Simple feature collection with 72739 features and 368 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -19942590 ymin: 2144435 xmax: 20012850 ymax: 11537130
## projected CRS:  WGS 84 / Pseudo-Mercator
socialvuln <- st_read("SVI2018_US_tract.shp")
## Reading layer `SVI2018_US_tract' from data source `SVI2018_US_tract.shp' using driver `ESRI Shapefile'
## Simple feature collection with 72837 features and 126 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -179.1489 ymin: 18.91036 xmax: 179.7785 ymax: 71.36516
## geographic CRS: NAD83
### We also need a few data layers for making figure 2. These are not mandatory for the map, but
### can be downloaded easily e.g., https://www.glc.org/greatlakesgis
### We need shapefiles for the lakes 
gl <- st_transform(st_read("GreatLakes.shp"),4326)
## Reading layer `GreatLakes' from data source `GreatLakes.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 11 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -92.19876 ymin: 41.37988 xmax: -70.4547 ymax: 49.01076
## geographic CRS: WGS 84
### We also need the boundaries of the States 
gl_s<- st_transform(st_read("GL_states.shp"),4326)
## Reading layer `GL_states' from data source `GL_states.shp' using driver `ESRI Shapefile'
## Simple feature collection with 8 features and 17 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -10824620 ymin: 4434968 xmax: -7990234 ymax: 6340334
## projected CRS:  Popular Visualisation CRS / Mercator

Combine datasets

We need to combine the datasets for our analysis. We drop the geometry of the data for faster processing

### We combine social environmental risk data
contRisk <- full_join(st_drop_geometry(risk), socialvuln, by = c("TRACTFIPS" = "FIPS"))
###We also create a list of the abbreviations of CONUS States for later filtering 
conus_states <- c("AL", "AZ", "AR", "CA", "CO", "CT", "DE", "FL", "GA", "ID", "IL", "IN", "IA", "KS", "KY", "LA", "ME", "MD", "MA", "MI", "MN", "MS", "MO", "MT", "NE", "NV","NH", "NJ", "NM", "NY", "NC", "ND", "OH", "OK", "OR", "PA", "RI", "SC", "SD", "TN", "TX", "UT", "VT", "VA", "WA", "WV", "WI", "WY")
###We also create a list of GLISA States
my_states <- c("IL", "IN","NY","MI", "MN", "PA", "OH", "WI")

Data processing

Before visualizing the data, we needed to correct several data anomalies, and convert to measures for easy comparison and visualization.

### RISK_SCORE is the combined score calculated for the National Risk Index for Natural Hazards (FEMA) 
### No data is represented as -9999.0 in this dataset. We convert these to NAs 
contRisk$RISK_SCORE[contRisk$RISK_SCORE == -9999.000] <- NA
### The variable that assesses "Total Social vulnerability" is RPL_THEMES in the CDC/ATSDR dataset
### We similarly give values lower than O an NA
contRisk$RPL_THEMES[contRisk$RPL_THEMES < 0.000] <- NA
##here we create a percentage measure for the overall tract summary ranking variable:RPL_THEMES 
contRisk$RPL_THEMES <-contRisk$RPL_THEMES * 100
### Finally, we create a variable indicating whether the State is in GLISA area or not "Other States" ,"GLISA"
contRisk <- contRisk %>%
mutate( GLISAcat =  ifelse( STATEFIPS != "26" & STATEFIPS != "17" & STATEFIPS != "36" & STATEFIPS != "55" & STATEFIPS != "39" & STATEFIPS != "42" & STATEFIPS != "18","Other States" ,"GLISA" ))

Boxplot (Figure 1)

For the study, we were interested in visualizing the risk scores for GLISA states compared to all other States. We created box plots of this risk for each census tract for all states for the FEMA data, and the SVI (CDC). The code below specifies these variables, and makes the data comparable across the states. We also indicate GLISA States in Orange and give the mean line across the State in blue

### National Risk Index for Natural Hazards (FEMA)
Fema <- contRisk %>%
  filter(STATEABBRV %in% conus_states) %>%
  ggplot(aes(x=reorder(STATE, RISK_SCORE, na.rm = TRUE), y = RISK_SCORE, fill = GLISAcat)) +
  geom_boxplot() +
  #facet_wrap(~GLISAcat,ncol = 2) +
  ylab("FEMA National Risk") +
  scale_fill_manual(values=c( "#E69F00", "#999999"))+
  theme_classic()+
  theme( axis.text.x=element_blank(),
        legend.title = element_blank(),axis.title.x=element_blank(),
        legend.position = c(0.1, 0.90)) +
  geom_hline(yintercept = mean(contRisk$RISK_SCORE, na.rm = T), color="#56B4E9", size = 0.6)

### CDC/ATSDR dataset
Social <- contRisk %>%
  filter(STATEABBRV %in% conus_states) %>%
  ggplot(aes(x=reorder(STATE, RISK_SCORE, na.rm = TRUE), y = RPL_THEMES, fill = GLISAcat)) +
  geom_boxplot() +
  xlab("Census tract per CONUS States") + ylab("CDC/ATSDR SVI") +
  scale_fill_manual(values=c( "#E69F00", "#999999"))+
  theme_classic()+
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 6),
        legend.position="none") +
  geom_hline(yintercept = mean(contRisk$RPL_THEMES, na.rm = T), color="#56B4E9", size = 0.6)

### This function combine the two box-plots
ggarrange(Fema, Social,
          #labels = c("A", "B"),
          nrow = 2)

Bivariate map (Figure 2)

Our goal in the paper was also map these scores to understand clusters and diversity of environmental stressors and social vulnerabilities across the States. We use a bivariate map approach, which simplifies these scores in to high, medium and low classes, and combines these.

### First we need to create 3 classes for each of the scores 
### We used a "quantile" approaches for defining the classes (high, medium and low classes)
data <- contRisk %>%
  filter(STATEABBRV %in% my_states) %>%
  ##here we combine the classes e.g. High-high for the FEMA and CDC scores 
  bi_class(x = RISK_SCORE, y = RPL_THEMES, style = "quantile", dim = 3) %>%
  st_sf() %>%
  st_transform(4326)

### This creates the bivariate legend 
legend <- bi_legend(pal = "DkCyan",
                    dim = 3,
                    xlab = "FEMA Risk ",
                    ylab = "CDC SVI ",
                    size = 10)

### Here we use ggplot to create the map. This incorperates different layers
### and map elements
map <- data %>%
  st_sf() %>%
  ggplot() +
  geom_sf(mapping = aes(fill = bi_class), color = "white", size = NA, show.legend = FALSE) +
  bi_scale_fill(pal = "DkCyan", dim = 3) +
  geom_sf(data = gl_s, fill  = NA) +
  geom_sf(data = gl, fill = "aliceblue") +
  annotation_scale(location = "bl", width_hint = 0.5) +
  annotation_north_arrow(location = "bl", which_north = "true", pad_x = unit(0.75, "in"), pad_y = unit(0.5, "in"), style = north_arrow_fancy_orienteering) +
  coord_sf() +
  xlab("Longitude") +
  ylab("Latitude") +
  ggtitle("Map of Social and Natural Risk: Great Lakes Region") +
  theme(panel.grid.major = element_line(color = gray(0.5), linetype = "dashed", size = 0.5), panel.background = element_rect(fill = "#fafafa"))

###finally, we combine the map with the bivariate legend
finalPlot <- ggdraw() +
  draw_plot(map, 0, 0, 1, 1) +
  draw_plot(legend, 0.62, 0.7, 0.2, 0.2)
## Scale on map varies by more than 10%, scale bar may be inaccurate
finalPlot