Evacuation_NewHampshirePDF.Rmd

---
title: "Evacuation Risks in New Hampshire"
author: "Marcos Luna and Neenah Estrella-Luna"
date: "`r Sys.Date()`"
output: 
  bookdown::pdf_document2:
    toc: true
    toc_depth: 3
    number_sections: true
    fig_caption: yes        
    includes:  
      in_header: my_header.tex
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo=FALSE, message=FALSE, warning=FALSE)

library(tidyverse)
library(sf)
library(tmap)
library(tmaptools)
library(lwgeom)
library(tigris)
options(tigris_use_cache = TRUE, tigris_class = "sf")
```
\pagebreak
# Analysis of evacuation-related risks in New Hampshire

This is an analysis of flood or storm surge-related evacuation risks for priority populations in New Hampshire In the event of significant inland or coastal flooding, evacuation may be required. For individuals and households with limited mobility, either due to inadequate access to transportation options or because of physical limitations, evacuation presents heightened risk. Evacuation may also prove especially difficult for individuals and households due to limited economic resources, difficulty understanding or accessing information, or low trust in official sources of information. 

## Analysis of Flood Hazard Exposure in New Hampshire

```{r NFHZAdata, include=FALSE, cache=TRUE}
# Analysis of flooding and hurricane evacuation risks for priority populations in New Hampshire
load("DATA/ne_layers.rds")

# Extract state census units and convert to projected local CRS EPSG:2823: NAD83(HARN) / New Hampshire See https://spatialreference.org/ref/epsg/2823/
nh_blkgrp_2823 <- ne_blkgrp_sf %>% 
  filter(STATE == "New Hampshire") %>% 
  st_transform(., crs = 2823)

nh_tracts_2823 <- ne_tracts_sf %>% 
  filter(STATE == "New Hampshire") %>% 
  st_transform(., crs = 2823)

# Get rid of empty geometries
empty_geo <- st_is_empty(nh_blkgrp_2823)
nh_blkgrp_2823 <- nh_blkgrp_2823[!empty_geo,]
empty_geo <- st_is_empty(nh_tracts_2823)
nh_tracts_2823 <- nh_tracts_2823[!empty_geo,]

# Extract state boundary
nh_state_sf <- ne_states_sf_cb %>% 
  filter(NAME == "New Hampshire") %>% 
  st_transform(., crs = 2823)

# Extract state countiies
nh_counties_sf <- counties(state = "NH", cb = TRUE) %>% 
  st_transform(., crs = 2823)

# #### Processing in ArcGIS #####
# # Write out block groups for processing in ArcGIS
# nh_blkgrp_2823 %>%
#   dplyr::select(GEOID) %>%
#   st_write(., "DATA/FEMA/NH/nh_blkgrp_2823.shp", delete_layer = TRUE)
# #
# # repeat for tracts
# nh_tracts_2823 %>%
#   dplyr::select(GEOID) %>%
#   st_write(., "DATA/FEMA/NH/nh_tracts_2823.shp", delete_layer = TRUE)
# 
# # repeat for state boundary
# nh_state_sf %>%
#   st_write(., "DATA/FEMA/NH/nh_state_2823.shp", delete_layer = TRUE)
# 
# # repeat for state counties
# nh_counties_sf %>% 
#   st_write(., "DATA/FEMA/NH/nh_counties_2823.shp", delete_layer = TRUE)

# Use dasymetric mapping to calculate populations within flood zones. Approach follows method used by Qiang (2019) to eliminate unpopulated areas of census polygons and then reallocate populations to developed areas as identified in National Land Cover Dataset (NLCD). 
# Perform NLCD raster-to-vector conversion, vector erase/difference, and vector intersections in ArcMap because it takes too long in R. 
# In ArcMap:
# Convert NLCD raster to shapefile. Clip to state.
# Isolate undeveloped areas (gridcode NOT 22 - 24). 
# Erase areas of nh_blkgrp_2823 and nh_tracts_2823 that overlap with undeveloped areas in NLCD shapefiles. Compute OldArea of erased polygons in sqm to identify area of developed polygons remaining. 
# Intersect erased nh_blkgrps and erased nh_tracts with NFHZA and Hurricane evacuation zones. Read back into R.

##### Return to working in R ######
# # read in processed nh_blkgrps and nh_tracts
# st_layers(dsn = "DATA/FEMA/NH")
nh_blkgrps_developed <- st_read(dsn = "DATA/FEMA/NH",
                                layer = "nh_blkgrps_developed") %>% 
  left_join(., as.data.frame(nh_blkgrp_2823), by = "GEOID") %>% 
  st_transform(., crs = 2823) %>% 
  mutate(NewArea = st_area(.)) %>%
  st_make_valid() %>% 
  mutate(Proportion = as.numeric(NewArea/OldArea),
         NewPop = totalpopE*Proportion,
         NewMinority = minorityE*Proportion,
         NewUnder5 = under5E*Proportion,
         NewOver64 = over64E*Proportion,
         NewUnder18 = under18E*Proportion,
         NewEng_limit = eng_limitE*Proportion,
         NewPov = num2povE*Proportion,
         NewLths = lthsE*Proportion)

nh_blkgrps_nfhza <- st_read(dsn = "DATA/FEMA/NH", 
                            layer = "nh_blkgrps_nfhza") %>% 
  left_join(., as.data.frame(nh_blkgrp_2823), by = "GEOID") %>% 
  st_transform(., crs = 2823) %>% 
  mutate(NewArea = st_area(.)) %>%
  st_make_valid()

nh_tracts_nfhza <- st_read(dsn = "DATA/FEMA/NH", 
                           layer = "nh_tracts_nfhza") %>% 
  left_join(., as.data.frame(nh_tracts_2823), by = "GEOID") %>% 
  st_transform(., crs = 2823) %>% 
  mutate(NewArea = st_area(.)) %>%
  st_make_valid()

# Apportion populations based on geographic proportion of intersect
nh_blkgrps_nfhza <- nh_blkgrps_nfhza %>%
  mutate(Proportion = as.numeric(NewArea/OldArea),
         NewPop = totalpopE*Proportion,
         NewMinority = minorityE*Proportion,
         NewUnder5 = under5E*Proportion,
         NewOver64 = over64E*Proportion,
         NewUnder18 = under18E*Proportion,
         NewEng_limit = eng_limitE*Proportion,
         NewPov = num2povE*Proportion,
         NewLths = lthsE*Proportion)

nh_tracts_nfhza <- nh_tracts_nfhza %>%
  mutate(Proportion = as.numeric(NewArea/OldArea),
         NewDisabled = disabledOver18E*Proportion,
         NewNoCar = HHnoCarE*Proportion)

# Compute total block group populations within flood zones
nh_flood_blkgrp_df <- nh_blkgrps_nfhza %>%
  as.data.frame() %>%
  summarize(`Total Pop` = as.integer(sum(NewPop)),
            Minority = as.integer(sum(NewMinority)),
            `Under 5` = as.integer(sum(NewUnder5)),
            `Over 64` = as.integer(sum(NewOver64)),
            `Under 18` = as.integer(sum(NewUnder18)),
            `Limited English HH` = as.integer(sum(NewEng_limit)),
            `Low Income` = as.integer(sum(NewPov)),
            `No HS Dip` = as.integer(sum(NewLths))) %>%
  gather(key = Group, value = FloodPop)

nh_flood_tracts_df <- nh_tracts_nfhza %>%
  as.data.frame() %>%
  summarize(`Disabled` = as.integer(sum(NewDisabled)),
            `No Car HH` = as.integer(sum(NewNoCar))) %>%
  gather(key = Group, value = FloodPop)

# Compute total tract populations within the state for same groups
nh_tract_flood_pops_df <- nh_tracts_2823 %>%
  as.data.frame() %>%
  summarize(`Disabled` = sum(disabledOver18E),
            `No Car HH` = sum(HHnoCarE)) %>%
  gather(key = Group, value = NHPop) %>%
  left_join(.,nh_flood_tracts_df, by = "Group")

# Compute populations for state,and join with flood pops
nh_FloodPops_df <- nh_blkgrp_2823 %>%
  as.data.frame() %>%
  summarize(`Total Pop` = sum(totalpopE),
            Minority = sum(minorityE),
            `Under 5` = sum(under5E),
            `Over 64` = sum(over64E),
            `Under 18` = sum(under18E),
            `Limited English HH` = sum(eng_limitE),
            `Low Income` = sum(num2povE),
            `No HS Dip` = sum(lthsE)) %>%
  gather(key = Group, value = NHPop) %>%
  left_join(., nh_flood_blkgrp_df, by = "Group") %>%
  rbind(.,nh_tract_flood_pops_df) %>%
  mutate(PctFlood = FloodPop/NHPop*100)

nh_towns_sf <- county_subdivisions(state = "NH", cb = TRUE) %>% 
  st_transform(., crs = 2823)

town_tracts_df <- nh_tracts_2823 %>% 
  dplyr::select(-GEOID) %>% 
  st_centroid(.) %>% 
  st_intersection(nh_towns_sf,.) %>% 
  as.data.frame() %>% 
  group_by(NAME) %>%
  summarize(TotalDisabled = sum(disabledOver18E, na.rm = TRUE),
            TotalOver18 = sum(Over18E, na.rm = TRUE),
            TotalNoCar = sum(HHnoCarE, na.rm = TRUE))

townpopsdf <- nh_blkgrp_2823 %>% 
  dplyr::select(-GEOID) %>% 
  st_centroid(.) %>% 
  st_intersection(nh_towns_sf,.) %>% 
  as.data.frame() %>% 
  group_by(NAME) %>% 
  summarize(GEOID = unique(GEOID),
            TotalPop = sum(totalpopE, na.rm = TRUE),
            TotalHH = sum(householdsE, na.rm = TRUE),
            TotalLEH = sum(eng_limitE, na.rm = TRUE),
            TotalMin = sum(minorityE, na.rm = TRUE),
            TotalLowInc = sum(num2povE, na.rm = TRUE),
            TotalNoHS = sum(lthsE, na.rm = TRUE),
            TotalOver64 = sum(over64E, na.rm = TRUE),
            TotalUnder5 = sum(under5E, na.rm = TRUE),
            TotalUnder18 = sum(under18E, na.rm = TRUE)) %>% 
  left_join(., town_tracts_df, by = "NAME")



nh_state_sf_cb <- ne_states_sf_cb %>% 
  filter(NAME == "New Hampshire") %>% 
  st_transform(., crs = st_crs(nh_blkgrp_2823))

# save data for later analysis
save(nh_blkgrps_nfhza,nh_tracts_nfhza,
     nh_FloodPops_df,
     file = "DATA/FEMA/NH/nfhza_census.Rds")
```

As a humid, coastal state, a significant portion of New Hampshire's land area and population are exposed to the risk of flooding from overbanking of inland water bodies (e.g., ponds and rivers) or from coastal storm surge and sea level rise. 

The analysis of flood exposure presented here is based on the Federal Emergency Management Agency's (FEMA) National Flood Hazard Layer (NFHL), a digital version of FEMA's most recent flood maps.[^NFHL]  FEMA identifies areas subject to varying levels of flood risk through its Flood Insurance Rate Maps (FIRMs), which are produced for most parts of the country to identify areas subject to flooding. FEMA's FIRMs are a national standard used by all federal agencies for the purposes of requiring and rating the purchase of flood insurance and regulating new development.  Population data comes from the American Community Survey (ACS) 2018 5-year estimates at the census tract and block group levels. 

A significant proportion of New Hampshire's population lives near the coast or near inland water bodies. As a result, a significant proportion of the population lives near or within known flood zones. Figure \@ref(fig:mapNFHZAPop) below shows the population distribution and areas subject to floods with an Annual Exceedance Probability (AEP) of 1% (also known as a '100-year' flood) and areas subject to 0.2% AEP (also known as a '500-year' flood). Areas within the 1% AEP flood zone are designated by FEMA as Special Flood Hazard Areas, and development within those zones must be covered by flood insurance. Areas within the 0.2% AEP are not currently regulated, but these areas are nevertheless subject to flood risk under more extreme, albeit less frequent, flooding circumstances. Coastal areas within the 0.2% AEP zone are a reasonable proxy measure for the risks of sea level rise. Note that portions of northwest and northeast New Hampshire are not currently mapped for flooding risk by FEMA. 
```{r mapNFHZAdata, include=FALSE, cache=TRUE}
# Process data for mapping 
# Create point layer of major cities for context
# Note cb=FALSE is necessary for extracting centroids from town polygons. Otherwise, if cb=TRUE, cannot extract centroids from multipolygon features.
# Read in NFHL for NH. Data comes from FEMA. 
# List available layers in geodatabase
# st_layers("DATA/FEMA/NH/NFHL_33_20190723.gdb")
# # Read in flood hazard areas
# nh_nfhza_2823 <- st_read(dsn = "DATA/FEMA/NH/NFHL_33_20190723.gdb", layer = "S_Fld_Haz_Ar", quiet = TRUE) %>%
#   filter(!FLD_ZONE %in% c("OPEN WATER","AREA NOT INCLUDED","D") &
#            !ZONE_SUBTY %in% c("AREA OF MINIMAL FLOOD HAZARD",
#                               "AREA WITH REDUCED FLOOD RISK DUE TO LEVEE")) %>%
#   mutate(FLD_ZONE = as.character(FLD_ZONE)) %>% # omit unused factor levels
#   st_transform(., crs = 2823) %>%
#   st_make_valid() %>%
#   group_by(FLD_ZONE) %>% # aggregate flood zone polygons
#   summarize(count = n()) %>%
#   mutate(Area = st_area(.),
#          Interval = case_when(
#            FLD_ZONE  == "A" ~ "100-year",
#            FLD_ZONE  == "AE" ~ "100-year",
#            FLD_ZONE  == "AH" ~ "100-year",
#            FLD_ZONE  == "AO" ~ "100-year",
#            FLD_ZONE  == "VE" ~ "100-year",
#            FLD_ZONE == "X" ~ "500-year"))
# 
# # get rid of empty geometries
# empty_geo <- st_is_empty(nh_nfhza_2823)
# nh_nfhza_2823 <- nh_nfhza_2823[!empty_geo,]
# 
# # save this layer since it takes more than 30 mins to load
# saveRDS(nh_nfhza_2823,"DATA/FEMA/NH/nh_nfhza_2823.Rds")
# load previously process flood layer
nh_nfhza_2823 <- readRDS("DATA/FEMA/NH/nh_nfhza_2823.Rds")


# download Census TIGERLine hydrography for NH
# # First, extract list of county names to use with tigris::water
# nh_counties <- counties("NH") %>%
#   pull(NAME)
# # Next, download water features for each county and rbind to one layer
# nh_awater_sf <- rbind_tigris(
#   lapply(
#     nh_counties, function(x) area_water(state = "NH", county = x)
#   )
# ) %>%
#   st_union() %>%
#   st_as_sf() %>%
#   st_transform(., crs = 2823)
# 
# # save time by loading file
# saveRDS(nh_awater_sf, "DATA/FEMA/NH/nh_awater_sf.Rds")
# load water features
nh_awater_sf <- readRDS("DATA/FEMA/NH/nh_awater_sf.Rds")

# # crop flood zones to land areas only
# start_time <- Sys.time()
# nh_nfhza_2823_land <- nh_nfhza_2823 %>%
#   crop_shape(., nh_state_sf, polygon = TRUE) %>%
#   st_difference(., nh_awater_sf) %>%
#   mutate(Area = st_area(.)) %>%
#   st_make_valid() # takes about 12 mins
# end_time <- Sys.time()
# end_time - start_time
# # save time by just loading file
# saveRDS(nh_nfhza_2823_land, file = "DATA/FEMA/NH/nh_nfhza_2823_land.Rds")
# load cropped flood zones
nh_nfhza_2823_land <- readRDS("DATA/FEMA/NH/nh_nfhza_2823_land.Rds")

# Create a dot density map of total populations and overlay on flood zones
# Create point layer of major cities for context
# Note cb=FALSE is necessary for extracting centroids from town polygons. Otherwise, if cb=TRUE, cannot extract centroids from multipolygon features.
nh_towns_sf_pts <- county_subdivisions(state = "NH", cb = TRUE) %>% 
  filter(NAME %in% c("Nashua",
                     "Portsmouth",
                     "Manchester",
                     "Concord",
                     "Laconia",
                     "Lebanon",
                     "Conway",
                     "Woodstock",
                     "Franconia",
                     "Berlin",
                     "Lancaster",
                     "Pittsburg",
                     "Keene")) %>% 
  st_transform(., crs = 2823) %>% 
  st_centroid(of_largest_polygon = TRUE)

# Create road layer for context
nh_highways <- primary_roads() %>% 
  filter(FULLNAME %in% c("I- 89","I- 91","I- 93","I- 95")) %>% 
  tmaptools::crop_shape(., ne_states_sf_cb) %>% 
  st_transform(., crs = 2823)

# nh_highways2nd <- primary_secondary_roads("NH") %>% 
#   filter(FULLNAME %in% c("US Hwy 202","Presidential Hwy","US Hwy 2")) %>% 
#   st_transform(., crs = 2823)

# Extract highway segments for labeling
I89roadSegment <- nh_highways %>% 
  filter(LINEARID == "1105281262324")

I91roadSegment <- nh_highways %>% 
  filter(LINEARID == "110373954766")

I95roadSegment <- nh_highways %>% 
  filter(LINEARID == "1105569136123")

I93roadSegment <- nh_highways %>% 
  filter(LINEARID == "1105598909781")

# Create custom icons of highway shields
I89 <- tmap_icons(file = "https://upload.wikimedia.org/wikipedia/commons/thumb/f/f4/I-89.svg/200px-I-89.svg.png")
I91 <- tmap_icons("https://upload.wikimedia.org/wikipedia/commons/thumb/9/90/I-91.svg/200px-I-91.svg.png")
I95 <- tmap_icons(file = "https://upload.wikimedia.org/wikipedia/commons/thumb/6/61/I-95.svg/200px-I-95.svg.png")
I93 <- tmap_icons("https://upload.wikimedia.org/wikipedia/commons/thumb/0/0d/I-93.svg/200px-I-93.svg.png")

# # Bring in outline of New York for western boundary
# ny_state_sf <- tigris::states(cb = TRUE) %>% 
#   filter(NAME == "New York")

# # Create polygons to show where there are no FEMA maps
# start_time <- Sys.time()
# nofemapolys <- st_read(dsn = "DATA/FEMA/NH/NFHL_33_20190723.gdb",
#                        layer = "S_Fld_Haz_Ar", quiet = TRUE) %>%
#   filter(FLD_ZONE != "AREA NOT INCLUDED") %>%
#   group_by() %>%
#   summarize() %>% # aggregate all flood zone polygons into one
#   st_transform(., crs = st_crs(nh_state_sf)) %>%
#   st_make_valid() %>%
#   st_difference(nh_state_sf,.) %>% # extract areas with no flood coverage
#   st_cast(., "POLYGON") %>% # disaggregate multipolygon feature
#   mutate(Area = st_area(.),
#          LABEL = "No FEMA\nFlood Data") %>% # calc area of each polygon
#   filter(as.numeric(Area)/10^6 > 1000) # remove polygons < 1000km2
# end_time <- Sys.time()
# end_time - start_time #takes about 2.4mins
# 
# # save time by loading previously processed file
# saveRDS(nofemapolys,"DATA/FEMA/NH/nofemapolys.Rds")
# load data
nofemapolys <- readRDS("DATA/FEMA/NH/nofemapolys.Rds")

# # Create random points, with 1 point for every 500 people
# nh_totalpop_pts <- nh_blkgrps_developed %>% 
#   select(NewPop) %>% 
#   filter(NewPop >= 500) %>% 
#   st_sample(., size = round(.$NewPop/500)) %>% # create 1 random point for every 500 people
#   st_sf(.) %>% 
#   group_by() %>% 
#   summarize() # reduce to one multipoint object

# Create a simplified version of NFHZA polygons to speed up mapping
# start_time <- Sys.time()
nh_nfhza_2823_land_simple <- nh_nfhza_2823_land %>% 
  select(Interval) %>% 
  st_cast(., "MULTIPOLYGON") %>% # homogenize type
  st_cast(., "POLYGON") %>% 
  st_simplify(., dTolerance = 100) %>% # reduce vertices
  st_make_valid() %>%
  group_by(Interval) %>% 
  summarize() %>% 
  mutate(Area = st_area(.))
# end_time <- Sys.time()
# end_time - start_time

# Create a simplified version of nh_awater_sf to speed up mapping
nh_awater_sf_simple <- nh_awater_sf %>% 
  st_cast(., "MULTIPOLYGON") %>% # homogenize type
  st_cast(., "POLYGON") %>% 
  st_simplify(., dTolerance = 100) %>% # reduce vertices
  st_make_valid() %>%
  group_by() %>% 
  summarize()

# read in Canadian province boundary files for background
canada <- st_read(dsn = "DATA/FEMA/NH", layer = "lpr_000b16a_e") %>% 
  filter(PRENAME == "Quebec") %>% 
  st_transform(., crs = 2823) %>% 
  st_simplify(., dTolerance = 100) %>%
  st_make_valid()
```

```{r mapNFHZAPop, echo=FALSE, message=FALSE, warning=FALSE, fig.align = "center", fig.cap="FEMA flood zones across New Hampshire."}
# Create a dot density map of total populations and overlay on flood zones
# Map totalpop and flood zones
m1 <- tm_layout(bg.color = "#e6f3f7") +
  tm_shape(nh_blkgrp_2823, unit = "mi") + tm_fill(col = "white") +
  tm_shape(ne_states_sf_cb) + tm_fill(col="white") +
  tm_shape(canada) + tm_fill(col = "white") +
  # tm_shape(ny_state_sf) + tm_fill(col="white") +
  tm_shape(nh_awater_sf_simple) + tm_fill(col = "#e6f3f7") +
  # tm_shape(ne_states_sf_cb) + tm_borders() +
  # tm_shape(ny_state_sf) + tm_borders() +
  # tm_shape(nh_totalpop_pts) + tm_dots(col = "darkgoldenrod3",
  #                                     labels = "1 dot = 500 people",
  #                                     alpha = 0.6) +
  tm_shape(nh_nfhza_2823_land_simple) +
  tm_fill(col = "Interval",
          palette = c("dodgerblue", "blueviolet"),
          labels = c("1% AEP (100-year)", "0.2% AEP (500-year)"),
          title = "FEMA Flood Zones",
          alpha = 0.6,
          border.alpha = 0) +
  tm_shape(ne_states_sf_cb) + tm_borders(lwd = 0.2, alpha = 0.8) +
  # tm_shape(ny_state_sf) + tm_borders(lwd = 0.8) +
  tm_shape(nh_highways) + tm_lines(col = "seashell4", lwd = 1.5) +
  # tm_shape(nh_highways2nd) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.15) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.15) +
  tm_shape(nh_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_shape(nofemapolys) + tm_fill(col = "gray", alpha = 0.4) +
  tm_text("LABEL", alpha = 0.5, size = .7) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5,
               position = c("left","TOP")) +
  # tm_add_legend(type = "fill", col = "darkgoldenrod3",
  #               alpha = 0.6,
  #               border.col = "white", border.alpha = 0,
  #               labels = "1 dot = 500 people",
  #               title = "Total Population") +
  tm_layout(title = "FEMA\nFlood Zones",
            frame = TRUE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m1, "DATA/FEMA/NH/nh_flood.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/FEMA/NH/nh_flood.png")
```

Approximately 5% of New Hampshire's land area falls within FEMA flood zones. The breakdown by type of flood zone is presented in Table \@ref(tab:tabNFHZAarea) below. 

```{r tabNFHZAarea, echo=FALSE, message=FALSE, warning=FALSE, fig.align = "center", fig.cap="New Hampshire Land Area within Flood Zones."}
# Total land in NH
nh_area <- as.numeric(nh_state_sf$ALAND)

# Total and percentage of land area of RI within flood zones
nh_nfhza_2823_land %>% 
  as.data.frame() %>% 
  group_by(Interval) %>% 
  summarize(SqKm = round(as.numeric(sum(Area)/10^6),1), 
            SqMi = round(as.numeric(SqKm/2.59),1),
            PctArea = paste0(as.character(round(as.numeric(sum(Area)/nh_area*100),1)),"%")) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','),
                    caption = "New Hampshire Land Area within Flood Zones",
                    digits = 0, align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header"))
```

These percentages approximate the proportion of the total population in the state that are living within these flood zones. However, exposure varies significantly by population subgroup. Table \@ref(tab:tabNFHZApop) below shows the total and percentages of the general population and various subgroups living within flood zones. Aside from the total population, the largest absolute numbers of subpopulations living within a flood zone are low income persons followed by people under 18 or over 64. The latter groups are priority populations who may have limited mobility in the event of an evacuation. 

```{r tabNFHZApop, echo=FALSE, message=FALSE, warning=FALSE, fig.align = "center", fig.cap="New Hampshire Populations Living within Flood Zones."}
# Show table of pops within flood zones
nh_FloodPops_df %>% 
  mutate(PctFlood = paste0(as.character(round(PctFlood,1)),"%")) %>% 
  arrange(-FloodPop) %>% 
  mutate(Group = recode(Group, "Minority" = "People of Color",
                        "No HS Dip" = "No HS Diploma")) %>%
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','),
                    caption = "New Hampshire Populations Living within Flood Zones",
                    digits = 0, align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header"))
```

Some groups are disproportionately exposed to these flood hazards. While approximately 5.8% of the general population are living within flood zones, households without access to a car, people with disabilities, people in low income households, and other priority populations, exceed the population average living within flood zones. These priority populations are disproportionately exposed to flood risk (see Figure \@ref(fig:chartNFHZApop) below).

```{r chartNFHZApop, echo=FALSE, message=FALSE, warning=FALSE, fig.align = "center", fig.cap="New Hampshire Populations Living within Flood Zones."}
# Create lollipop plot of pops within flood zones
nh_FloodPops_df %>% 
  mutate(Group = recode(Group, "Minority" = "People of Color",
                        "No HS Dip" = "No HS Diploma")) %>%
  ggplot(aes(x = reorder(Group,PctFlood), 
             y = PctFlood)) +
  geom_segment(aes(x = reorder(Group,PctFlood), 
                   xend = reorder(Group,PctFlood),
                   y = nh_FloodPops_df[1,4], yend = PctFlood), 
               color = "skyblue") +
  geom_point(color = "blue", size = 4, alpha = 0.8) +
  coord_flip() + xlab("") + ylab("") + 
  ggtitle("New Hampshire Populations Living within Flood Zones") + theme_light() +
  theme(panel.grid.major.y = element_blank(),
        panel.border = element_blank(),
        axis.ticks.y = element_blank()) +
  geom_text(aes(x = Group, y = PctFlood + 0.2 * sign(PctFlood), 
                label = paste0(round(PctFlood,1),"%")), 
            hjust = 1.8, vjust = -1, size = 3,
            color=rgb(100,100,100, maxColorValue=255)) +
  scale_y_continuous(labels = function(x) paste0(x, "%")) +
  geom_hline(yintercept = nh_FloodPops_df[1,4], linetype = "dashed") +
  geom_text(aes(x = "Over 64", y = 7, label = "Above state\naverage"),
            color = "gray48") +
  geom_text(aes(x = "Over 64", y = 5, label = "Below state\naverage"),
            color = "gray48") +
  expand_limits(y=c(4.5,8))

ggsave("images/NH_FEMA_graph.png")
  # geom_segment(aes(x = "Total Pop", xend = "Total Pop", y = 5.9, yend = 6.6), 
  #              arrow = arrow(length = unit(0.3,"cm")))
```

The maps below (Figures \@ref(fig:mapNFHZAPopConcern1) and \@ref(fig:mapNFHZAPopConcern2)) show households without a car and disabled persons who are disproportionately exposed to flood risk. Maps for other priority populations can be found in Figures \@ref(fig:mapNFHZAPopConcern3) to \@ref(fig:mapNFHZAPopConcern9) and breakdowns by municipality in Tables \@ref(tab:tabNFHZANoCar) to \@ref(tab:tabNFHZAU5) in Appendix B.

```{r mapNFHZAPopConcern, echo=FALSE, message=FALSE, warning=FALSE, fig.align = "center", fig.cap="New Hampshire Populations Living within Flood Zones."}
# Create a dot density map of transit-dependent populations in flood zone
# # create a function to select variable, filter records, and create sample points
# pop2points <- function(sf, x, group){
#   x <- enquo(x)
#   sf %>%
#     select(!!x) %>%
#     filter(!!x >= 10) %>%
#     st_sample(., size = round(.[[as_label(x)]]/10)) %>%
#     st_sf(.) %>%
#     mutate(Group = group)
# }
# 
# LimitEngHH_nfhza_pts <- pop2points(nh_blkgrps_nfhza,NewEng_limit,"Limited English HH")
# 
# NoCarHH_nfhza_pts <- pop2points(nh_tracts_nfhza,NewNoCar,"No Car HH")
# 
# NoHSdip_nfhza_pts <- pop2points(nh_blkgrps_nfhza,NewLths,"No HS Dip")
# 
# LowInc_nfhza_pts <- pop2points(nh_blkgrps_nfhza,NewPov,"Low Income")
# 
# Disabled_nfhza_pts <- pop2points(nh_tracts_nfhza,NewDisabled,"Disabled")
# 
# Minor_nfhza_pts <- pop2points(nh_blkgrps_nfhza,NewMinority,"Minority")
# 
# Under5_nfhza_pts <- pop2points(nh_blkgrps_nfhza,NewUnder5,"Under 5")
# 
# Over64_nfhza_pts <- pop2points(nh_blkgrps_nfhza,NewOver64,"Over 64")
# 
# Under18_nfhza_pts <- pop2points(nh_blkgrps_nfhza,NewUnder18,"Under 18")
# 
# # Bring them together
# nh_nfhza_vulnerable <- rbind(NoCarHH_nfhza_pts,
#                              LimitEngHH_nfhza_pts,
#                              Disabled_nfhza_pts,
#                              LowInc_nfhza_pts,
#                              Minor_nfhza_pts,
#                              Over64_nfhza_pts,
#                              Under5_nfhza_pts,
#                              Under18_nfhza_pts,
#                              NoHSdip_nfhza_pts) %>%
#   # slice(sample(1:n())) %>% # randomise order to avoid bias in plotting order
#   group_by(Group) %>%
#   summarize()
# 
# # clean up pts layers
# rm(list = ls(pattern = "_nfhza_pts"))
# 
# # save time by just loading file
# saveRDS(nh_nfhza_vulnerable, file = "DATA/FEMA/NH/nh_nfhza_vulnerable.Rds")

# load vulnerable point file
nh_nfhza_vulnerable <- readRDS("DATA/FEMA/NH/nh_nfhza_vulnerable.Rds")
```

```{r mapNFHZAPopConcern1, fig.align = "center", fig.cap="Households without a Car Living within Flood Zones."}
# create pop subset
nfhza_sub <- nh_nfhza_vulnerable %>% 
  filter(Group == "No Car HH")

# Map over-represented pops and flood zones
m <- tm_layout(bg.color = "#e6f3f7") +
  tm_shape(nh_blkgrp_2823, unit = "mi") + tm_fill(col = "white") +
  tm_shape(ne_states_sf_cb) + tm_fill(col="white") +
  tm_shape(canada) + tm_fill(col = "white") +
  # tm_shape(ny_state_sf) + tm_fill(col="white") +
  tm_shape(nh_awater_sf_simple) + tm_fill(col = "#e6f3f7") +
  tm_shape(nfhza_sub) +
  tm_dots(col = "darkgoldenrod3",
          labels = "1 dot = 10 people",
          alpha = 0.8) +
  # tm_facets(by = "Group",
  #           ncol = 2,
  #           free.coords = FALSE) +
  tm_shape(nh_nfhza_2823_land_simple) +
  tm_fill(col = "Interval",
          palette = c("dodgerblue", "blueviolet"),
          labels = c("1% AEP (100-year)", "0.2% AEP (500-year)"),
          title = "FEMA Flood Zones",
          alpha = 0.6,
          border.alpha = 0) +
  tm_shape(ne_states_sf_cb) + tm_borders(lwd = 0.6) +
  # tm_shape(ny_state_sf) + tm_borders(lwd = 0.6) +
  tm_shape(nh_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(nh_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_shape(nofemapolys) +
  tm_fill(col = "gray", alpha = 0.4) +
  tm_text("LABEL", alpha = 0.5, size = 0.6) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5,
               position = c("left","TOP")) +
  tm_add_legend(type = "fill", col = "darkgoldenrod3",
                alpha = 0.6,
                border.col = "white", border.alpha = 0,
                labels = "1 dot = 10 people") +
  tm_layout(title = "Households\nwithout a Car\nLiving within\nFlood Zones",
            frame = TRUE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/FEMA/NH/NH_nhfza_vulnerable1.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/FEMA/NH/NH_nhfza_vulnerable1.png")
```

```{r mapNFHZAPopConcern2, fig.align = "center", fig.cap="Disabled Persons Living within Flood Zones."}
# create pop subset
nfhza_sub <- nh_nfhza_vulnerable %>% 
  filter(Group == "Disabled")

# Map over-represented pops and flood zones
m <- tm_layout(bg.color = "#e6f3f7") +
  tm_shape(nh_blkgrp_2823, unit = "mi") + tm_fill(col = "white") +
  tm_shape(ne_states_sf_cb) + tm_fill(col="white") +
  tm_shape(canada) + tm_fill(col = "white") +
  # tm_shape(ny_state_sf) + tm_fill(col="white") +
  tm_shape(nh_awater_sf_simple) + tm_fill(col = "#e6f3f7") +
  tm_shape(nfhza_sub) +
  tm_dots(col = "darkgoldenrod3",
          labels = "1 dot = 10 people",
          alpha = 0.8) +
  # tm_facets(by = "Group",
  #           ncol = 2,
  #           free.coords = FALSE) +
  tm_shape(nh_nfhza_2823_land_simple) +
  tm_fill(col = "Interval",
          palette = c("dodgerblue", "blueviolet"),
          labels = c("1% AEP (100-year)", "0.2% AEP (500-year)"),
          title = "FEMA Flood Zones",
          alpha = 0.6,
          border.alpha = 0) +
  tm_shape(ne_states_sf_cb) + tm_borders(lwd = 0.6) +
  # tm_shape(ny_state_sf) + tm_borders(lwd = 0.6) +
  tm_shape(nh_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(nh_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_shape(nofemapolys) +
  tm_fill(col = "gray", alpha = 0.4) +
  tm_text("LABEL", alpha = 0.5, size = 0.6) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5,
               position = c("left","TOP")) +
  tm_add_legend(type = "fill", col = "darkgoldenrod3",
                alpha = 0.6,
                border.col = "white", border.alpha = 0,
                labels = "1 dot = 10 people") +
  tm_layout(title = "Disabled\nPersons\nLiving within\nFlood Zones",
            frame = TRUE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/FEMA/NH/NH_nhfza_vulnerable2.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/FEMA/NH/NH_nhfza_vulnerable2.png")
```

\pagebreak
## Analysis of Hurricane Evacuation Risk
```{r HEVACdata, include=FALSE}
# read in hurricane evacuation zone layer
nh_hea_sf <- st_read(dsn = "DATA/FEMA/NH",layer = "nh_hevac") %>% 
  mutate(Zone = as.character(Zone)) %>% 
  filter(Zone %in% c("A","B")) %>% 
  st_transform(., crs = 2823) %>% 
  st_make_valid() %>% 
  mutate(Area = st_area(.))

# read in processed nh_blkgrps and nh_tracts
nh_blkgrps_hevac <- st_read(dsn = "DATA/FEMA/NH", 
                            layer = "nh_blkgrps_hevac") %>% 
  left_join(., as.data.frame(nh_blkgrp_2823), by = "GEOID") %>%
  st_transform(., crs = 2823) %>% 
  mutate(NewArea = st_area(.)) %>%
  st_make_valid()

nh_tracts_hevac <- st_read(dsn = "DATA/FEMA/NH", 
                           layer = "nh_tracts_hevac") %>% 
  left_join(., as.data.frame(nh_tracts_2823), by = "GEOID") %>%
  st_transform(., crs = 2823) %>% 
  mutate(NewArea = st_area(.)) %>% 
  st_make_valid()

# Apportion populations based on geographic proportion of intersect
nh_blkgrps_hevac <- nh_blkgrps_hevac %>%
  mutate(Proportion = as.numeric(NewArea/OldArea),
         NewPop = totalpopE*Proportion,
         NewMinority = minorityE*Proportion,
         NewUnder5 = under5E*Proportion,
         NewOver64 = over64E*Proportion,
         NewUnder18 = under18E*Proportion,
         NewEng_limit = eng_limitE*Proportion,
         NewPov = num2povE*Proportion,
         NewLths = lthsE*Proportion)

nh_tracts_hevac <- nh_tracts_hevac %>%
  mutate(Proportion = as.numeric(NewArea/OldArea),
         NewDisabled = disabledOver18E*Proportion,
         NewNoCar = HHnoCarE*Proportion)

# Compute total populations within hurricane evac zones
nh_hevac_blkgrp_df <- nh_blkgrps_hevac %>%
  as.data.frame() %>%
  summarize(`Total Pop` = as.integer(sum(NewPop)),
            Minority = as.integer(sum(NewMinority)),
            `Under 5` = as.integer(sum(NewUnder5)),
            `Over 64` = as.integer(sum(NewOver64)),
            `Under 18` = as.integer(sum(NewUnder18)),
            `Limited English HH` = as.integer(sum(NewEng_limit)),
            `Low Income` = as.integer(sum(NewPov)),
            `No HS Dip` = as.integer(sum(NewLths))) %>%
  gather(key = Group, value = HevacPop)

nh_hevac_tracts_df <- nh_tracts_hevac %>%
  as.data.frame() %>%
  summarize(`Disabled` = as.integer(sum(NewDisabled)),
            `No Car HH` = as.integer(sum(NewNoCar))) %>%
  gather(key = Group, value = HevacPop)

# Compute total tract populations within the state for same groups
nh_tranh_hevac_pops_df <- nh_tracts_2823 %>%
  as.data.frame() %>%
  summarize(`Disabled` = sum(disabledOver18E),
            `No Car HH` = sum(HHnoCarE)) %>%
  gather(key = Group, value = NHPop) %>%
  left_join(.,nh_hevac_tracts_df, by = "Group")

# Compute populations for state,  and join with hurricane evac pops
nh_HevacPops_df <- nh_blkgrp_2823 %>%
  as.data.frame() %>%
  summarize(`Total Pop` = sum(totalpopE),
            Minority = sum(minorityE),
            `Under 5` = sum(under5E),
            `Over 64` = sum(over64E),
            `Under 18` = sum(under18E),
            `Limited English HH` = sum(eng_limitE),
            `Low Income` = sum(num2povE),
            `No HS Dip` = sum(lthsE)) %>%
  gather(key = Group, value = NHPop) %>%
  left_join(., nh_hevac_blkgrp_df, by = "Group") %>%
  rbind(.,nh_tranh_hevac_pops_df) %>%
  mutate(PctHevac = HevacPop/NHPop*100)

# save data for later analysis
save(nh_blkgrps_hevac, nh_tracts_hevac, 
     nh_HevacPops_df,
     file = "DATA/FEMA/NH/hevac_census.Rds")
```

New Hampshire is subject to hurricane risk. Since 1900, New Hampshire has been struck by hurricanes twice. The most recent hurricane to hit New Hampshire directly was [Hurricane Bob in 1991](https://www.mass.gov/service-details/the-worst-New Hampshire-hurricanes-of-the-20th-century),[^HurricaneBob] a Category 2 storm when it struck the  south shore of New Hampshire. The latter caused over \$2.3 million in damages in New Hampshire alone and left 20% of the state without power. Hurricane Sandy in 2012 did not hit New Hampshire directly, but nevertheless resulted in flooding of New Hampshire coastal communities that resulted in over \$1.6 million in damanges. In general, New Hampshire is subject to a [hurricane return period of approximately 30 years](https://www.nhc.noaa.gov/climo/#returns).[^HurricaneReturn] 

Hurricane intensity is ranked from 1 to 5 on the Saffir-Simpson scale, which is based on sustained wind speeds. However the most destructive and deadly aspects of hurricanes are due to storm surge - ocean water pushed inland by the high winds. This coastal storm surge may be combined with intense rainfall which can exacerbate flooding. As a result, coastal evacuation may be necessary as a hurricane approaches.

The analysis of hurricane evacuation risk presented here is based on [storm surge modeling by the National Weather Service's National Hurricane Center](https://www.nhc.noaa.gov/nationalsurge/) using the hydrodynamic Sea, Lake, and Overland Surges from Hurricanes (SLOSH) model to simulate storm surge from tropical cyclones. These data products are provided to federal, state, and local government agencies to assist in planning, risk assessment studies, and operational decision-making.[^SurgeModel] Population data comes from the American Community Survey (ACS) 2018 5-year estimates at the census tract and block group levels. 

A significant proportion of New Hampshire's population lives near the coast and is therefore within areas subject to inundation due to hurricane storm surge. Figure \@ref(fig:mapHEVACPop) below shows the population distribution and areas subject to hurricane inundation. Zone “A” is recommended to be evacuated prior to an expected category 1 or 2 hurricane.  Zone “B” is recommended to be evacuated prior to an expected category 3 or 4 hurricane.

```{r mapHEVACPop, echo=FALSE, message=FALSE, warning=FALSE, fig.align = "center", fig.cap="Hurricane flood zones across New Hampshire."}
# Create a simplified version of hea polygons to speed up mapping
empty_geo <- st_is_empty(nh_hea_sf)
nh_hea_sf <- nh_hea_sf[!empty_geo,]

nh_hea_sf_simple <- nh_hea_sf %>% 
  select(Zone) %>% 
  # st_buffer(., dist = 0) %>%
  # st_set_precision(1e6) %>%
  # st_cast(., "MULTIPOLYGON") %>% # homogenize type
  # st_cast(., "POLYGON") %>% 
  st_simplify(., dTolerance = 100) %>% # reduce vertices
  group_by(Zone) %>% 
  summarize() %>% 
  st_make_valid() %>%
  mutate(Area = st_area(.))

# create tailored labels for close-up of hurricane areas
I95roadSegment2 <- nh_highways %>% 
  filter(LINEARID == "1105569136123")

# add labels for towns for close-up of hurricane areas
nh_towns_sf_pts2 <- county_subdivisions(state = "NH", cb = TRUE) %>% 
  filter(NAME %in% c("Seabrook",
                     "Hampton",
                     "Hampton Falls",
                     "North Hampton",
                     "Rye",
                     "Portsmouth",
                     "Newington",
                     "Greenland",
                     "Durham",
                     "Dover")) %>% 
  st_transform(., crs = 2823) %>% 
  st_centroid(of_largest_polygon = TRUE)

# clip water to area of focus for faster mapping
coastalWater <- nh_counties_sf %>% 
  filter(NAME %in% c("Rockingham","Strafford")) %>% 
  crop_shape(nh_awater_sf,.,polygon = TRUE)

# Map totalpop and hurricane evacuation zones
m <- tm_layout(bg.color = "#e6f3f7") +
  tm_shape(nh_hea_sf_simple, unit = "mi") + tm_fill(col = "white") +
  tm_shape(nh_blkgrp_2823, unit = "mi") + tm_fill(col = "white") +
  tm_shape(ne_states_sf_cb) + tm_fill(col="white") +
  # tm_shape(ny_state_sf) + tm_fill(col="white") +
  tm_shape(coastalWater) + tm_fill(col = "#e6f3f7") + 
  tm_shape(ne_states_sf_cb) + tm_borders() +
  # tm_shape(ny_state_sf) + tm_borders() +
  tm_shape(nh_hea_sf) + 
  tm_fill(col = "Zone", 
          palette = c("deeppink", "purple3"), 
          labels = c("A: Category 1 - 2", 
                     "B: Category 3 - 4"), 
          title = "Inundation Zone and\nHurricane Category",
          alpha = 0.6,
          border.alpha = 0) +
  # tm_shape(nh_totalpop_pts) + 
  # tm_dots(col = "springgreen3",
  #         labels = "1 dot = 500 people",
  #         alpha = 0.8) +
  tm_shape(ne_states_sf_cb) + tm_borders(lwd = 0.8) +
  # tm_shape(ny_state_sf) + tm_borders(lwd = 0.8) +
  tm_shape(nh_highways) + tm_lines(col = "seashell4", lwd = 1) +
  # tm_shape(nh_highways2nd) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(I89roadSegment) + 
  tm_symbols(shape = I89, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment) + 
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) + 
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(nh_towns_sf_pts2) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c("left","TOP")) +
  # tm_add_legend(type = "fill", col = "springgreen3",
  #               alpha = 0.6,
  #               border.col = "white", border.alpha = 0,
  #               labels = "1 dot = 500 people",
  #               title = "Total Population") +
  tm_layout(title = "Hurricane\nInundation Zones", 
            frame = TRUE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/FEMA/NH/NH_hevac.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/FEMA/NH/NH_hevac.png")
```

Approximately 0.3% of New Hampshire's land area falls within areas subject to hurricane inundation. The breakdown by hurricane inundation zone is presented in Table \@ref(tab:tabHEVACarea) below.

```{r tabHEVACarea, echo=FALSE, message=FALSE, warning=FALSE, fig.align = "center", fig.cap="New Hampshire Populations within Hurricane Inundation Zones."}
# Total and percentage of land area within hurricane evacuation zones
nh_hea_sf %>% 
  as.data.frame() %>% 
  group_by(Zone) %>% 
  summarize(SqKm = round(as.numeric(sum(Area)/10^6),1), 
            SqMi = round(as.numeric(SqKm/2.59),1),
            PctArea = paste0(as.character(round(as.numeric(sum(Area)/nh_area*100),1)),"%")) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','),
                    caption = "New Hampshire Land Area within Hurricane Inundation Zones",
                    digits = 0, align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header"))
```

These percentages belie the proportion of the population in the state that are living within these inundation zones. Exposure varies significantly by population subgroup. Table \@ref(tab:tabHEVACpop) below shows the total and percentages of the general population and various subgroups living within hurricane inundation zones. Aside from the total population, the largest subpopulation living within an inundation zone are individuals over age 64 and low income individuals. The latter groups are priority populations who may have limited mobility in the event of an evacuation. 

```{r tabHEVACpop, echo=FALSE, message=FALSE, warning=FALSE, fig.align = "center", fig.cap="New Hampshire Populations Living within Hurricane Evacuation Zones."}
# Show table of pops within hurricane evacuation zones
nh_HevacPops_df %>% 
  mutate(PctHevac = paste0(as.character(round(PctHevac,1)),"%")) %>% 
  arrange(-HevacPop) %>% 
  mutate(Group = recode(Group, "Minority" = "People of Color",
                        "No HS Dip" = "No HS Diploma")) %>%
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','),
                    caption = "New Hampshire Populations Living within Hurricane Inundation Zones",
                    digits = 0, align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header"))
```

Some groups are disproportionately exposed to these hurricane evacuation hazards. While approximately 1.1% of the general population are living within hurricane inundation zones, individuals over age 64 and households without access to a car exceed the population average living within hurricane inundation zones. These priority populations are disproportionately exposed to the risk of hurricane inundation (see Figures \@ref(fig:chartHEVACpop) and \@ref(fig:mapHEVACPopConcern1) below).

```{r chartHEVACpop, echo=FALSE, message=FALSE, warning=FALSE, fig.align = "center", fig.cap="New Hampshire Populations Living within Hurricane Evacuation Zones."}
# Create lollipop plot of pops within flood zones
nh_HevacPops_df %>% 
  mutate(Group = recode(Group, "Minority" = "People of Color",
                        "No HS Dip" = "No HS Diploma")) %>%
  ggplot(aes(x = reorder(Group,PctHevac), 
             y = PctHevac)) +
  geom_segment(aes(x = reorder(Group,PctHevac), 
                   xend = reorder(Group,PctHevac),
                   y = nh_HevacPops_df[1,4], yend = PctHevac), 
               color = "skyblue") +
  geom_point(color = "blue", size = 4, alpha = 0.8) +
  coord_flip() + xlab("") + ylab("") + 
  ggtitle("New Hampshire Populations Living within Hurricane\nInundation Zones") + theme_light() +
  theme(panel.grid.major.y = element_blank(),
        panel.border = element_blank(),
        axis.ticks.y = element_blank()) +
  geom_text(aes(x = Group, y = PctHevac + 0.2 * sign(PctHevac), 
                label = paste0(round(PctHevac,1),"%")), 
            hjust = 1.7, vjust = -0.7, size = 3,
            color=rgb(100,100,100, maxColorValue=255)) +
  scale_y_continuous(labels = function(x) paste0(x, "%")) +
  geom_hline(yintercept = nh_HevacPops_df[1,4], linetype = "dashed") +
  geom_text(aes(x = "Total Pop", y = 1.5, label = "Above state\naverage"),
            color = "gray48") +
  geom_text(aes(x = "Total Pop", y = .8, label = "Below state\naverage"),
            color = "gray48")

ggsave("images/NH_HEVAC_graph.png")
  # geom_segment(aes(x = "Total Pop", xend = "Total Pop", y = 1.3, yend = 2), 
  #              arrow = arrow(length = unit(0.3,"cm")))
```

The maps below (Figures \@ref(fig:mapHEVACPopConcern1) and \@ref(fig:mapHEVACPopConcern2)) show persons over 64 and households without a car who are disproportionately exposed to hurricane evacuation risk. Maps for other priority populations can be found in Figures \@ref(fig:mapHEVACPopConcern3) to \@ref(fig:mapHEVACPopConcern7) and breakdowns by municipality in Tables \@ref(tab:tabHEVACO64) to \@ref(tab:tabHEVACLEH) in Appendix B.

```{r mapHEVACPopConcern, echo=FALSE, message=FALSE, warning=FALSE, fig.align = "center", fig.cap="New Hampshire Populations Living within Hurricane Evacuation Zones."}
# # Create a dot density map of over-represented populations in hurricane evacuation zones
# 
# # LimitEngHH_hevac_pts <- pop2points(nh_blkgrps_hevac,NewEng_limit,"Limited English HH") # no blkgrps >= 10
# 
# NoCarHH_hevac_pts <- pop2points(nh_tracts_hevac,NewNoCar,"No Car HH")
# 
# NoHSdip_hevac_pts <- pop2points(nh_blkgrps_hevac,NewLths,"No HS Dip")
# 
# LowInc_hevac_pts <- pop2points(nh_blkgrps_hevac,NewPov,"Low Income")
# 
# Disabled_hevac_pts <- pop2points(nh_tracts_hevac,NewDisabled,"Disabled")
# 
# Minor_hevac_pts <- pop2points(nh_blkgrps_hevac,NewMinority,"Minority")
# 
# Under5_hevac_pts <- pop2points(nh_blkgrps_hevac,NewUnder5,"Under 5")
# 
# Over64_hevac_pts <- pop2points(nh_blkgrps_hevac,NewOver64,"Over 64")
# 
# Under18_hevac_pts <- pop2points(nh_blkgrps_hevac,NewUnder18,"Under 18")
# 
# # Bring them together
# nh_hevac_vulnerable <- rbind(NoCarHH_hevac_pts,
#                              Disabled_hevac_pts,
#                              LowInc_hevac_pts,
#                              Minor_hevac_pts,
#                              Over64_hevac_pts,
#                              Under5_hevac_pts,
#                              Under18_hevac_pts,
#                              NoHSdip_hevac_pts) %>%
#   # slice(sample(1:n())) %>% # randomise order to avoid bias in plotting order
#   group_by(Group) %>%
#   summarize()
# 
# # clean up pts layers
# rm(list = ls(pattern = "_hevac_pts"))
# 
# # save time by just loading file
# saveRDS(nh_hevac_vulnerable, file = "DATA/FEMA/NH/nh_hevac_vulnerable.Rds")

# load vulnerable point file
nh_hevac_vulnerable <- readRDS("DATA/FEMA/NH/nh_hevac_vulnerable.Rds")
```

```{r mapHEVACPopConcern1, fig.align = "center", fig.cap="Persons Over 64 Living within Hurricane Inundation Zones."}
# Create data subset
hevac_sub <- nh_hevac_vulnerable %>% 
  filter(Group == "Over 64")

# Map over-represented pops and hurricane evacuation zones
m2 <- tm_layout(bg.color = "#e6f3f7") +
  tm_shape(nh_hea_sf_simple, unit = "mi") + tm_fill(col = "white") +
  tm_shape(nh_blkgrp_2823, unit = "mi") + tm_fill(col = "white") +
  tm_shape(ne_states_sf_cb) + tm_fill(col="white") +
  # tm_shape(ny_state_sf) + tm_fill(col="white") +
  tm_shape(coastalWater) + tm_fill(col = "#e6f3f7") + 
  tm_shape(ne_states_sf_cb) + tm_borders() +
  # tm_shape(ny_state_sf) + tm_borders() +
  tm_shape(nh_hea_sf) + 
  tm_fill(col = "Zone", 
          palette = c("deeppink", "purple3"), 
          labels = c("A: Category 1 - 2", 
                     "B: Category 3 - 4"), 
          title = "Inundation Zone and\nHurricane Category",
          alpha = 0.6,
          border.alpha = 0) +
  tm_shape(nh_hevac_vulnerable) +
  tm_dots(col = "springgreen3",
          labels = "1 dot = 10 people",
          alpha = 1) +
  # tm_facets(by = "Group",
  #           ncol = 2,
  #           free.coords = FALSE) +
  tm_shape(ne_states_sf_cb) + tm_borders(lwd = 0.2) +
  # tm_shape(ny_state_sf) + tm_borders(lwd = 0.8) +
  tm_shape(nh_highways) + tm_lines(col = "seashell4", lwd = 1) +
  # tm_shape(nh_highways2nd) + tm_lines(col = "seashell4", lwd = 1.5) +
  # tm_shape(I89roadSegment) + 
  # tm_symbols(shape = I89, border.lwd = NA, size = 0.2) +
  # tm_shape(I91roadSegment) + 
  # tm_symbols(shape = I91, border.lwd = NA, size = 0.2) +
  # tm_shape(I93roadSegment) +
  # tm_symbols(shape = I93, border.lwd = NA, size = 0.2) +
  tm_shape(I95roadSegment2) + 
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(nh_towns_sf_pts2) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c("left","TOP")) +
  tm_add_legend(type = "fill", col = "springgreen3",
                alpha = 0.6,
                border.col = "white", border.alpha = 0,
                labels = "1 dot = 10 people") +
  tm_layout(title = "Persons\nOver 64\nLiving within\nHurricane\nInundation Zones", frame = TRUE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m2, "DATA/FEMA/NH/NH_hevac_vulnerable1.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/FEMA/NH/NH_hevac_vulnerable1.png")
```

```{r mapHEVACPopConcern2, fig.align = "center", fig.cap="Households without a Car Living within Hurricane Inundation Zones."}
# Create data subset
hevac_sub <- nh_hevac_vulnerable %>% 
  filter(Group == "No Car HH")

# Map over-represented pops and hurricane evacuation zones
m2 <- tm_layout(bg.color = "#e6f3f7") +
  tm_shape(nh_hea_sf_simple, unit = "mi") + tm_fill(col = "white") +
  tm_shape(nh_blkgrp_2823, unit = "mi") + tm_fill(col = "white") +
  tm_shape(ne_states_sf_cb) + tm_fill(col="white") +
  # tm_shape(ny_state_sf) + tm_fill(col="white") +
  tm_shape(coastalWater) + tm_fill(col = "#e6f3f7") + 
  tm_shape(ne_states_sf_cb) + tm_borders() +
  # tm_shape(ny_state_sf) + tm_borders() +
  tm_shape(nh_hea_sf) + 
  tm_fill(col = "Zone", 
          palette = c("deeppink", "purple3"), 
          labels = c("A: Category 1 - 2", 
                     "B: Category 3 - 4"), 
          title = "Inundation Zone and\nHurricane Category",
          alpha = 0.6,
          border.alpha = 0) +
  tm_shape(nh_hevac_vulnerable) +
  tm_dots(col = "springgreen3",
          labels = "1 dot = 10 people",
          alpha = 1) +
  # tm_facets(by = "Group",
  #           ncol = 2,
  #           free.coords = FALSE) +
  tm_shape(ne_states_sf_cb) + tm_borders(lwd = 0.2) +
  # tm_shape(ny_state_sf) + tm_borders(lwd = 0.8) +
  tm_shape(nh_highways) + tm_lines(col = "seashell4", lwd = 1) +
  # tm_shape(nh_highways2nd) + tm_lines(col = "seashell4", lwd = 1.5) +
  # tm_shape(I89roadSegment) + 
  # tm_symbols(shape = I89, border.lwd = NA, size = 0.2) +
  # tm_shape(I91roadSegment) + 
  # tm_symbols(shape = I91, border.lwd = NA, size = 0.2) +
  # tm_shape(I93roadSegment) +
  # tm_symbols(shape = I93, border.lwd = NA, size = 0.2) +
  tm_shape(I95roadSegment2) + 
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(nh_towns_sf_pts2) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c("left","TOP")) +
  tm_add_legend(type = "fill", col = "springgreen3",
                alpha = 0.6,
                border.col = "white", border.alpha = 0,
                labels = "1 dot = 10 people") +
  tm_layout(title = "Households\nwithout a Car\nLiving within\nHurricane\nInundation Zones", frame = TRUE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m2, "DATA/FEMA/NH/NH_hevac_vulnerable2.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/FEMA/NH/NH_hevac_vulnerable2.png")
```


[^NFHL]: See FEMA National Flood Hazard Layer (NFHL). https://www.fema.gov/national-flood-hazard-layer-nfhl.

[^NHCsloshMaps]: The National Hurricane Center provides static and interactive maps of storm surge flooding vulnerability across the country based on SLOSH modeling. See NOAA National Hurricane Center and Central Pacific Hurricane Center. National Storm Surge Hazard Maps - Version 2. https://www.nhc.noaa.gov/nationalsurge/

[^HurricaneBob]: See Mass.gov. "The Worst Massachusetts Hurricanes of the 20th Century." https://www.mass.gov/service-details/the-worst-massachusetts-hurricanes-of-the-20th-century. See also The Northeast States Emergency Consortium (NESEC). "Hurricanes." http://nesec.org/hurricanes/.

[^BobMaine]: See Maine Emergency Management Agency. "Hurricanes." https://www.maine.gov/mema/hazards/natural-hazards/hurricanes.

[^BobRI]: See https://www.nhc.noaa.gov/archive/storm_wallets/atlantic/atl1991/bob/preloc/fema02.gif

[^IreneVT]: See Wislon Ring and Jennifer Peltz. 2011. "Vermont Towns Battle Historic Floods From Irene As Death Toll Tops 40." *HuffPost*. Oct 29. https://www.huffpost.com/entry/vermont-towns-battle-flood-irene_n_941207.

[^SandyCT]: See Ken Dixon. 2012. "Sandy storm damage tops $360M in state." *ctpost*. Nov 14. https://www.ctpost.com/news/article/Sandy-storm-damage-tops-360M-in-state-4037538.php.

[^SandyME]: See Seth Koenig. 2012. "Sandy knocks out power for thousands of Mainers; flood watch inland." *BDN Portland*. Oct 29. https://bangordailynews.com/2012/10/29/news/portland/hurricane-sandy-remains-unpredictable-for-maine-especially-islanders/.

[^HurricaneReturn]: See NOAA National Hurricane Center and Central Pacific Hurricane Center. "Tropical Cyclone Climatology - U.S. Hurricane Return Periods." https://www.nhc.noaa.gov/climo/#returns.

[^SurgeModel]: See NOAA National Hurricane Center and Central Pacific Hurricane Center. "National Storm Surge Hazard Maps - Version 2." https://www.nhc.noaa.gov/nationalsurge/.

[^Bruggers2018]: James Bruggers. Nov 1, 2018. "FEMA Flood Maps Ignore Climate Change, and Homeowners Are Paying the Price." *Insideclimatenews.* https://insideclimatenews.org/news/01112018/fema-flood-map-climate-change-hurricane-mexico-beach-florida-sea-level-rise 

[^Dasymetric]: Areal apportionment of populations for mapping and analysis, also known as dasymetric mapping, is commonly used method to improve risk assessment of flooding. See for example Yi Qiang. 2019. "Disparities of population exposed to flood hazards in the United States." *Journal of Environmental Management* 232:295-304; Yi Qiang, Nina S. N. Lam, Heng Cai, and Lei Zou. 2017. "Changes in Exposure to Flood Hazards in the United States." *Annals of the American Association of Geographers* 107(6):1332-1350. https://doi.org/10.1080/24694452.2017.1320214

[^NLCD2016]: See MRLC NLCD 2016 Land Cover (CONUS). https://www.mrlc.gov/data/nlcd-2016-land-cover-conus.

[^NLCDdetails]: See Collin Homer, Jon Dewitz, Suming Jin, George Xian, Catherine Costello, Patrick Danielson, Leila Gass, Michelle Funk, James Wickham, Stephen Stehman, Roger Auch, and Kurt Riitters. 2020. "Conterminous United States land cover change patterns 2001–2016 from the 2016 National Land Cover Database." *ISPRS Journal of Photogrammetry and Remote Sensing* 162: 184-199. https://doi.org/10.1016/j.isprsjprs.2020.02.019.


\pagebreak
# Appendix A: Data and Methodology {-}

The analyses presented here are based on data from five sources:

* FEMA National Flood Hazard Layer (NFHL)
* MRLC National Land Cover Database (NLCD)
* State Hurricane Evacuation Zones (CT,MA,RI)
* NOAA Sea, Lake, and Overland Surges from Hurricanes (SLOSH) model
* American Community Survey 5-year Estimates
  + Population demographics

## FEMA National Flood Hazard Layer (NFHL) {-}

As part of the National Flood Insurance Program, the U.S. Federal Emergency Management Agency (FEMA) produces maps of flood risk across the nation. These flood risk maps are based on . 

The National Flood Hazard Layer (NFHL) is a digital geospatial database that contains current effective flood hazard mapping data from the U.S. Federal Emergency Management Agnecy's (FEMA) National Flood Insurance Program (NFIP). NFHL data is available for download at https://msc.fema.gov/portal/advanceSearch. This geodatabase identifies areas subject to various levels of flood risk as determined by FEMA, consistent with Flood Insurance Rate Maps (FIRMs), based upon hydrodynamic models as well as historic flooding. 

This analysis focused on identifying populations at risk from floods with a 1% or 0.2% annual exceedance probability. The 1% annual exceedance probability flood, sometimes referred to as a '100-year' flood, is a flood with a 1% chance of being equaled or exceeded in a given year. Areas subject to the 1% flood risk are designated as Special Flood Hazard Areas. Properties within these these flood risk zones are required to carry flood insurance in order to participate in the National Flood Insurance Program. The 0.2% exceedance probability flood, sometimes referred to as a '500-year' flood, is a flood with a 0.2% chance of being equaled or exceeded in a given year. There is no federal regulatory requirement for properties within areas subject to 0.2% annual exceedance probability flood risk. While lower risk, these areas are frequently identified as areas subject to flooding under extreme or worst case scenarios. Note that FEMA flood risk designations do not take into consideration increasing or projected flood risk due to sea level rise, more extreme rainfall events, or other climate change-enhanced risks.[^Bruggers2018] 

In order to identify populations at risk from floods with a 1% or 0.2% annual exceedance probability, NFHL polygons pertaining to the 1% flood risk (zones A, AE, AH, AO, and VE) and the 0.2% flood risk (zone X), and excluding water or areas of "minimal flood hazard," were extracted and intersected with the developed portions of Census Block Groups. Developed portions of Census Block Groups were identified based upon the National Land Cover Database (NLCD), using a process of areal apportionment (see below). The population at risk from flooding was calculated as the product of the  areal proportion of the intersecting flood and developed Block Group polygons:

*Population at risk = Proportion of developed Block Group Intersection x Population of developed Block Group*

For example, if 10% of the developed area of a Census Block Group intersected/overlapped with a flood polygon, it was assumed that 10% of the population is exposed to that flood risk. Assuming a population of 100 people in the developed portion of the Block Group, this would mean *100 x .10 = 10* people would be subject to flood risk. FEMA flood risk modeling was done in R.


## MRLC National Land Cover Database (NLCD) {-}

The National Land Cover Database (NLCD) provides nationwide data on land cover and land cover change at a 30m resolution with a 16-class legend based on a modified Anderson Level II classification system:

```{r NLCDtab}
gridCodes <- c(11,12,21,22,23,24,31,41,42,43,51,52,71,72,73,74,81,82,90,95)

landUse <- c("Open Water", "Perennial Ice/Snow", "Developed, Open Space", "Developed, Low Intensity", "Developed, Medium Intensity", "Developed, High Intensity", "Barren Land (Rock/Sand/Clay)", "Deciduous Forest", "Evergreen Forest", "Mixed Forest", "Dwarf Scrub", "Shrub/Scrub", "Grassland/Herbaceous", "Sedge/Herbaceous", "Lichens", "Moss", "Pasture/Hay", "Cultivated Crops", "Woody Wetlands", "Emergent Herbaceous Wetlands")

data.frame(gridCodes,landUse) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    col.names = c("GRID Code","Land Use Description"),
                    caption = "NLCD 2016 Land Cover Codes and Classifications", align = "l") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header"))
```

The NLCD is produced by the US Geological Survey (USGS) in partnership with several federal agencies as part of the Multi-Resolution Land Characteristics Consortium (MRLC).[^NLCD2016] NLCD identified land cover types based on semi-automated land use calssification algorithms applied to LANDSAT satellite imagery and supplemented with ancillary data.[^NLCDdetails] NLCD is distributed as a raster GeoTIFF at https://www.mrlc.gov/data. 

NLCD data was used to identify developed areas of Census Block Groups in order to better estimate potential flood risk in areas where people are most likely to reside.[^Dasymetric] The NLCD 2016 was downloaded and converted to a vector polygon layer and clipped to New England states. Undeveloped polygon areas were extracted (all GRID codes excluding 22 - 24; see Table \@ref(tab:NLCDtab)). These undeveloped areas of the NLCD polygon layer were then used to erase overlapping areas of Census Block Groups. The latter step served to eliminate areas of Census Block Groups that were classified as undeveloped and thus unlikely to be occupied by residences. The remaining developed areas of Census Block Groups were assumed to contain all populations assigned to that respective Block Group. Processing of NLCD and developed Census Block Groups was done in ArcGIS Desktop 10.7. These developed Census Block Group polygons were used for identifying populations at risk from flooding or hurricane storm surge.

## State Hurricane Evacuation Zones (CT,MA,RI) {-}

Connecticut, Massachusetts, and Rhode Island provide maps and geospatial layers identifying evacuation zones in their respective states in the event of a hurricane. These hurricane evaucation zone maps were developed in cooperation between state emergency management agencies and the U.S. Army Corps of Enginneers. Hurricane storm surge risk was based on the SLOSH (Sea, Lake, and Overland Surges from Hurricanes) model, which is a numerical model used by the National Weather Service to compute storm surge. Hurricane evaucation zones generally encompass an area larger than the specific areas subject to inundation. 
Connecticut Hurricane Evacuation Zones: https://data.ct.gov/Public-Safety/Hurricane-Evacuation-Zones-2014-/kdsj-af66
Massachusetts Hurricane Evacuation Zones: https://www.nae.usace.army.mil/Missions/Projects-Topics/Massachusetts-Hurricane-Studies/
Rhode Island Hurricane Evacuation Zones: http://www.riema.ri.gov/resources/citizens/prepare/evacuation.php 

In order to identify populations at risk from hurricane evacuation, hurricane evacuation zone polygons were intersected with the developed portions of Census Block Groups. Developed portions of Census Block Groups were identified based upon the National Land Cover Database (NLCD), using a process of areal apportionment (see above). The population at risk from hurricane evacuation risk was calculated as the product of the  areal proportion of the intersecting evacuation zone and developed Block Group polygons:

*Population at risk = Proportion of developed Block Group Intersection x Population of developed Block Group*

For example, if 10% of the developed area of a Census Block Group intersected/overlapped with a hurricane evacuation zone, it was assumed that 10% of the population is exposed to that evacuation risk. Assuming a population of 100 people in the developed portion of the Block Group, this would mean *100 x .10 = 10* people would be subject to hurricane evacuation risk. 

## NOAA Sea, Lake, and Overland Surges from Hurricanes (SLOSH) model (ME,NH) {-}

The National Oceanic and Atmospheric Administration (NOAA), specifically the National Weather Service's (NWS) National Hurricane Center (NHC), utilizes the hydrodynamic Sea, Lake, and Overland Surges from Hurricanes (SLOSH) model to simulate storm surge from tropical cyclones, including hurricanes. SLOSH models of coastal storm surge are based on tens of thousands of climatology-based hypothetical tropical cyclone simulations.[^SurgeModel]  The composite results of these simulations are used to assess and visualize storm surge risk under varying conditions. Maps and estimates of risk are provided to federal, state, and local partners to assist in a range of planning processes, risk assessment studies, and operational decision-making. SLOSH models and maps are the basis of state-level hurricane evacuation zone designations. However, only three of the six New England states (Connecticut, Massachusetts, and Rhode Island) have developed hurricane evacuation zone maps. Maine and New Hampshire, which are subject to hurricane risk, have not developed, or do not publicly offer, hurricane evacuation zone maps. In order to model hurricane evacuation risk in those states, SLOSH inundation zones were used instead. The use of storm surge inundation risk, rather than hurricane evacuation zones, likely underestimates the population that would be subject to evacuation in the event of an approaching hurricane.

In order to identify populations at risk from hurricane evacuation in Maine and New Hampshire, SLOSH surge hazard maps were downloaded from the National Hurricane Center in the form of GeoTIFF raster files.[^SurgeModel] Specifically, scenarios for Category 2 and Category 4 hurricanes at high tide were used. The raster files were clipped to the land areas of Maine and New Hampshire and converted to a vector polygon layer and merged. Hurricane storm surge polygons were intersected with the developed portions of Census Block Groups. Developed portions of Census Block Groups were identified based upon the National Land Cover Database (NLCD), using a process of areal apportionment (see above). The population at risk from hurricane storm surge risk was calculated as the product of the areal proportion of the intersecting storm surge polygons and developed Block Group polygons:

*Population at risk = Proportion of developed Block Group Intersection x Population of developed Block Group*

For example, if 10% of the developed area of a Census Block Group intersected/overlapped with a hurricane storm surge polygon, it was assumed that 10% of the population is exposed to that flood risk. Assuming a population of 100 people in the developed portion of the Block Group, this would mean *100 x .10 = 10* people would be subject to hurricane storm surge risk. Processing of SLOSH GeoTIFF files and polygon intersection were performed in ArcGIS 10.7. 


## American Community Survey 5-year Estimates {-}

The American Community Survey (ACS) is an ongoing survey by the U.S. Census Bureau that provides information on a yearly basis about the U.S. and its people. The ACS provides detailed information on economic, housing, and demographic characteristics about the population that are not captured by the decennial Census. 

The ACS provides greater demographic detail and temporal resolution than the decennial Census, but its geographic resolution is more limited. While the decennial Census is based on an enumeration (i.e., a total count) of everyone in the U.S. once every decade, the ACS is based on a statistical sample of approximately 3.5 million addresses across the country each year. As a result of this sampling approach, the ACS estimates must be pooled, or combined, across multiple years in order to provide reliable estimates for smaller areas (i.e. areas with less than 20,000 people), such as at the Tract or Block Group levels. While it is possible to know the number of low income households across the U.S. annually, one may only know this about a Tract or Block Group based on 5-year estimates. Since 2010, the ACS has published 5-year data (beginning with 2005–2009 estimates) for all geographic areas down to the census Tract and Block Group levels. For more detail on the ACS, see Understanding and Using American Community Survey Data: What All Data Users Need to Know at https://www.census.gov/programs-surveys/acs/guidance/handbooks/general.html.

All data was analyzed or aggregated geographically by Census Tract and Block Group. A Census Tract is a small, relatively permanent statistical subdivision of a county that contains between 1,200 and 8,000 people. The entire area of a county is covered by Census Tracts, just as the entire area of a state is covered by counties or county equivalents. Census Tracts range in areal size depending on the population density; areal sizes are smaller in denser areas and larger in less densely populated areas. Census Block Groups are subdivisions of Census Tracts that contain between 600 and 3,000 people. Like Tracts, Block Groups range in areal size depending on the population density of the area. Block Groups are the smallest geographic unit at which detailed demographic and household data from the American Community Survey is made available by the U.S. Census Bureau. 

For the purposes of this analysis ACS 5-year estimates for the period 2014 - 2018 for Census Tracts and Block Groups in New England, as well as their associated TIGER/Line spatial files, were downloaded from the Census Bureau via API using the `tigris` package in R. Demographic variables consistent with those used by the EPA in EJSCREEN were chosen, as well as environmental justice population thresholds used by states where available. Table \@ref(tab:ACSvariables) below lists the demographic variables that were downloaded directly or computed from ACS variables:

```{r ACSvariables}
ACSvDF <- data.frame(Variable = c("Total Population",
                                  "People of Color",
                                  "Low Income",
                                  "Limited English Household",
                                  "Less than High School Education",
                                  "Under 5",
                                  "Under 18",
                                  "Over 64",
                                  "Disabled",
                                  "No Car Household",
                                  "CT Income",
                                  "MA Income",
                                  "MA Limited English Household",
                                  "MA POC",
                                  "RI Income",
                                  "RI POC"),
                     Description = c("Total population",
                                     "Persons of Hispanic or Latino origin or persons who are not White",
                                     "People in households where the household income is less than or equal to twice the federal poverty level",
                                     "People in households where all adults speak English less than \"very well\"",
                                     "Adults 25 years+ with less than a high school education",
                                     "Persons under 5 years of age",
                                     "Persons under 18 years of age",
                                     "Persons over 64 years of age",
                                     "Persons 18 years+ with a disability",
                                     "Households with no vehicle available",
                                     "Connecticut Low Income threshold: 30% or more people in households where the household income is less than or equal to twice the federal poverty level",
                                     "Massachusetts Low Income threshold: 25% or more of households with median household incomes below 65% statewide median",
                                     "Massachusetts Language Isolation threshold: 25% or more people in households where all adults speak English less than \"very well\"",
                                     "Massachusetts POC threshold: 40% or more are persons of Hispanic or Latino origin or persons who are not White OR 25% or more are persons of Hispanic or Latino origin or persons who are not White AND the median household income of the municipality is less than 150% of the statewide median household income",
                                     "Rhode Island Low Income threshold: highest 15% of block groups with people in households where the household income is less than or equal to twice the federal poverty level",
                                     "Rhode Island POC threshold: highest 15% of block groups with persons of Hispanic or Latino origin or persons who are not White"),
                     'ACS Table ID' = c("B03002: Hispanic or Latino Origin by Race",
                                "B03002: Hispanic or Latino Origin by Race",
                                "C17002: Ratio of Income to Poverty Level",
                                "C16002: Household Language by Household Limited English Speaking Status",
                                "B15002: Sex by Educational Attainment",
                                "B01001: Sex by Age",
                                "B01001: Sex by Age",
                                "B01001: Sex by Age",
                                "B18101: Sex by Age by Disability Status",
                                "B08201: Household Size by Vehicles Available",
                                "C17002: Ratio of Income to Poverty Level",
                                "B19001: Household Income",
                                "C16002: Household Language by Household Limited English Speaking Status",
                                "B03002: Hispanic or Latino Origin by Race",
                                "C17002: Ratio of Income to Poverty Level",
                                "B03002: Hispanic or Latino Origin by Race"),
                     Geography = c("Block Group",
                                   "Block Group",
                                   "Block Group",
                                   "Block Group",
                                   "Block Group",
                                   "Block Group",
                                   "Block Group",
                                   "Block Group",
                                   "Tract",
                                   "Tract",
                                   "Block Group",
                                   "Block Group",
                                   "Block Group",
                                   "Block Group",
                                   "Block Group",
                                   "Block Group"))

ACSvDF %>% 
  kableExtra::kable(longtable = T, booktabs = T, 
                    caption = "Demographic Variables", align = "l",
                    col.names = c("Variable","Description","ACS Table ID","Geography")) %>% 
  kableExtra::column_spec(2:3, width = "3cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
  
```


\pagebreak
# Appendix B: Supplementary Figures {-}

## Priority Populations in FEMA Flood Zones {-}
```{r tabNFHZANoCar, fig.align = "center", fig.cap="Households without a car living within FEMA Flood Zones."}
nh_tracts_nfhza %>% 
  st_centroid(.) %>% 
  st_intersection(nh_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  filter(TotalNoCar > 0) %>% 
  group_by(NAME) %>%
  summarize(`No Car HH` = sum(NewNoCar,na.rm = TRUE),
            `Pct of No Car HH` = round(sum(NewNoCar)/max(TotalNoCar)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  mutate(`No Car HH` = as.integer(`No Car HH`),
         `Pct of No Car HH` = paste0(`Pct of No Car HH`,"%")) %>% 
  filter(`No Car HH` > 0) %>% 
  arrange(desc(`No Car HH`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Households without a Car Living in FEMA Flood Zones", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r tabNFHZADIS, fig.align = "center", fig.cap="Disabled persons living within FEMA Flood Zones."}
nh_tracts_nfhza %>% 
  st_centroid(.) %>% 
  st_intersection(nh_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  filter(TotalDisabled > 0) %>% 
  group_by(NAME) %>%
  summarize(`Disabled` = sum(NewDisabled,na.rm = TRUE),
            `Pct of Disabled` = round(sum(NewDisabled)/max(TotalDisabled)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  mutate(`Disabled` = as.integer(`Disabled`),
         `Pct of Disabled` = paste0(`Pct of Disabled`,"%")) %>% 
  filter(`Disabled` > 0) %>% 
  arrange(desc(`Disabled`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Disabled Persons Living in FEMA Flood Zones", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapNFHZAPopConcern3, fig.align = "center", fig.cap="Low Income Persons Living within Flood Zones."}
# create pop subset
nfhza_sub <- nh_nfhza_vulnerable %>% 
  filter(Group == "Low Income")

# Map over-represented pops and flood zones
m <- tm_layout(bg.color = "#e6f3f7") +
  tm_shape(nh_blkgrp_2823, unit = "mi") + tm_fill(col = "white") +
  tm_shape(ne_states_sf_cb) + tm_fill(col="white") +
  tm_shape(canada) + tm_fill(col = "white") +
  # tm_shape(ny_state_sf) + tm_fill(col="white") +
  tm_shape(nh_awater_sf_simple) + tm_fill(col = "#e6f3f7") +
  tm_shape(nfhza_sub) +
  tm_dots(col = "darkgoldenrod3",
          labels = "1 dot = 10 people",
          alpha = 0.8) +
  # tm_facets(by = "Group",
  #           ncol = 2,
  #           free.coords = FALSE) +
  tm_shape(nh_nfhza_2823_land_simple) +
  tm_fill(col = "Interval",
          palette = c("dodgerblue", "blueviolet"),
          labels = c("1% AEP (100-year)", "0.2% AEP (500-year)"),
          title = "FEMA Flood Zones",
          alpha = 0.6,
          border.alpha = 0) +
  tm_shape(ne_states_sf_cb) + tm_borders(lwd = 0.6) +
  # tm_shape(ny_state_sf) + tm_borders(lwd = 0.6) +
  tm_shape(nh_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(nh_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_shape(nofemapolys) +
  tm_fill(col = "gray", alpha = 0.4) +
  tm_text("LABEL", alpha = 0.5, size = 0.6) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5,
               position = c("left","TOP")) +
  tm_add_legend(type = "fill", col = "darkgoldenrod3",
                alpha = 0.6,
                border.col = "white", border.alpha = 0,
                labels = "1 dot = 10 people") +
  tm_layout(title = "Low Income\nPersons\nLiving within\nFlood Zones",
            frame = TRUE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/FEMA/NH/NH_nhfza_vulnerable3.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/FEMA/NH/NH_nhfza_vulnerable3.png")
```

```{r tabNFHZALINC, fig.align = "center", fig.cap="Low income persons living within FEMA Flood Zones."}
nh_blkgrps_nfhza %>% 
  st_centroid(.) %>% 
  st_intersection(nh_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  filter(TotalLowInc > 0) %>% 
  group_by(NAME) %>%
  summarize(`Low Income` = sum(NewPov,na.rm = TRUE),
            `Pct of Low Income` = round(sum(NewPov)/max(TotalLowInc)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  mutate(`Low Income` = as.integer(`Low Income`),
         `Pct of Low Income` = paste0(`Pct of Low Income`,"%")) %>% 
  filter(`Low Income` > 0) %>% 
  arrange(desc(`Low Income`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Low Income Persons Living in FEMA Flood Zones", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapNFHZAPopConcern4, fig.align = "center", fig.cap="Adults without a High School Diploma Living within Flood Zones."}
# create pop subset
nfhza_sub <- nh_nfhza_vulnerable %>% 
  filter(Group == "No HS Dip")

# Map over-represented pops and flood zones
m <- tm_layout(bg.color = "#e6f3f7") +
  tm_shape(nh_blkgrp_2823, unit = "mi") + tm_fill(col = "white") +
  tm_shape(ne_states_sf_cb) + tm_fill(col="white") +
  tm_shape(canada) + tm_fill(col = "white") +
  # tm_shape(ny_state_sf) + tm_fill(col="white") +
  tm_shape(nh_awater_sf_simple) + tm_fill(col = "#e6f3f7") +
  tm_shape(nfhza_sub) +
  tm_dots(col = "darkgoldenrod3",
          labels = "1 dot = 10 people",
          alpha = 0.8) +
  # tm_facets(by = "Group",
  #           ncol = 2,
  #           free.coords = FALSE) +
  tm_shape(nh_nfhza_2823_land_simple) +
  tm_fill(col = "Interval",
          palette = c("dodgerblue", "blueviolet"),
          labels = c("1% AEP (100-year)", "0.2% AEP (500-year)"),
          title = "FEMA Flood Zones",
          alpha = 0.6,
          border.alpha = 0) +
  tm_shape(ne_states_sf_cb) + tm_borders(lwd = 0.6) +
  # tm_shape(ny_state_sf) + tm_borders(lwd = 0.6) +
  tm_shape(nh_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(nh_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_shape(nofemapolys) +
  tm_fill(col = "gray", alpha = 0.4) +
  tm_text("LABEL", alpha = 0.5, size = 0.6) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5,
               position = c("left","TOP")) +
  tm_add_legend(type = "fill", col = "darkgoldenrod3",
                alpha = 0.6,
                border.col = "white", border.alpha = 0,
                labels = "1 dot = 10 people") +
  tm_layout(title = "Adults without a\nHigh School\nDiploma\nLiving within\nFlood Zones",
            frame = TRUE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/FEMA/NH/NH_nhfza_vulnerable4.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/FEMA/NH/NH_nhfza_vulnerable4.png")
```

```{r tabNFHZALTHS, fig.align = "center", fig.cap="Adults without a high school diploma living within FEMA Flood Zones."}
nh_blkgrps_nfhza %>% 
  st_centroid(.) %>% 
  st_intersection(nh_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  filter(TotalNoHS > 0) %>% 
  group_by(NAME) %>%
  summarize(`No HS Dip` = sum(NewLths,na.rm = TRUE),
            `Pct of No HS Dip` = round(sum(NewLths)/max(TotalNoHS)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  mutate(`No HS Dip` = as.integer(`No HS Dip`),
         `Pct of No HS Dip` = paste0(`Pct of No HS Dip`,"%")) %>% 
  filter(`No HS Dip` > 0) %>% 
  arrange(desc(`No HS Dip`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Adults without a High School Diploma Living in FEMA Flood Zones", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapNFHZAPopConcern5, fig.align = "center", fig.cap="Persons Over 64 Living within Flood Zones."}
# create pop subset
nfhza_sub <- nh_nfhza_vulnerable %>% 
  filter(Group == "Over 64")

# Map over-represented pops and flood zones
m <- tm_layout(bg.color = "#e6f3f7") +
  tm_shape(nh_blkgrp_2823, unit = "mi") + tm_fill(col = "white") +
  tm_shape(ne_states_sf_cb) + tm_fill(col="white") +
  tm_shape(canada) + tm_fill(col = "white") +
  # tm_shape(ny_state_sf) + tm_fill(col="white") +
  tm_shape(nh_awater_sf_simple) + tm_fill(col = "#e6f3f7") +
  tm_shape(nfhza_sub) +
  tm_dots(col = "darkgoldenrod3",
          labels = "1 dot = 10 people",
          alpha = 0.8) +
  # tm_facets(by = "Group",
  #           ncol = 2,
  #           free.coords = FALSE) +
  tm_shape(nh_nfhza_2823_land_simple) +
  tm_fill(col = "Interval",
          palette = c("dodgerblue", "blueviolet"),
          labels = c("1% AEP (100-year)", "0.2% AEP (500-year)"),
          title = "FEMA Flood Zones",
          alpha = 0.6,
          border.alpha = 0) +
  tm_shape(ne_states_sf_cb) + tm_borders(lwd = 0.6) +
  # tm_shape(ny_state_sf) + tm_borders(lwd = 0.6) +
  tm_shape(nh_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(nh_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_shape(nofemapolys) +
  tm_fill(col = "gray", alpha = 0.4) +
  tm_text("LABEL", alpha = 0.5, size = 0.6) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5,
               position = c("left","TOP")) +
  tm_add_legend(type = "fill", col = "darkgoldenrod3",
                alpha = 0.6,
                border.col = "white", border.alpha = 0,
                labels = "1 dot = 10 people") +
  tm_layout(title = "Persons\nOver 64\nLiving within\nFlood Zones",
            frame = TRUE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/FEMA/NH/NH_nhfza_vulnerable5.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/FEMA/NH/NH_nhfza_vulnerable5.png")
```

```{r tabNFHZAO64, fig.align = "center", fig.cap="Persons over 64 living within FEMA Flood Zones."}
nh_blkgrps_nfhza %>% 
  st_centroid(.) %>% 
  st_intersection(nh_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  filter(TotalOver64 > 0) %>% 
  group_by(NAME) %>%
  summarize(`Over 64` = sum(NewOver64,na.rm = TRUE),
            `Pct of Over 64` = round(sum(NewOver64)/max(TotalOver64)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  mutate(`Over 64` = as.integer(`Over 64`),
         `Pct of Over 64` = paste0(`Pct of Over 64`,"%")) %>% 
  filter(`Over 64` > 0) %>% 
  arrange(desc(`Over 64`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Persons Over 64 Living in FEMA Flood Zones", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapNFHZAPopConcern6, fig.align = "center", fig.cap="People of Color Living within Flood Zones."}
# create pop subset
nfhza_sub <- nh_nfhza_vulnerable %>% 
  filter(Group == "Minority")

# Map over-represented pops and flood zones
m <- tm_layout(bg.color = "#e6f3f7") +
  tm_shape(nh_blkgrp_2823, unit = "mi") + tm_fill(col = "white") +
  tm_shape(ne_states_sf_cb) + tm_fill(col="white") +
  tm_shape(canada) + tm_fill(col = "white") +
  # tm_shape(ny_state_sf) + tm_fill(col="white") +
  tm_shape(nh_awater_sf_simple) + tm_fill(col = "#e6f3f7") +
  tm_shape(nfhza_sub) +
  tm_dots(col = "darkgoldenrod3",
          labels = "1 dot = 10 people",
          alpha = 0.8) +
  # tm_facets(by = "Group",
  #           ncol = 2,
  #           free.coords = FALSE) +
  tm_shape(nh_nfhza_2823_land_simple) +
  tm_fill(col = "Interval",
          palette = c("dodgerblue", "blueviolet"),
          labels = c("1% AEP (100-year)", "0.2% AEP (500-year)"),
          title = "FEMA Flood Zones",
          alpha = 0.6,
          border.alpha = 0) +
  tm_shape(ne_states_sf_cb) + tm_borders(lwd = 0.6) +
  # tm_shape(ny_state_sf) + tm_borders(lwd = 0.6) +
  tm_shape(nh_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(nh_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_shape(nofemapolys) +
  tm_fill(col = "gray", alpha = 0.4) +
  tm_text("LABEL", alpha = 0.5, size = 0.6) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5,
               position = c("left","TOP")) +
  tm_add_legend(type = "fill", col = "darkgoldenrod3",
                alpha = 0.6,
                border.col = "white", border.alpha = 0,
                labels = "1 dot = 10 people") +
  tm_layout(title = "People\nof Color\nLiving within\nFlood Zones",
            frame = TRUE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/FEMA/NH/NH_nhfza_vulnerable6.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/FEMA/NH/NH_nhfza_vulnerable6.png")
```

```{r tabNFHZAMIN, fig.align = "center", fig.cap="People of Color living within FEMA Flood Zones."}
nh_blkgrps_nfhza %>% 
  st_centroid(.) %>% 
  st_intersection(nh_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  filter(TotalMin > 0) %>% 
  group_by(NAME) %>%
  summarize(`Minority` = sum(NewMinority,na.rm = TRUE),
            `Pct of Minority` = round(sum(NewMinority)/max(TotalMin)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  mutate(`Minority` = as.integer(`Minority`),
         `Pct of Minority` = paste0(`Pct of Minority`,"%")) %>% 
  filter(`Minority` > 0) %>% 
  arrange(desc(`Minority`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "People of Color Living in FEMA Flood Zones", align = "r", col.names = c("City/Town", "People of Color", "Pct of POC")) %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapNFHZAPopConcern7, fig.align = "center", fig.cap="Limited English Speaking Households Living within Flood Zones."}
# create pop subset
nfhza_sub <- nh_nfhza_vulnerable %>% 
  filter(Group == "Limited English HH")

# Map over-represented pops and flood zones
m <- tm_layout(bg.color = "#e6f3f7") +
  tm_shape(nh_blkgrp_2823, unit = "mi") + tm_fill(col = "white") +
  tm_shape(ne_states_sf_cb) + tm_fill(col="white") +
  tm_shape(canada) + tm_fill(col = "white") +
  # tm_shape(ny_state_sf) + tm_fill(col="white") +
  tm_shape(nh_awater_sf_simple) + tm_fill(col = "#e6f3f7") +
  tm_shape(nfhza_sub) +
  tm_dots(col = "darkgoldenrod3",
          labels = "1 dot = 10 people",
          alpha = 0.8) +
  # tm_facets(by = "Group",
  #           ncol = 2,
  #           free.coords = FALSE) +
  tm_shape(nh_nfhza_2823_land_simple) +
  tm_fill(col = "Interval",
          palette = c("dodgerblue", "blueviolet"),
          labels = c("1% AEP (100-year)", "0.2% AEP (500-year)"),
          title = "FEMA Flood Zones",
          alpha = 0.6,
          border.alpha = 0) +
  tm_shape(ne_states_sf_cb) + tm_borders(lwd = 0.6) +
  # tm_shape(ny_state_sf) + tm_borders(lwd = 0.6) +
  tm_shape(nh_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(nh_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_shape(nofemapolys) +
  tm_fill(col = "gray", alpha = 0.4) +
  tm_text("LABEL", alpha = 0.5, size = 0.6) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5,
               position = c("left","TOP")) +
  tm_add_legend(type = "fill", col = "darkgoldenrod3",
                alpha = 0.6,
                border.col = "white", border.alpha = 0,
                labels = "1 dot = 10 people") +
  tm_layout(title = "Limited\nEnglish\nSpeaking\nHouseholds\nLiving within\nFlood Zones",
            frame = TRUE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/FEMA/NH/NH_nhfza_vulnerable7.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/FEMA/NH/NH_nhfza_vulnerable7.png")
```

```{r tabNFHZALEH, fig.align = "center", fig.cap="Limited English speaking households living within FEMA Flood Zones."}
nh_blkgrps_nfhza %>% 
  st_centroid(.) %>% 
  st_intersection(nh_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  filter(TotalLEH > 0) %>% 
  group_by(NAME) %>%
  summarize(`Limited English HH` = sum(NewEng_limit,na.rm = TRUE),
            `Pct of Limited English HH` = round(sum(NewEng_limit)/max(TotalLEH)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  mutate(`Limited English HH` = as.integer(`Limited English HH`),
         `Pct of Limited English HH` = paste0(`Pct of Limited English HH`,"%")) %>% 
  filter(`Limited English HH` > 0) %>% 
  arrange(desc(`Limited English HH`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Limited English Speaking Households Living in FEMA Flood Zones", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapNFHZAPopConcern8, fig.align = "center", fig.cap="Persons Under 18 Living within Flood Zones."}
# create pop subset
nfhza_sub <- nh_nfhza_vulnerable %>% 
  filter(Group == "Under 18")

# Map over-represented pops and flood zones
m <- tm_layout(bg.color = "#e6f3f7") +
  tm_shape(nh_blkgrp_2823, unit = "mi") + tm_fill(col = "white") +
  tm_shape(ne_states_sf_cb) + tm_fill(col="white") +
  tm_shape(canada) + tm_fill(col = "white") +
  # tm_shape(ny_state_sf) + tm_fill(col="white") +
  tm_shape(nh_awater_sf_simple) + tm_fill(col = "#e6f3f7") +
  tm_shape(nfhza_sub) +
  tm_dots(col = "darkgoldenrod3",
          labels = "1 dot = 10 people",
          alpha = 0.8) +
  # tm_facets(by = "Group",
  #           ncol = 2,
  #           free.coords = FALSE) +
  tm_shape(nh_nfhza_2823_land_simple) +
  tm_fill(col = "Interval",
          palette = c("dodgerblue", "blueviolet"),
          labels = c("1% AEP (100-year)", "0.2% AEP (500-year)"),
          title = "FEMA Flood Zones",
          alpha = 0.6,
          border.alpha = 0) +
  tm_shape(ne_states_sf_cb) + tm_borders(lwd = 0.6) +
  # tm_shape(ny_state_sf) + tm_borders(lwd = 0.6) +
  tm_shape(nh_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(nh_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_shape(nofemapolys) +
  tm_fill(col = "gray", alpha = 0.4) +
  tm_text("LABEL", alpha = 0.5, size = 0.6) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5,
               position = c("left","TOP")) +
  tm_add_legend(type = "fill", col = "darkgoldenrod3",
                alpha = 0.6,
                border.col = "white", border.alpha = 0,
                labels = "1 dot = 10 people") +
  tm_layout(title = "Persons\nUnder 18\nLiving within\nFlood Zones",
            frame = TRUE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/FEMA/NH/NH_nhfza_vulnerable8.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/FEMA/NH/NH_nhfza_vulnerable8.png")
```

```{r tabNFHZAU18, fig.align = "center", fig.cap="Persons under 18 living within FEMA Flood Zones."}
nh_blkgrps_nfhza %>% 
  st_centroid(.) %>% 
  st_intersection(nh_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  filter(TotalUnder18 > 0) %>% 
  group_by(NAME) %>%
  summarize(`Under 18` = sum(NewUnder18,na.rm = TRUE),
            `Pct of Under 18` = round(sum(NewUnder18)/max(TotalUnder18)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  mutate(`Under 18` = as.integer(`Under 18`),
         `Pct of Under 18` = paste0(`Pct of Under 18`,"%")) %>% 
  filter(`Under 18` > 0) %>% 
  arrange(desc(`Under 18`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Persons Under 18 Living in FEMA Flood Zones", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapNFHZAPopConcern9, fig.align = "center", fig.cap="Children Under 5 Living within Flood Zones."}
# create pop subset
nfhza_sub <- nh_nfhza_vulnerable %>% 
  filter(Group == "Under 5")

# Map over-represented pops and flood zones
m <- tm_layout(bg.color = "#e6f3f7") +
  tm_shape(nh_blkgrp_2823, unit = "mi") + tm_fill(col = "white") +
  tm_shape(ne_states_sf_cb) + tm_fill(col="white") +
  tm_shape(canada) + tm_fill(col = "white") +
  # tm_shape(ny_state_sf) + tm_fill(col="white") +
  tm_shape(nh_awater_sf_simple) + tm_fill(col = "#e6f3f7") +
  tm_shape(nfhza_sub) +
  tm_dots(col = "darkgoldenrod3",
          labels = "1 dot = 10 people",
          alpha = 0.8) +
  # tm_facets(by = "Group",
  #           ncol = 2,
  #           free.coords = FALSE) +
  tm_shape(nh_nfhza_2823_land_simple) +
  tm_fill(col = "Interval",
          palette = c("dodgerblue", "blueviolet"),
          labels = c("1% AEP (100-year)", "0.2% AEP (500-year)"),
          title = "FEMA Flood Zones",
          alpha = 0.6,
          border.alpha = 0) +
  tm_shape(ne_states_sf_cb) + tm_borders(lwd = 0.6) +
  # tm_shape(ny_state_sf) + tm_borders(lwd = 0.6) +
  tm_shape(nh_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(nh_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_shape(nofemapolys) +
  tm_fill(col = "gray", alpha = 0.4) +
  tm_text("LABEL", alpha = 0.5, size = 0.6) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5,
               position = c("left","TOP")) +
  tm_add_legend(type = "fill", col = "darkgoldenrod3",
                alpha = 0.6,
                border.col = "white", border.alpha = 0,
                labels = "1 dot = 10 people") +
  tm_layout(title = "Children\nUnder 5\nLiving within\nFlood Zones",
            frame = TRUE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/FEMA/NH/NH_nhfza_vulnerable9.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/FEMA/NH/NH_nhfza_vulnerable9.png")
```

```{r tabNFHZAU5, fig.align = "center", fig.cap="Children under 5 living within FEMA Flood Zones."}
nh_blkgrps_nfhza %>% 
  st_centroid(.) %>% 
  st_intersection(nh_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  filter(TotalUnder5 > 0) %>% 
  group_by(NAME) %>%
  summarize(`Under 5` = sum(NewUnder5,na.rm = TRUE),
            `Pct of Under 5` = round(sum(NewUnder5)/max(TotalUnder5)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  mutate(`Under 5` = as.integer(`Under 5`),
         `Pct of Under 5` = paste0(`Pct of Under 5`,"%")) %>% 
  filter(`Under 5` > 0) %>% 
  arrange(desc(`Under 5`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Children Under 5 Living in FEMA Flood Zones", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```


\pagebreak
## Priority Populations in Hurricane Evacuation Zones {-}

```{r tabHEVACO64, fig.align = "center", fig.cap="Persons over 64 living within Hurricane Inundation Zones."}
nh_blkgrps_hevac %>% 
  st_centroid(.) %>% 
  st_intersection(nh_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  filter(TotalOver64 > 0) %>% 
  group_by(NAME) %>%
  summarize(`Over 64` = sum(NewOver64,na.rm = TRUE),
            `Pct of Over 64` = round(sum(NewOver64)/max(TotalOver64)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  mutate(`Over 64` = as.integer(`Over 64`),
         `Pct of Over 64` = paste0(`Pct of Over 64`,"%")) %>% 
  filter(`Over 64` > 0) %>% 
  arrange(desc(`Over 64`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Persons Over 64 Living in Hurricane Inundation Zones", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r tabHEVACNoCar, fig.align = "center", fig.cap="Households without a car living within Hurricane Inundation Zones."}
nh_tracts_hevac %>% 
  st_centroid(.) %>% 
  st_intersection(nh_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  filter(TotalNoCar > 0) %>% 
  group_by(NAME) %>%
  summarize(`No Car HH` = sum(NewNoCar,na.rm = TRUE),
            `Pct of No Car HH` = round(sum(NewNoCar)/max(TotalNoCar)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  mutate(`No Car HH` = as.integer(`No Car HH`),
         `Pct of No Car HH` = paste0(`Pct of No Car HH`,"%")) %>% 
  filter(`No Car HH` > 0) %>% 
  arrange(desc(`No Car HH`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Households without a Car Living in Hurricane Inundation Zones", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapHEVACPopConcern3, fig.align = "center", fig.cap="Disabled Persons Living within Hurricane Inundation Zones."}
# Create data subset
hevac_sub <- nh_hevac_vulnerable %>% 
  filter(Group == "Disabled")

# Map over-represented pops and hurricane evacuation zones
m2 <- tm_layout(bg.color = "#e6f3f7") +
  tm_shape(nh_hea_sf_simple, unit = "mi") + tm_fill(col = "white") +
  tm_shape(nh_blkgrp_2823, unit = "mi") + tm_fill(col = "white") +
  tm_shape(ne_states_sf_cb) + tm_fill(col="white") +
  # tm_shape(ny_state_sf) + tm_fill(col="white") +
  tm_shape(coastalWater) + tm_fill(col = "#e6f3f7") + 
  tm_shape(ne_states_sf_cb) + tm_borders() +
  # tm_shape(ny_state_sf) + tm_borders() +
  tm_shape(nh_hea_sf) + 
  tm_fill(col = "Zone", 
          palette = c("deeppink", "purple3"), 
          labels = c("A: Category 1 - 2", 
                     "B: Category 3 - 4"), 
          title = "Inundation Zone and\nHurricane Category",
          alpha = 0.6,
          border.alpha = 0) +
  tm_shape(nh_hevac_vulnerable) +
  tm_dots(col = "springgreen3",
          labels = "1 dot = 10 people",
          alpha = 1) +
  # tm_facets(by = "Group",
  #           ncol = 2,
  #           free.coords = FALSE) +
  tm_shape(ne_states_sf_cb) + tm_borders(lwd = 0.2) +
  # tm_shape(ny_state_sf) + tm_borders(lwd = 0.8) +
  tm_shape(nh_highways) + tm_lines(col = "seashell4", lwd = 1) +
  # tm_shape(nh_highways2nd) + tm_lines(col = "seashell4", lwd = 1.5) +
  # tm_shape(I89roadSegment) + 
  # tm_symbols(shape = I89, border.lwd = NA, size = 0.2) +
  # tm_shape(I91roadSegment) + 
  # tm_symbols(shape = I91, border.lwd = NA, size = 0.2) +
  # tm_shape(I93roadSegment) +
  # tm_symbols(shape = I93, border.lwd = NA, size = 0.2) +
  tm_shape(I95roadSegment2) + 
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(nh_towns_sf_pts2) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c("left","TOP")) +
  tm_add_legend(type = "fill", col = "springgreen3",
                alpha = 0.6,
                border.col = "white", border.alpha = 0,
                labels = "1 dot = 10 people") +
  tm_layout(title = "Disabled\nPersons\nLiving within\nHurricane\nInundation Zones", frame = TRUE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m2, "DATA/FEMA/NH/NH_hevac_vulnerable3.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/FEMA/NH/NH_hevac_vulnerable3.png")
```

```{r tabHEVACDIS, fig.align = "center", fig.cap="Disabled persons living within Hurricane Inundation Zones."}
nh_tracts_hevac %>% 
  st_centroid(.) %>% 
  st_intersection(nh_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  filter(TotalDisabled > 0) %>% 
  group_by(NAME) %>%
  summarize(`Disabled` = sum(NewDisabled,na.rm = TRUE),
            `Pct of Disabled` = round(sum(NewDisabled)/max(TotalDisabled)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  mutate(`Disabled` = as.integer(`Disabled`),
         `Pct of Disabled` = paste0(`Pct of Disabled`,"%")) %>% 
  filter(`Disabled` > 0) %>% 
  arrange(desc(`Disabled`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Disabled Persons Living in Hurricane Inundation Zones", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapHEVACPopConcern4, fig.align = "center", fig.cap="Low Income Persons Living within Hurricane Inundation Zones."}
# Create data subset
hevac_sub <- nh_hevac_vulnerable %>% 
  filter(Group == "Low Income")

# Map over-represented pops and hurricane evacuation zones
m2 <- tm_layout(bg.color = "#e6f3f7") +
  tm_shape(nh_hea_sf_simple, unit = "mi") + tm_fill(col = "white") +
  tm_shape(nh_blkgrp_2823, unit = "mi") + tm_fill(col = "white") +
  tm_shape(ne_states_sf_cb) + tm_fill(col="white") +
  # tm_shape(ny_state_sf) + tm_fill(col="white") +
  tm_shape(coastalWater) + tm_fill(col = "#e6f3f7") + 
  tm_shape(ne_states_sf_cb) + tm_borders() +
  # tm_shape(ny_state_sf) + tm_borders() +
  tm_shape(nh_hea_sf) + 
  tm_fill(col = "Zone", 
          palette = c("deeppink", "purple3"), 
          labels = c("A: Category 1 - 2", 
                     "B: Category 3 - 4"), 
          title = "Inundation Zone and\nHurricane Category",
          alpha = 0.6,
          border.alpha = 0) +
  tm_shape(nh_hevac_vulnerable) +
  tm_dots(col = "springgreen3",
          labels = "1 dot = 10 people",
          alpha = 1) +
  # tm_facets(by = "Group",
  #           ncol = 2,
  #           free.coords = FALSE) +
  tm_shape(ne_states_sf_cb) + tm_borders(lwd = 0.2) +
  # tm_shape(ny_state_sf) + tm_borders(lwd = 0.8) +
  tm_shape(nh_highways) + tm_lines(col = "seashell4", lwd = 1) +
  # tm_shape(nh_highways2nd) + tm_lines(col = "seashell4", lwd = 1.5) +
  # tm_shape(I89roadSegment) + 
  # tm_symbols(shape = I89, border.lwd = NA, size = 0.2) +
  # tm_shape(I91roadSegment) + 
  # tm_symbols(shape = I91, border.lwd = NA, size = 0.2) +
  # tm_shape(I93roadSegment) +
  # tm_symbols(shape = I93, border.lwd = NA, size = 0.2) +
  tm_shape(I95roadSegment2) + 
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(nh_towns_sf_pts2) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c("left","TOP")) +
  tm_add_legend(type = "fill", col = "springgreen3",
                alpha = 0.6,
                border.col = "white", border.alpha = 0,
                labels = "1 dot = 10 people") +
  tm_layout(title = "Low Income\nPersons\nLiving within\nHurricane\nInundation Zones", frame = TRUE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m2, "DATA/FEMA/NH/NH_hevac_vulnerable4.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/FEMA/NH/NH_hevac_vulnerable4.png")
```

```{r tabHEVACLINC, fig.align = "center", fig.cap="Low income persons living within Hurricane Inundation Zones."}
nh_blkgrps_hevac %>% 
  st_centroid(.) %>% 
  st_intersection(nh_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  filter(TotalLowInc > 0) %>% 
  group_by(NAME) %>%
  summarize(`Low Income` = sum(NewPov,na.rm = TRUE),
            `Pct of Low Income` = round(sum(NewPov)/max(TotalLowInc)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  mutate(`Low Income` = as.integer(`Low Income`),
         `Pct of Low Income` = paste0(`Pct of Low Income`,"%")) %>% 
  filter(`Low Income` > 0) %>% 
  arrange(desc(`Low Income`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Low Income Persons Living in Hurricane Inundation Zones", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapHEVACPopConcern5, fig.align = "center", fig.cap="Persons Under 18 Living within Hurricane Inundation Zones."}
# Create data subset
hevac_sub <- nh_hevac_vulnerable %>% 
  filter(Group == "Under 18")

# Map over-represented pops and hurricane evacuation zones
m2 <- tm_layout(bg.color = "#e6f3f7") +
  tm_shape(nh_hea_sf_simple, unit = "mi") + tm_fill(col = "white") +
  tm_shape(nh_blkgrp_2823, unit = "mi") + tm_fill(col = "white") +
  tm_shape(ne_states_sf_cb) + tm_fill(col="white") +
  # tm_shape(ny_state_sf) + tm_fill(col="white") +
  tm_shape(coastalWater) + tm_fill(col = "#e6f3f7") + 
  tm_shape(ne_states_sf_cb) + tm_borders() +
  # tm_shape(ny_state_sf) + tm_borders() +
  tm_shape(nh_hea_sf) + 
  tm_fill(col = "Zone", 
          palette = c("deeppink", "purple3"), 
          labels = c("A: Category 1 - 2", 
                     "B: Category 3 - 4"), 
          title = "Inundation Zone and\nHurricane Category",
          alpha = 0.6,
          border.alpha = 0) +
  tm_shape(nh_hevac_vulnerable) +
  tm_dots(col = "springgreen3",
          labels = "1 dot = 10 people",
          alpha = 1) +
  # tm_facets(by = "Group",
  #           ncol = 2,
  #           free.coords = FALSE) +
  tm_shape(ne_states_sf_cb) + tm_borders(lwd = 0.2) +
  # tm_shape(ny_state_sf) + tm_borders(lwd = 0.8) +
  tm_shape(nh_highways) + tm_lines(col = "seashell4", lwd = 1) +
  # tm_shape(nh_highways2nd) + tm_lines(col = "seashell4", lwd = 1.5) +
  # tm_shape(I89roadSegment) + 
  # tm_symbols(shape = I89, border.lwd = NA, size = 0.2) +
  # tm_shape(I91roadSegment) + 
  # tm_symbols(shape = I91, border.lwd = NA, size = 0.2) +
  # tm_shape(I93roadSegment) +
  # tm_symbols(shape = I93, border.lwd = NA, size = 0.2) +
  tm_shape(I95roadSegment2) + 
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(nh_towns_sf_pts2) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c("left","TOP")) +
  tm_add_legend(type = "fill", col = "springgreen3",
                alpha = 0.6,
                border.col = "white", border.alpha = 0,
                labels = "1 dot = 10 people") +
  tm_layout(title = "Persons\nUnder 18\nLiving within\nHurricane\nInundation Zones", frame = TRUE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m2, "DATA/FEMA/NH/NH_hevac_vulnerable5.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/FEMA/NH/NH_hevac_vulnerable5.png")
```

```{r tabHEVACU18, fig.align = "center", fig.cap="Persons under 18 living within Hurricane Inundation Zones."}
nh_blkgrps_hevac %>% 
  st_centroid(.) %>% 
  st_intersection(nh_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  filter(TotalUnder18 > 0) %>% 
  group_by(NAME) %>%
  summarize(`Under 18` = sum(NewUnder18,na.rm = TRUE),
            `Pct of Under 18` = round(sum(NewUnder18)/max(TotalUnder18)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  mutate(`Under 18` = as.integer(`Under 18`),
         `Pct of Under 18` = paste0(`Pct of Under 18`,"%")) %>% 
  filter(`Under 18` > 0) %>% 
  arrange(desc(`Under 18`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Persons Under 18 Living in Hurricane Inundation Zones", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapHEVACPopConcern6, fig.align = "center", fig.cap="Adults without a High School Living within Hurricane Inundation Zones."}
# Create data subset
hevac_sub <- nh_hevac_vulnerable %>% 
  filter(Group == "No HS Dip")

# Map over-represented pops and hurricane evacuation zones
m2 <- tm_layout(bg.color = "#e6f3f7") +
  tm_shape(nh_hea_sf_simple, unit = "mi") + tm_fill(col = "white") +
  tm_shape(nh_blkgrp_2823, unit = "mi") + tm_fill(col = "white") +
  tm_shape(ne_states_sf_cb) + tm_fill(col="white") +
  # tm_shape(ny_state_sf) + tm_fill(col="white") +
  tm_shape(coastalWater) + tm_fill(col = "#e6f3f7") + 
  tm_shape(ne_states_sf_cb) + tm_borders() +
  # tm_shape(ny_state_sf) + tm_borders() +
  tm_shape(nh_hea_sf) + 
  tm_fill(col = "Zone", 
          palette = c("deeppink", "purple3"), 
          labels = c("A: Category 1 - 2", 
                     "B: Category 3 - 4"), 
          title = "Inundation Zone and\nHurricane Category",
          alpha = 0.6,
          border.alpha = 0) +
  tm_shape(nh_hevac_vulnerable) +
  tm_dots(col = "springgreen3",
          labels = "1 dot = 10 people",
          alpha = 1) +
  # tm_facets(by = "Group",
  #           ncol = 2,
  #           free.coords = FALSE) +
  tm_shape(ne_states_sf_cb) + tm_borders(lwd = 0.2) +
  # tm_shape(ny_state_sf) + tm_borders(lwd = 0.8) +
  tm_shape(nh_highways) + tm_lines(col = "seashell4", lwd = 1) +
  # tm_shape(nh_highways2nd) + tm_lines(col = "seashell4", lwd = 1.5) +
  # tm_shape(I89roadSegment) + 
  # tm_symbols(shape = I89, border.lwd = NA, size = 0.2) +
  # tm_shape(I91roadSegment) + 
  # tm_symbols(shape = I91, border.lwd = NA, size = 0.2) +
  # tm_shape(I93roadSegment) +
  # tm_symbols(shape = I93, border.lwd = NA, size = 0.2) +
  tm_shape(I95roadSegment2) + 
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(nh_towns_sf_pts2) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c("left","TOP")) +
  tm_add_legend(type = "fill", col = "springgreen3",
                alpha = 0.6,
                border.col = "white", border.alpha = 0,
                labels = "1 dot = 10 people") +
  tm_layout(title = "Adults without a\nHigh School\nDiploma\nLiving within\nHurricane\nInundation Zones", frame = TRUE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m2, "DATA/FEMA/NH/NH_hevac_vulnerable6.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/FEMA/NH/NH_hevac_vulnerable6.png")
```

```{r tabHEVACLTHS, fig.align = "center", fig.cap="Adults without a high school diploma living within Hurricane Inundation Zones."}
nh_blkgrps_hevac %>% 
  st_centroid(.) %>% 
  st_intersection(nh_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  filter(TotalNoHS > 0) %>% 
  group_by(NAME) %>%
  summarize(`No HS Dip` = sum(NewLths,na.rm = TRUE),
            `Pct of No HS Dip` = round(sum(NewLths)/max(TotalNoHS)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  mutate(`No HS Dip` = as.integer(`No HS Dip`),
         `Pct of No HS Dip` = paste0(`Pct of No HS Dip`,"%")) %>% 
  filter(`No HS Dip` > 0) %>% 
  arrange(desc(`No HS Dip`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Adults wihtout a High School Diploma Living in Hurricane Inundation Zones", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapHEVACPopConcern7, fig.align = "center", fig.cap="People of Color Living within Hurricane Inundation Zones."}
# Create data subset
hevac_sub <- nh_hevac_vulnerable %>% 
  filter(Group == "Minority")

# Map over-represented pops and hurricane evacuation zones
m2 <- tm_layout(bg.color = "#e6f3f7") +
  tm_shape(nh_hea_sf_simple, unit = "mi") + tm_fill(col = "white") +
  tm_shape(nh_blkgrp_2823, unit = "mi") + tm_fill(col = "white") +
  tm_shape(ne_states_sf_cb) + tm_fill(col="white") +
  # tm_shape(ny_state_sf) + tm_fill(col="white") +
  tm_shape(coastalWater) + tm_fill(col = "#e6f3f7") + 
  tm_shape(ne_states_sf_cb) + tm_borders() +
  # tm_shape(ny_state_sf) + tm_borders() +
  tm_shape(nh_hea_sf) + 
  tm_fill(col = "Zone", 
          palette = c("deeppink", "purple3"), 
          labels = c("A: Category 1 - 2", 
                     "B: Category 3 - 4"), 
          title = "Inundation Zone and\nHurricane Category",
          alpha = 0.6,
          border.alpha = 0) +
  tm_shape(nh_hevac_vulnerable) +
  tm_dots(col = "springgreen3",
          labels = "1 dot = 10 people",
          alpha = 1) +
  # tm_facets(by = "Group",
  #           ncol = 2,
  #           free.coords = FALSE) +
  tm_shape(ne_states_sf_cb) + tm_borders(lwd = 0.2) +
  # tm_shape(ny_state_sf) + tm_borders(lwd = 0.8) +
  tm_shape(nh_highways) + tm_lines(col = "seashell4", lwd = 1) +
  # tm_shape(nh_highways2nd) + tm_lines(col = "seashell4", lwd = 1.5) +
  # tm_shape(I89roadSegment) + 
  # tm_symbols(shape = I89, border.lwd = NA, size = 0.2) +
  # tm_shape(I91roadSegment) + 
  # tm_symbols(shape = I91, border.lwd = NA, size = 0.2) +
  # tm_shape(I93roadSegment) +
  # tm_symbols(shape = I93, border.lwd = NA, size = 0.2) +
  tm_shape(I95roadSegment2) + 
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(nh_towns_sf_pts2) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c("left","TOP")) +
  tm_add_legend(type = "fill", col = "springgreen3",
                alpha = 0.6,
                border.col = "white", border.alpha = 0,
                labels = "1 dot = 10 people") +
  tm_layout(title = "People\nof Color\nLiving within\nHurricane\nInundation Zones", frame = TRUE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m2, "DATA/FEMA/NH/NH_hevac_vulnerable7.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/FEMA/NH/NH_hevac_vulnerable7.png")
```

```{r tabHEVACMIN, fig.align = "center", fig.cap="People of Color living within Hurricane Inundation Zones."}
nh_blkgrps_hevac %>% 
  st_centroid(.) %>% 
  st_intersection(nh_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  filter(TotalMin > 0) %>% 
  group_by(NAME) %>%
  summarize(`Minority` = sum(NewMinority,na.rm = TRUE),
            `Pct of Minority` = round(sum(NewMinority)/max(TotalMin)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  mutate(`Minority` = as.integer(`Minority`),
         `Pct of Minority` = paste0(`Pct of Minority`,"%")) %>% 
  filter(`Minority` > 0) %>% 
  arrange(desc(`Minority`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "People of Color Living in Hurricane Inundation Zones", align = "r", col.names = c("City/Town", "People of Color", "Pct of POC")) %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header","striped")) 
```

```{r tabHEVACLEH, fig.align = "center", fig.cap="Limited English speaking households living within Hurricane Inundation Zones."}
nh_blkgrps_hevac %>% 
  st_centroid(.) %>% 
  st_intersection(nh_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  filter(TotalLEH > 0) %>% 
  group_by(NAME) %>%
  summarize(`Limited English HH` = sum(NewEng_limit,na.rm = TRUE),
            `Pct of Limited English HH` = round(sum(NewEng_limit)/max(TotalLEH)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  mutate(`Limited English HH` = as.integer(`Limited English HH`),
         `Pct of Limited English HH` = paste0(`Pct of Limited English HH`,"%")) %>% 
  filter(`Limited English HH` > 0) %>% 
  arrange(desc(`Limited English HH`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Limited English Speaking Households Living in Hurricane Inundation Zones", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header","striped")) 
```