Evacuation_VermontPDF.Rmd

---
title: "Evacuation Risks in Vermont"
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(rmapshaper)
library(parallel)
library(tigris)
options(tigris_use_cache = TRUE, tigris_class = "sf")
library(kableExtra)
```
\pagebreak
# Analysis of evacuation-related risks in Vermont

This is an analysis of flood evacuation risks for priority populations in Vermont. In the event of significant inland 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 Vermont

```{r FloodData, include=FALSE}
# Analysis of flooding  evacuation risks for priority populations in Vermont
# Estimate proportion of area occupied by FEMA flood zones and river corridors
# # Read in flood hazard areas
# vt_nfhza_2852 <- st_read(dsn = "DATA/FEMA/VT/NFHL_50_20161207.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 = 2852) %>%
#   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")) %>% 
#   filter(!st_is_empty())

# # save this layer to save time
# saveRDS(vt_nfhza_2852,"DATA/FEMA/VT/vt_nfhza_2852.Rds")
# # load previously process flood layer
# vt_nfhza_2852 <- readRDS("DATA/FEMA/VT/vt_nfhza_2852.Rds")


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

# # crop flood zones to land areas only
# start_time <- Sys.time()
# vt_nfhza_2852_land <- vt_nfhza_2852 %>%
#   crop_shape(., vt_state_sf, polygon = TRUE) %>%
#   st_difference(., vt_awater_sf) %>%
#   mutate(Area = st_area(.)) %>%
#   st_make_valid() # takes about 1.5 hours
# end_time <- Sys.time()
# end_time - start_time
# # save time by just loading file
# saveRDS(vt_nfhza_2852_land, file = "DATA/FEMA/VT/vt_nfhza_2852_land.Rds")
# load cropped flood zones
vt_nfhza_2852_land <- readRDS("DATA/FEMA/VT/vt_nfhza_2852_land.Rds") %>% 
  mutate(Area = as.numeric(Area))

# read in river corridors layer. downloaded from VT Open Geodata Portal (https://geodata.vermont.gov/datasets/VTANR::river-corridors-august-2019) in December 2020. 
rc <- st_read(dsn = "DATA/FEMA/VT/WaterHydro_RiverCorridors",
              layer = "WaterHydro_RiverCorridors_poly") %>% 
  st_transform(., crs = 2852) %>% 
  st_make_valid()

# erase areas of rc that overlap with nfhza. 51.8 secs
rc_erased <- ms_erase(rc, vt_nfhza_2852_land) %>% 
  filter(!st_is_empty(.)) %>%
  st_make_valid()

# append erased rc to FEMA flood zones to get total area at flood risk. NOTE that order of rbind matters! first, need to harmonize columns so that they match.
rc_erased2 <- transmute(rc_erased, TYPE = "RC")
# harmonize nfhza variables with rc and rbind to create one combined layer
vt_flood <- vt_nfhza_2852_land %>% 
  transmute(TYPE = Interval) %>% 
  rename(geometry = SHAPE) %>% 
  rbind(rc_erased2, .) %>% 
  filter(!st_is_empty(.)) %>%
  st_make_valid()

# erase areas of FEMA flood zones that overlap with river corridors. 68 secs.
fema_erased <- ms_erase(vt_nfhza_2852_land, rc) %>% 
  filter(!st_is_empty(.)) %>%
  st_make_valid()


# Estimate populations exposed to non-overlapping areas of nfhza and rc and to combined areas
load("DATA/ne_layers.rds")

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

vt_tracts_2852 <- ne_tracts_sf %>% 
  filter(STATE == "Vermont",
         !st_is_empty(.)) %>% 
  st_transform(., crs = 2852)

# Extract state boundary
vt_state_sf <- ne_states_sf_cb %>% 
  filter(NAME == "Vermont") %>% 
  st_transform(., crs = 2852)

#### Processing in ArcGIS #####
# # Write out block groups for processing in ArcGIS
# vt_blkgrp_2852 %>%
#   dplyr::select(GEOID) %>%
#   st_write(., "DATA/FEMA/VT/vt_blkgrp_2852.shp", delete_layer = TRUE)
# #
# # repeat for tracts
# vt_tracts_2852 %>%
#   dplyr::select(GEOID) %>%
#   st_write(., "DATA/FEMA/VT/vt_tracts_2852.shp", delete_layer = TRUE)
# 
# # repeat for state boundary
# vt_state_sf %>%
#   st_write(., "DATA/FEMA/VT/vt_state_2852.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 vt_blkgrp_2852 and vt_tracts_2852 that overlap with undeveloped areas in NLCD shapefiles. Compute OldArea of erased polygons in sqm to identify area of developed polygons remaining. 
# Intersect erased vt_blkgrps and erased vt_tracts with NFHZA and Hurricane evacuation zones. Read back into R.

##### Return to working in R ######
# # read in processed vt_blkgrps and vt_tracts
# st_layers(dsn = "DATA/FEMA/VT")

vt_blkgrps_developed <- st_read(dsn = "DATA/FEMA/VT",
                                layer = "vt_blkgrps_developed") %>% 
  st_transform(., crs = 2852) %>%
  st_make_valid()

vt_tracts_developed <- st_read(dsn = "DATA/FEMA/VT",
                               layer = "vt_tracts_developed") %>%  
  st_transform(., crs = 2852) %>%
  st_make_valid()

# # calculate population within total flood area
# # # intersect total flood area with developed block groups. 24 mins.
# # system.time(vt_blkgrps_flood <- vt_blkgrps_developed %>% 
# #   transmute(GEOID = GEOID, 
# #             OldArea = as.numeric(st_area(.))) %>% 
# #   st_intersection(., vt_flood) %>% 
# #   mutate(NewArea = as.numeric(st_area(.))))
# # # save the output to save time
# # save(vt_blkgrps_flood, file = "DATA/FEMA/VT/vt_blkgrps_flood.Rds")
# load("DATA/FEMA/VT/vt_blkgrps_flood.Rds")

# calculate populations within flood risk areas
# USE PARALLEL PROCESSING
# determine number of cores to use for parallel processing
n.cores <- detectCores()-1

# create parallel cluster nodes
clust <- makeCluster(n.cores)

# use ClusterEvalQ to load needed packages in each cluster
clusterEvalQ(clust, {
  library(tidyverse)
  library(sf)
  library(tmaptools)
})

# export variables to each node for use in processing
clusterExport(clust, varlist = c("vt_flood", "vt_blkgrps_developed"))

# break up data by county for faster processing
# download counties layer
vt_counties <- counties(state = "VT", cb = TRUE) %>% 
  st_transform(.,crs=2852)

# break up counties into a list to run with parLapply
county_list <- list()
for (i in 1:(nrow(vt_counties))){
  county_list[[i]] <- vt_counties[i,]
}

# run in parallel. takes about 1.8 mins; roughly 13x faster than non-parallel (24 mins)!
vt_blkgrps_flood_list <- parLapply(clust, county_list, function(x){
    crop_shape(vt_flood, x, polygon = T) %>% 
      st_intersection(vt_blkgrps_developed, .) %>% 
      mutate(NewArea = as.numeric(st_area(.)))
  } )

stopCluster(clust)

# bring it all together
vt_blkgrps_flood <- do.call(rbind, vt_blkgrps_flood_list) %>% 
  st_make_valid()


# repeat for tracts
# create parallel cluster nodes
clust <- makeCluster(n.cores)

# use ClusterEvalQ to load needed packages in each cluster
clusterEvalQ(clust, {
  library(tidyverse)
  library(sf)
  library(tmaptools)
})

# export variables to each node for use in processing
clusterExport(clust, varlist = c("vt_flood", "vt_tracts_developed"))

# run in parallel. takes about 6.2 mins; about 8x faster than non-parallel (48 mins)
vt_tracts_flood_list <- parLapply(clust, county_list, function(x){
    crop_shape(vt_flood, x, polygon = T) %>% 
      st_intersection(vt_tracts_developed, .) %>% 
      mutate(NewArea = as.numeric(st_area(.)))
  } )

stopCluster(clust)

# bring it all together
vt_tracts_flood <- do.call(rbind, vt_tracts_flood_list) %>% 
  st_make_valid()

# allocate populations to intersected areas
vt_blkgrps_flood <- vt_blkgrps_flood %>% 
  left_join(., st_drop_geometry(vt_blkgrp_2852), by = "GEOID") %>% 
  filter(!st_is_empty(.)) %>% 
  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)

vt_tracts_flood <- vt_tracts_flood %>%
  left_join(., st_drop_geometry(vt_tracts_2852), by = "GEOID") %>% 
  mutate(Proportion = as.numeric(NewArea/OldArea),
         NewDisabled = disabledOver18E*Proportion,
         NewNoCar = HHnoCarE*Proportion)

# Compute total block group populations within river corridors
vt_flood_blkgrp_df <- vt_blkgrps_flood %>%
  st_drop_geometry() %>%
  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)

vt_flood_tracts_df <- vt_tracts_flood %>% 
  st_drop_geometry() %>% 
  summarize(`Disabled` = as.integer(sum(NewDisabled)),
            `No Car HH` = as.integer(sum(NewNoCar))) %>%
  gather(key = Group, value = FloodPop)


# calculate population within river corridors, excluding areas overlapping with FEMA flood zones
# intersect river corridors with developed block groups. about 1 min for block groups; about 3.8 mins for tracts.
vt_blkgrps_rc <- vt_blkgrps_developed %>% 
  transmute(GEOID = GEOID, 
            OldArea = as.numeric(st_area(.))) %>% 
  st_intersection(., rc_erased) %>% 
  mutate(NewArea = as.numeric(st_area(.)))

vt_tracts_rc <- vt_tracts_developed %>% 
  transmute(GEOID = GEOID, 
            OldArea = as.numeric(st_area(.))) %>% 
  st_intersection(., rc_erased) %>% 
  mutate(NewArea = as.numeric(st_area(.)))

# allocate populations to intersected areas
vt_blkgrps_rc <- vt_blkgrps_rc %>% 
  left_join(., st_drop_geometry(vt_blkgrp_2852), by = "GEOID") %>% 
  filter(!st_is_empty(.)) %>% 
  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)

vt_tracts_rc <- vt_tracts_rc %>% 
  left_join(., st_drop_geometry(vt_tracts_2852), by = "GEOID") %>% 
  mutate(Proportion = as.numeric(NewArea/OldArea),
         NewDisabled = disabledOver18E*Proportion,
         NewNoCar = HHnoCarE*Proportion)

# Compute total block group populations within river corridors
vt_rc_blkgrp_df <- vt_blkgrps_rc %>%
  st_drop_geometry() %>%
  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 = RCPop)

vt_rc_tracts_df <- vt_tracts_rc %>% 
  st_drop_geometry() %>% 
  summarize(`Disabled` = as.integer(sum(NewDisabled)),
            `No Car HH` = as.integer(sum(NewNoCar))) %>%
  gather(key = Group, value = RCPop)


# calculate population within FEMA flood zones, excluding areas overlapping with river corridors
# intersect fema flood zones with developed block groups.
# create parallel cluster nodes
clust <- makeCluster(n.cores)

# use ClusterEvalQ to load needed packages in each cluster
clusterEvalQ(clust, {
  library(tidyverse)
  library(sf)
  library(tmaptools)
})

# export variables to each node for use in processing
clusterExport(clust, varlist = c("fema_erased", "vt_blkgrps_developed"))

# run in parallel. takes about 1 min; roughly 13x faster than non-parallel (13.7 mins)!
vt_blkgrps_nfhza_list <- parLapply(clust, county_list, function(x){
    crop_shape(fema_erased, x, polygon = T) %>% 
      st_intersection(vt_blkgrps_developed, .) %>% 
      mutate(NewArea = as.numeric(st_area(.)))
  } )

stopCluster(clust)

# bring it all together
vt_blkgrps_nfhza <- do.call(rbind, vt_blkgrps_nfhza_list) %>% 
  st_make_valid()

# repeat for tracts
# create parallel cluster nodes
clust <- makeCluster(n.cores)

# use ClusterEvalQ to load needed packages in each cluster
clusterEvalQ(clust, {
  library(tidyverse)
  library(sf)
  library(tmaptools)
})

# export variables to each node for use in processing
clusterExport(clust, varlist = c("fema_erased", "vt_tracts_developed"))

# run in parallel. takes about 1.9 mins
vt_tracts_nfhza_list <- parLapply(clust, county_list, function(x){
    crop_shape(fema_erased, x, polygon = T) %>% 
      st_intersection(vt_tracts_developed, .) %>% 
      mutate(NewArea = as.numeric(st_area(.)))
  } )

stopCluster(clust)

# bring it all together
vt_tracts_nfhza <- do.call(rbind, vt_tracts_nfhza_list) %>% 
  st_make_valid()


# allocate populations to intersected areas
vt_blkgrps_nfhza <- vt_blkgrps_nfhza %>% 
  left_join(., st_drop_geometry(vt_blkgrp_2852), by = "GEOID") %>% 
  filter(!st_is_empty(.)) %>% 
  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)

vt_tracts_nfhza <- vt_tracts_nfhza %>% 
  left_join(., st_drop_geometry(vt_tracts_2852), by = "GEOID") %>% 
  filter(!st_is_empty(.)) %>% 
  mutate(Proportion = as.numeric(NewArea/OldArea),
         NewDisabled = disabledOver18E*Proportion,
         NewNoCar = HHnoCarE*Proportion)

# Compute total block group populations within fema flood zones
vt_nfhza_blkgrp_df <- vt_blkgrps_nfhza %>%
  st_drop_geometry() %>%
  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 = FEMAPop)

vt_nfhza_tracts_df <- vt_tracts_nfhza %>% 
  st_drop_geometry() %>% 
  summarize(`Disabled` = as.integer(sum(NewDisabled)),
            `No Car HH` = as.integer(sum(NewNoCar))) %>%
  gather(key = Group, value = FEMAPop)


# join the dfs and calculate populations living in areas that are both exposed to FEMA flood zones AND to river corridors and join with total state populations to calculate percentages
# Compute total tract populations within the state for same groups
vt_tract_flood_pops_df <- vt_tracts_2852 %>%
  as.data.frame() %>%
  summarize(`Disabled` = sum(disabledOver18E),
            `No Car HH` = sum(HHnoCarE)) %>%
  gather(key = Group, value = VTPop) %>%
  left_join(.,vt_nfhza_tracts_df, by = "Group") %>% 
  left_join(., vt_rc_tracts_df, by = "Group") %>% 
  left_join(., vt_flood_tracts_df, by = "Group")

# Compute populations for state,and join with flood pops
vt_FloodPops_df <- vt_blkgrp_2852 %>%
  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 = VTPop) %>%
  left_join(.,vt_nfhza_blkgrp_df, by = "Group") %>% 
  left_join(., vt_rc_blkgrp_df, by = "Group") %>% 
  left_join(., vt_flood_blkgrp_df, by = "Group") %>%
  rbind(.,vt_tract_flood_pops_df) %>%
  mutate(FEMA_RCpop = FloodPop - (FEMAPop + RCPop), .after = RCPop) %>% 
  mutate(across(FEMAPop:FloodPop, ~ .x / VTPop * 100, .names = "Pct{.col}"))

vt_towns_sf <- county_subdivisions(state = "VT", cb = TRUE) %>% 
  st_transform(., crs = 2852)

town_tracts_df <- vt_tracts_2852 %>% 
  dplyr::select(-GEOID) %>% 
  st_centroid(.) %>% 
  st_intersection(vt_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 <- vt_blkgrp_2852 %>% 
  dplyr::select(-GEOID) %>% 
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  group_by(NAME) %>% 
  summarize(GEOID = first(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")

vt_state_sf_cb <- ne_states_sf_cb %>% 
  filter(NAME == "Vermont") %>% 
  st_transform(., crs = st_crs(vt_blkgrp_2852))

# clean up
rm(list = ls(pattern = "_list"))

# # save data for later analysis
save(vt_blkgrps_flood, vt_blkgrps_nfhza, vt_blkgrps_rc, vt_tracts_nfhza, 
     vt_tracts_flood, vt_tracts_rc, vt_FloodPops_df, 
     file = "DATA/FEMA/VT/nfhza_census.Rds")
```

In comparison to other part's of New England, relatively little of Vermont's land area is within flood zones. However, a significant proportion of the population are nevertheless exposed to the risk of flooding from overbanking of inland water bodies (e.g., ponds and rivers) as well as fluvial erosion from streams and rivers. Indeed, damage surveys in Vermont have shown that fluvial erosion, not inundation, is the most common natural hazard type in Vermont.  

The analysis of flood exposure presented here is based on a combination of the Federal Emergency Management Agency's (FEMA) National Flood Hazard Layer (NFHL), a digital version of FEMA's most recent flood maps,[^NFHL] and river corridors mapped by the Vermont Agency of Natural Resources (ANR).[^ANR]  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.  The ANR's river corridors represent areas surrounding a river that allow for the meandering, floodplain, and the riparian functions necessary to restore and maintain the naturally stable or least erosive form of a river in order to minimize erosion hazards over time. Lands within and immediately abutting a river corridor are at higher risk to fluvial erosion.  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 Vermont's population lives near inland water bodies. As a result, a significant proportion of the population lives near or within known flood zones as well as areas subject to fluvial erosion. Figure \@ref(fig:mapNFHZAPop) below shows 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), as well as areas within river corridors. 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. Note that portions of northwest and northeast Vermont are not currently mapped for flooding risk by FEMA. 
```{r mapNFHZAdata, include=FALSE}
# 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.
vt_towns_sf_pts <- county_subdivisions(state = "VT", cb = TRUE) %>% 
  filter(NAME %in% c("Brattleboro",
                     "Manchester",
                     "Rutland",
                     "Middlebury",
                     "Montpelier",
                     "Burlington",
                     "Stowe",
                     "Greensboro"),
         GEOID != "5002161225") %>% # remove duplicate Rutland point
  st_transform(., crs = 2852) %>% 
  st_centroid(of_largest_polygon = TRUE)

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

vt_highways2nd <- primary_secondary_roads("VT") %>% 
  filter(FULLNAME %in% c("US Hwy 2","US Hwy 4","US Hwy 5","US Hwy 7")) %>%
  st_transform(., crs = 2852)

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

I91roadSegment <- vt_highways %>% 
  filter(LINEARID == "110344465713")

I91roadSegment2 <- vt_highways %>% 
  filter(LINEARID == "1103059124170")

I93roadSegment <- vt_highways %>% 
  filter(LINEARID == "1105281295268")

USHwy2Segment <- vt_highways %>% 
  filter(LINEARID == "1105317260813")

USHwy4Segment <- vt_highways %>% 
  filter(LINEARID == "1104258038927")

USHwy7Segment <- vt_highways %>% 
  filter(LINEARID == "1106087319624")

# 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")
I93 <- tmap_icons("https://upload.wikimedia.org/wikipedia/commons/thumb/0/0d/I-93.svg/200px-I-93.svg.png")
Hwy2 <- tmap_icons("https://upload.wikimedia.org/wikipedia/commons/thumb/2/27/US_2.svg/200px-US_2.svg.png")
Hwy4 <- tmap_icons("https://upload.wikimedia.org/wikipedia/commons/thumb/7/71/US_4.svg/200px-US_4.svg.png")
Hwy7 <- tmap_icons("https://upload.wikimedia.org/wikipedia/commons/thumb/9/95/US_7.svg/200px-US_7.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/VT/NFHL_50_20161207.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(vt_state_sf)) %>%
#   st_make_valid() %>%
#   st_difference(vt_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 9.2mins
# 
# # save time by loading previously processed file
# saveRDS(nofemapolys,"DATA/FEMA/VT/nofemapolys.Rds")
# load data
nofemapolys <- readRDS("DATA/FEMA/VT/nofemapolys.Rds")

# # Create random points, with 1 point for every 100 people
# vt_totalpop_pts <- vt_blkgrps_developed %>% 
#   select(NewPop) %>% 
#   filter(NewPop >= 100) %>% 
#   st_sample(., size = round(.$NewPop/100)) %>% # create 1 random point for every 100 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()
vt_nfhza_2852_land_simple <- vt_nfhza_2852_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 vt_awater_sf to speed up mapping
vt_awater_sf_simple <- vt_awater_sf %>% 
  st_cast(., "MULTIPOLYGON") %>% # homogenize type
  st_cast(., "POLYGON") %>% 
  st_simplify(., dTolerance = 100) %>% # reduce vertices
  st_make_valid() %>%
  group_by() %>% 
  summarize()

rc_simple <- ms_simplify(rc)
```

```{r mapNFHZAPop, fig.align = "center", fig.cap="Flood Zones and River Corridors across Vermont."}
# Map flood risk zones
m1 <- tm_shape(vt_blkgrp_2852, 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(vt_awater_sf_simple) + tm_fill(col = "#e6f3f7") +
  # tm_shape(ne_states_sf_cb) + tm_borders() +
  # tm_shape(ny_state_sf) + tm_borders() +
  # tm_shape(vt_totalpop_pts) + tm_dots(col = "darkgoldenrod3",
  #                                     labels = "1 dot = 100 people",
  #                                     alpha = 0.6) +
  tm_shape(vt_nfhza_2852_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(rc_simple) + tm_fill(col = "forestgreen", alpha = 0.5) +
  tm_shape(ne_states_sf_cb) + tm_borders(lwd = 0.2, alpha = 0.8) +
  tm_shape(ny_state_sf) + tm_borders(lwd = 0.2, alpha = 0.8) +
  tm_shape(vt_highways) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(vt_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(I91roadSegment2) +
  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(USHwy2Segment) +
  # tm_symbols(shape = Hwy2, border.lwd = NA, size = 0.2) +
  # tm_shape(USHwy4Segment) +
  # tm_symbols(shape = Hwy4, border.lwd = NA, size = 0.2) +
  # tm_shape(USHwy7Segment) +
  # tm_symbols(shape = Hwy7, border.lwd = NA, size = 0.2) +
  tm_shape(vt_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) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5,
               position = c(0.6,0.005)) +
  # tm_add_legend(type = "fill", col = "darkgoldenrod3",
  #               alpha = 0.6,
  #               border.col = "white", border.alpha = 0,
  #               labels = "1 dot = 100 people",
  #               title = "Total Population") +
  tm_add_legend(type = "fill", col = "gray",
                alpha = 0.4,
                border.col = "white", border.alpha = 0,
                labels = "No FEMA Data",
                title = "FEMA Coverage") +
  tm_add_legend(type = "fill", col = "forestgreen",
                alpha = 0.4,
                border.col = "white", border.alpha = 0,
                labels = "River Corridor",
                title = "VT Agency of Natural Resources") +
  tm_layout(title = "Flood Risk",
            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/VT/vt_flood.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/FEMA/VT/vt_flood.png")
```

Approximately 1.8% of Vermont's land area falls within FEMA flood zones, while 3.4% of the land area falls within river corridors. Approximately 50% of river corridor areas overlap with FEMA flood hazard areas. A total of 4.4% of land area is prone to FEMA flood risk or fluvial erosion within river corridors. The breakdown by type of risk is presented in Table \@ref(tab:tabNFHZAarea) below. 

```{r tabNFHZAarea, echo=FALSE, message=FALSE, warning=FALSE, fig.align = "center", fig.cap="Vermont Land Area within Flood Risk Areas."}
# Total land in VT
vt_area <- as.numeric(vt_state_sf$ALAND)

# Total and percentage of land area at flood risk - non-overlapping FEMA and RC
flood_area_df <- vt_flood %>% 
  mutate(Area = st_area(.)) %>% 
  as.data.frame() %>% 
  group_by() %>% 
  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)/vt_area*100),1)),"%")) %>% 
  mutate(`Risk Type` = "TOTAL", .before = "SqKm")

# Total and percentage of land area within river corridors
rc_area_df <- rc %>% 
  mutate(Area = st_area(.)) %>% 
  as.data.frame() %>% 
  group_by() %>% 
  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)/vt_area*100),1)),"%")) %>% 
  mutate(`Risk Type` = "River Corridor", .before = "SqKm")

# tally total area within river corridors
rc_area <- rc %>% mutate(Area = as.numeric(st_area(.))/10^6) %>% as.data.frame() %>% summarize(sum(Area)) %>% pull()
# tally area within river corridors, but only for areas where there is also FEMA data. erase river corridors from areas missing fema data and tally remaining river corridor area
rc_infema_area <- nofemapolys %>% 
  select(-Area) %>% 
  ms_erase(rc, .) %>% 
  mutate(Area = as.numeric(st_area(.))/10^6) %>% as.data.frame() %>% summarize(sum(Area)) %>% pull()
# tally area of river corridors not within FEMA flood zones, but only for areas where there is FEMA data
rc_erased_infema_area <-  nofemapolys %>% 
  select(-Area) %>% 
  ms_erase(rc_erased, .) %>% 
  mutate(Area = as.numeric(st_area(.))/10^6) %>% as.data.frame() %>% summarize(sum(Area)) %>% pull()
# calculate percentage of river corridors that overlap with FEMA flood zones for areas that contain both FEMA data and river corridor data
rc_overlap_pct <- (rc_infema_area - rc_erased_infema_area)/rc_infema_area*100

# do the inverse
# tally total area within FEMA flood zones
fema_area <- vt_nfhza_2852_land %>% mutate(AreaSqKm = as.numeric(st_area(.))/10^6) %>% as.data.frame() %>% summarize(sum(AreaSqKm)) %>% pull()
# tally area of FEMA flood zones not within river corridors
fema_erased_area <- fema_erased %>% mutate(Area = as.numeric(st_area(.))/10^6) %>% as.data.frame() %>% summarize(sum(Area)) %>% pull()
# calculate percentage of FEMA flood zones within river corridors
fema_overlap_pct <- (fema_area - fema_erased_area)/fema_area*100

# Total and percentage of land area within FEMA flood zones and river corridors
vt_nfhza_2852_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)/vt_area*100),1)),"%")) %>% 
  rename(`Risk Type` = Interval) %>% 
  rbind(rc_area_df) %>% 
  rbind(flood_area_df) %>% 
  kable(longtable = T, booktabs = T, format.args = list(big.mark = ','),
        caption = "Vermont Land within Flood Areas", digits = 0, 
        align = "r") %>% 
  kable_styling(latex_options = c("repeat_header")) %>% 
  row_spec(4, bold=T, hline_after = F) %>% 
  footnote(general = "Approximately 50% of river corridor areas overlap with FEMA flood zones.")
```

The proportion of the total state population living within areas prone to flooding or fluvial erosion exceeds the percentage of the land area within  flood- or erosion-prone areas. Moreover, exposure varies by population subgroup. Table \@ref(tab:tabNFHZApop) below shows the total and percentages of the general population and various subgroups living within flood-prone areas. Aside from the total population, the largest absolute numbers of subpopulations living within a flood-prone area are low income persons followed by those under 18 and 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="Vermont Populations Living within Flood Zones or River Corridors."}
# Show table of pops within flood zones
vt_FloodPops_df %>% 
  mutate(across(PctFEMAPop:PctFloodPop, 
                ~ paste0(as.character(round(.x,1)),"%"))) %>% 
  arrange(-FloodPop) %>% 
  mutate(Group = recode(Group, "Minority" = "People of Color",
                        "No HS Dip" = "No HS Diploma")) %>%
  select(Group, VTPop, FEMAPop, PctFEMAPop, RCPop, PctRCPop, FEMA_RCpop, 
         PctFEMA_RCpop, FloodPop, PctFloodPop) %>% 
  kable(longtable = T, booktabs = T, format.args = list(big.mark = ','),
        caption = "Vermont Populations Living within Flood Zones or River Corridors",
        digits = 0, align = "r", col.names = c("Group", "VTPop", "Total", "Pct", "Total", "Pct", "Total", "Pct", "Total", "Pct")) %>% 
  add_header_above(c(" " = 2, "FEMA Only" = 2, "RC Only" = 2,
                                 "FEMA & RC" = 2, 
                                 "All Pops at Risk" = 2)) %>% 
  kable_styling(latex_options = c("repeat_header"))
```

Some groups are disproportionately exposed to these flood hazards. While approximately 9.9% of the general population are living within FEMA flood zones or river corridors, people in low income households, adults with a disability, and other priority populations, exceed the population average living within these flood-prone areas. 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="Vermont Populations Living within Flood Zones or River Corridors."}
# Create lollipop plot of pops within flood zones
vt_FloodPops_df %>% 
  mutate(Group = recode(Group, "Minority" = "People of Color",
                        "No HS Dip" = "No HS Diploma")) %>%
  ggplot(aes(x = reorder(Group,PctFloodPop), 
             y = PctFloodPop)) +
  geom_segment(aes(x = reorder(Group,PctFloodPop), 
                   xend = reorder(Group,PctFloodPop),
                   y = vt_FloodPops_df[1,10], yend = PctFloodPop), 
               color = "skyblue") +
  geom_point(color = "blue", size = 4, alpha = 0.8) +
  coord_flip() + xlab("") + ylab("") + 
  ggtitle("Vermont Populations Living within FEMA Flood Zones\nor River Corridors") + 
  theme_light() +
  theme(panel.grid.major.y = element_blank(),
        panel.border = element_blank(),
        axis.ticks.y = element_blank()) +
  geom_text(aes(x = Group, y = PctFloodPop + 0.2 * sign(PctFloodPop), 
                label = paste0(round(PctFloodPop,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 = vt_FloodPops_df[1,10], linetype = "dashed") +
  geom_text(aes(x = "Under 18", y = 11, label = "Above state\naverage"),
            color = "gray48") +
  geom_text(aes(x = "Under 18", y = 8.5, label = "Below state\naverage"),
            color = "gray48") +
  expand_limits(y = c(7.5,11.5))

ggsave("images/VT_FEMA_graph.png")
  # geom_segment(aes(x = "Total Pop", xend = "Total Pop", y = 5.8, 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 low income persons who are disproportionately exposed to flood risk in FEMA flood zones or river corridors. Maps for other priority populations can be found in Figures \@ref(fig:mapNFHZAPopConcern3) to \@ref(fig:mapNFHZAPopConcern8) and breakdowns by municipality in Tables \@ref(tab:tabNFHZANoCar) to \@ref(tab:tabNFHZALEH) in Appendix B.

```{r mapNFHZAPopConcern, include=FALSE}
# # 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_flood_pts <- pop2points(vt_blkgrps_flood,NewEng_limit,"Limited English HH") # no blkgrp with at least 10
# 
# NoCarHH_flood_pts <- pop2points(vt_tracts_flood,NewNoCar,"No Car HH")
# 
# NoHSdip_flood_pts <- pop2points(vt_blkgrps_flood,NewLths,"No HS Dip")
# 
# LowInc_flood_pts <- pop2points(vt_blkgrps_flood,NewPov,"Low Income")
# 
# Disabled_flood_pts <- pop2points(vt_tracts_flood,NewDisabled,"Disabled")
# 
# Minor_flood_pts <- pop2points(vt_blkgrps_flood,NewMinority,"Minority")
# 
# Under5_flood_pts <- pop2points(vt_blkgrps_flood,NewUnder5,"Under 5")
# 
# Over64_flood_pts <- pop2points(vt_blkgrps_flood,NewOver64,"Over 64")
# 
# Under18_flood_pts <- pop2points(vt_blkgrps_flood,NewUnder18,"Under 18")
# 
# # Bring them together
# vt_flood_vulnerable <- rbind(NoCarHH_flood_pts,
#                              Disabled_flood_pts,
#                              LowInc_flood_pts,
#                              Minor_flood_pts,
#                              Over64_flood_pts,
#                              Under5_flood_pts,
#                              Under18_flood_pts,
#                              NoHSdip_flood_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 = "_flood_pts"))
# 
# # save time by just loading file
# saveRDS(vt_flood_vulnerable, file = "DATA/FEMA/VT/vt_flood_vulnerable.Rds")

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

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

# Map over-represented pops and flood zones
m <- tm_layout(bg.color = "#e6f3f7") +
  tm_shape(vt_blkgrp_2852, 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(vt_awater_sf_simple) + tm_fill(col = "#e6f3f7") +
  tm_shape(flood_sub) +
  tm_dots(col = "darkgoldenrod3",
          labels = "1 dot = 10 people",
          alpha = 0.8) +
  # tm_facets(by = "Group",
  #           ncol = 2,
  #           free.coords = FALSE) +
  tm_shape(vt_nfhza_2852_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(rc_simple) + tm_fill(col = "forestgreen", alpha = 0.5) +
  tm_shape(ne_states_sf_cb) + tm_borders(lwd = 0.2) +
  tm_shape(ny_state_sf) + tm_borders(lwd = 0.2) +
  tm_shape(vt_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(vt_highways2nd) + 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(I91roadSegment2) +
  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(vt_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(0.6,0.005)) +
  tm_add_legend(type = "fill", col = "darkgoldenrod3",
                alpha = 0.6,
                border.col = "white", border.alpha = 0,
                labels = "1 dot = 10 people") +
  tm_add_legend(type = "fill", col = "gray",
                alpha = 0.4,
                border.col = "white", border.alpha = 0,
                labels = "No FEMA Data",
                title = "FEMA Coverage") +
  tm_add_legend(type = "fill", col = "forestgreen",
                alpha = 0.4,
                border.col = "white", border.alpha = 0,
                labels = "River Corridor",
                title = "VT Agency of Natural Resources") +
  tm_layout(title = "Households\nwithout a Car\nLiving within\nFlood Zones or\nRiver Corridors",
            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/VT/VT_flood_vulnerable1.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/FEMA/VT/VT_flood_vulnerable1.png")
```

```{r mapNFHZAPopConcern2, fig.align = "center", fig.cap="Low Income Persons Living within Flood Zones or River Corridors."}
# create pop subset
flood_sub <- vt_flood_vulnerable %>% 
  filter(Group == "Low Income")

# Map over-represented pops and flood zones
m <- tm_layout(bg.color = "#e6f3f7") +
  tm_shape(vt_blkgrp_2852, 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(vt_awater_sf_simple) + tm_fill(col = "#e6f3f7") +
  tm_shape(flood_sub) +
  tm_dots(col = "darkgoldenrod3",
          labels = "1 dot = 10 people",
          alpha = 0.8) +
  # tm_facets(by = "Group",
  #           ncol = 2,
  #           free.coords = FALSE) +
  tm_shape(vt_nfhza_2852_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(rc_simple) + tm_fill(col = "forestgreen", alpha = 0.5) +
  tm_shape(ne_states_sf_cb) + tm_borders(lwd = 0.2) +
  tm_shape(ny_state_sf) + tm_borders(lwd = 0.2) +
  tm_shape(vt_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(vt_highways2nd) + 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(I91roadSegment2) +
  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(vt_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(0.6,0.005)) +
  tm_add_legend(type = "fill", col = "darkgoldenrod3",
                alpha = 0.6,
                border.col = "white", border.alpha = 0,
                labels = "1 dot = 10 people") +
  tm_add_legend(type = "fill", col = "gray",
                alpha = 0.4,
                border.col = "white", border.alpha = 0,
                labels = "No FEMA Data",
                title = "FEMA Coverage") +
  tm_add_legend(type = "fill", col = "forestgreen",
                alpha = 0.4,
                border.col = "white", border.alpha = 0,
                labels = "River Corridor",
                title = "VT Agency of Natural Resources") +
  tm_layout(title = "Low Income\nPersons\nLiving within\nFlood Zones or\nRiver Corridors",
            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/VT/VT_flood_vulnerable2.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/FEMA/VT/VT_flood_vulnerable2.png")
```


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

[^ANR]: See Flood Ready Vermont: River Corridors. https://floodready.vermont.gov/flood_protection/river_corridors_floodplains/river_corridors.

[^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 based on the hydrologic and topographic characteristics of an area as well as data from historic flooding.  

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.


## Vermont Agency of Natural Resources (ANR) River Corridors {-}

River corridors are areas that encompass both the channel and adjacent land of rivers and streams. The adjacent land in a river corridor is the minimal area needed to contain the meander of a river in its least erosive form. These corridors were delineated to identify areas "with the expectation that new and existing structures outside the corridor may be protected from lateral channel migration using bank stabilization practices without creating new or additional hazards. Within a river corridor, existing infrastructure and improved property may be at a heightened risk from erosion and be more likely to require river management to protect over time."[^ANR] Indeed, state officials have found that most flood-related damage in Vermont occurs outside of FEMA's Special Flood Hazard Areas and inside of river corridors, primarily as a result of fluvial erosion. 

For the purposes of this analysis, a statewide geospatial data layer of river corridors was downloaded from the Vermont Open Geodata Portal (https://geodata.vermont.gov/datasets/VTANR::river-corridors-august-2019), last updated in August 2019. 

In order to identify populations at risk within river corridors, river corridor 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 below). The population at risk within river corridors was calculated as the product of the  areal proportion of the intersecting river corridor 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 river corridor polygon, it was assumed that 10% of the population is exposed to that river corridor 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 river corridor risk. River corridor 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) %>% 
  kable(longtable = T, booktabs = T,
                    col.names = c("GRID Code","Land Use Description"),
                    caption = "NLCD 2016 Land Cover Codes and Classifications", align = "l") %>% 
  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.


## 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 %>% 
  kable(longtable = T, booktabs = T, 
                    caption = "Demographic Variables", align = "l",
                    col.names = c("Variable","Description","ACS Table ID","Geography")) %>% 
  column_spec(2:3, width = "3cm") %>%
  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 or River Corridors."}
vt_tracts_flood %>% 
  st_centroid(.) %>% 
  st_intersection(vt_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`)) %>% 
  kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Households without a Car Living in FEMA Flood Zones or River Corridors", align = "r") %>% 
  kable_styling(latex_options = c("repeat_header","striped"))
```

```{r tabNFHZALINC, fig.align = "center", fig.cap="Low income persons living within FEMA Flood Zones or River Corridors."}
vt_blkgrps_flood %>% 
  st_centroid(.) %>% 
  st_intersection(vt_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`)) %>% 
  kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Low Income Persons Living in FEMA Flood Zones or River Corridors", align = "r") %>% 
  kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapNFHZAPopConcern3, fig.align = "center", fig.cap="Disabled Persons Living within Flood Zones or River Corridors."}
# create pop subset
flood_sub <- vt_flood_vulnerable %>% 
  filter(Group == "Disabled")

# Map over-represented pops and flood zones
m <- tm_layout(bg.color = "#e6f3f7") +
  tm_shape(vt_blkgrp_2852, 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(vt_awater_sf_simple) + tm_fill(col = "#e6f3f7") +
  tm_shape(flood_sub) +
  tm_dots(col = "darkgoldenrod3",
          labels = "1 dot = 10 people",
          alpha = 0.8) +
  # tm_facets(by = "Group",
  #           ncol = 2,
  #           free.coords = FALSE) +
  tm_shape(vt_nfhza_2852_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(rc_simple) + tm_fill(col = "forestgreen", alpha = 0.5) +
  tm_shape(ne_states_sf_cb) + tm_borders(lwd = 0.2) +
  tm_shape(ny_state_sf) + tm_borders(lwd = 0.2) +
  tm_shape(vt_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(vt_highways2nd) + 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(I91roadSegment2) +
  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(vt_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(0.6,0.005)) +
  tm_add_legend(type = "fill", col = "darkgoldenrod3",
                alpha = 0.6,
                border.col = "white", border.alpha = 0,
                labels = "1 dot = 10 people") +
  tm_add_legend(type = "fill", col = "gray",
                alpha = 0.4,
                border.col = "white", border.alpha = 0,
                labels = "No FEMA Data",
                title = "FEMA Coverage") +
  tm_add_legend(type = "fill", col = "forestgreen",
                alpha = 0.4,
                border.col = "white", border.alpha = 0,
                labels = "River Corridor",
                title = "VT Agency of Natural Resources") +
  tm_layout(title = "Disabled\nPersons\nLiving within\nFlood Zones or\nRiver Corridors",
            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/VT/VT_flood_vulnerable3.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/FEMA/VT/VT_flood_vulnerable3.png")
```

```{r tabNFHZADIS, fig.align = "center", fig.cap="Disabled persons living within FEMA Flood Zones or River Corridors."}
vt_tracts_flood %>% 
  st_centroid(.) %>% 
  st_intersection(vt_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`)) %>% 
  kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Disabled Persons Living in FEMA Flood Zones or River Corridors", align = "r") %>% 
  kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapNFHZAPopConcern4, fig.align = "center", fig.cap="Persons Over 64 Living within Flood Zones or River Corridors."}
# create pop subset
flood_sub <- vt_flood_vulnerable %>% 
  filter(Group == "Over 64")

# Map over-represented pops and flood zones
m <- tm_layout(bg.color = "#e6f3f7") +
  tm_shape(vt_blkgrp_2852, 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(vt_awater_sf_simple) + tm_fill(col = "#e6f3f7") +
  tm_shape(flood_sub) +
  tm_dots(col = "darkgoldenrod3",
          labels = "1 dot = 10 people",
          alpha = 0.8) +
  # tm_facets(by = "Group",
  #           ncol = 2,
  #           free.coords = FALSE) +
  tm_shape(vt_nfhza_2852_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(rc_simple) + tm_fill(col = "forestgreen", alpha = 0.5) +
  tm_shape(ne_states_sf_cb) + tm_borders(lwd = 0.2) +
  tm_shape(ny_state_sf) + tm_borders(lwd = 0.2) +
  tm_shape(vt_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(vt_highways2nd) + 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(I91roadSegment2) +
  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(vt_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(0.6,0.005)) +
  tm_add_legend(type = "fill", col = "darkgoldenrod3",
                alpha = 0.6,
                border.col = "white", border.alpha = 0,
                labels = "1 dot = 10 people") +
  tm_add_legend(type = "fill", col = "gray",
                alpha = 0.4,
                border.col = "white", border.alpha = 0,
                labels = "No FEMA Data",
                title = "FEMA Coverage") +
  tm_add_legend(type = "fill", col = "forestgreen",
                alpha = 0.4,
                border.col = "white", border.alpha = 0,
                labels = "River Corridor",
                title = "VT Agency of Natural Resources") +
  tm_layout(title = "Persons\nOver 64\nLiving within\nFlood Zones or\nRiver Corridors",
            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/VT/VT_flood_vulnerable4.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/FEMA/VT/VT_flood_vulnerable4.png")
```

```{r tabNFHZAO64, fig.align = "center", fig.cap="Persons over 64 living within FEMA Flood Zones or River Corridors."}
vt_blkgrps_flood %>% 
  st_centroid(.) %>% 
  st_intersection(vt_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`)) %>% 
  kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Persons Over 64 Living in FEMA Flood Zones or River Corridors", align = "r") %>% 
  kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapNFHZAPopConcern5, fig.align = "center", fig.cap="People of Color Living within Flood Zones or River Corridors."}
# create pop subset
flood_sub <- vt_flood_vulnerable %>% 
  filter(Group == "Minority")

# Map over-represented pops and flood zones
m <- tm_layout(bg.color = "#e6f3f7") +
  tm_shape(vt_blkgrp_2852, 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(vt_awater_sf_simple) + tm_fill(col = "#e6f3f7") +
  tm_shape(flood_sub) +
  tm_dots(col = "darkgoldenrod3",
          labels = "1 dot = 10 people",
          alpha = 0.8) +
  # tm_facets(by = "Group",
  #           ncol = 2,
  #           free.coords = FALSE) +
  tm_shape(vt_nfhza_2852_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(rc_simple) + tm_fill(col = "forestgreen", alpha = 0.5) +
  tm_shape(ne_states_sf_cb) + tm_borders(lwd = 0.2) +
  tm_shape(ny_state_sf) + tm_borders(lwd = 0.2) +
  tm_shape(vt_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(vt_highways2nd) + 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(I91roadSegment2) +
  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(vt_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(0.6,0.005)) +
  tm_add_legend(type = "fill", col = "darkgoldenrod3",
                alpha = 0.6,
                border.col = "white", border.alpha = 0,
                labels = "1 dot = 10 people") +
  tm_add_legend(type = "fill", col = "gray",
                alpha = 0.4,
                border.col = "white", border.alpha = 0,
                labels = "No FEMA Data",
                title = "FEMA Coverage") +
  tm_add_legend(type = "fill", col = "forestgreen",
                alpha = 0.4,
                border.col = "white", border.alpha = 0,
                labels = "River Corridor",
                title = "VT Agency of Natural Resources") +
  tm_layout(title = "People\nof Color\nLiving within\nFlood Zones or\nRiver Corridors",
            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/VT/VT_flood_vulnerable5.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/FEMA/VT/VT_flood_vulnerable5.png")
```

```{r tabNFHZAMIN, fig.align = "center", fig.cap="People of Color living within FEMA Flood Zones or River Corridors."}
vt_blkgrps_flood %>% 
  st_centroid(.) %>% 
  st_intersection(vt_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`)) %>% 
  kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "People of Color Living in FEMA Flood Zones or River Corridors", align = "r", col.names = c("City/Town", "People of Color", "Pct of POC")) %>% 
  kable_styling(latex_options = c("repeat_header","striped")) 
```

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

# Map over-represented pops and flood zones
m <- tm_layout(bg.color = "#e6f3f7") +
  tm_shape(vt_blkgrp_2852, 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(vt_awater_sf_simple) + tm_fill(col = "#e6f3f7") +
  tm_shape(flood_sub) +
  tm_dots(col = "darkgoldenrod3",
          labels = "1 dot = 10 people",
          alpha = 0.8) +
  # tm_facets(by = "Group",
  #           ncol = 2,
  #           free.coords = FALSE) +
  tm_shape(vt_nfhza_2852_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(rc_simple) + tm_fill(col = "forestgreen", alpha = 0.5) +
  tm_shape(ne_states_sf_cb) + tm_borders(lwd = 0.2) +
  tm_shape(ny_state_sf) + tm_borders(lwd = 0.2) +
  tm_shape(vt_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(vt_highways2nd) + 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(I91roadSegment2) +
  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(vt_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(0.6,0.005)) +
  tm_add_legend(type = "fill", col = "darkgoldenrod3",
                alpha = 0.6,
                border.col = "white", border.alpha = 0,
                labels = "1 dot = 10 people") +
  tm_add_legend(type = "fill", col = "gray",
                alpha = 0.4,
                border.col = "white", border.alpha = 0,
                labels = "No FEMA Data",
                title = "FEMA Coverage") +
  tm_add_legend(type = "fill", col = "forestgreen",
                alpha = 0.4,
                border.col = "white", border.alpha = 0,
                labels = "River Corridor",
                title = "VT Agency of Natural Resources") +
  tm_layout(title = "Adults without a\nHigh School\nDiploma\nLiving within\nFlood Zones or\nRiver Corridors",
            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/VT/VT_flood_vulnerable6.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/FEMA/VT/VT_flood_vulnerable6.png")
```

```{r tabNFHZALTHS, fig.align = "center", fig.cap="Adults without a high school diploma living within FEMA Flood Zones or River Corridors."}
vt_blkgrps_flood %>% 
  st_centroid(.) %>% 
  st_intersection(vt_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`)) %>% 
  kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Adults wihtout a High School Diploma Living in FEMA Flood Zones or River Corridors", align = "r") %>% 
  kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapNFHZAPopConcern7, fig.align = "center", fig.cap="Persons Under 18 Living within Flood Zones or River Corridors."}
# create pop subset
flood_sub <- vt_flood_vulnerable %>% 
  filter(Group == "Under 18")

# Map over-represented pops and flood zones
m <- tm_layout(bg.color = "#e6f3f7") +
  tm_shape(vt_blkgrp_2852, 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(vt_awater_sf_simple) + tm_fill(col = "#e6f3f7") +
  tm_shape(flood_sub) +
  tm_dots(col = "darkgoldenrod3",
          labels = "1 dot = 10 people",
          alpha = 0.8) +
  # tm_facets(by = "Group",
  #           ncol = 2,
  #           free.coords = FALSE) +
  tm_shape(vt_nfhza_2852_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(rc_simple) + tm_fill(col = "forestgreen", alpha = 0.5) +
  tm_shape(ne_states_sf_cb) + tm_borders(lwd = 0.2) +
  tm_shape(ny_state_sf) + tm_borders(lwd = 0.2) +
  tm_shape(vt_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(vt_highways2nd) + 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(I91roadSegment2) +
  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(vt_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(0.6,0.005)) +
  tm_add_legend(type = "fill", col = "darkgoldenrod3",
                alpha = 0.6,
                border.col = "white", border.alpha = 0,
                labels = "1 dot = 10 people") +
  tm_add_legend(type = "fill", col = "gray",
                alpha = 0.4,
                border.col = "white", border.alpha = 0,
                labels = "No FEMA Data",
                title = "FEMA Coverage") +
  tm_add_legend(type = "fill", col = "forestgreen",
                alpha = 0.4,
                border.col = "white", border.alpha = 0,
                labels = "River Corridor",
                title = "VT Agency of Natural Resources") +
  tm_layout(title = "Persons\nUnder 18\nLiving within\nFlood Zones or\nRiver Corridors",
            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/VT/VT_flood_vulnerable7.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/FEMA/VT/VT_flood_vulnerable7.png")
```

```{r tabNFHZAU18, fig.align = "center", fig.cap="Persons under 18 living within FEMA Flood Zones or River Corridors."}
vt_blkgrps_flood %>% 
  st_centroid(.) %>% 
  st_intersection(vt_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`)) %>% 
  kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Persons Under 18 Living in FEMA Flood Zones or River Corridors", align = "r") %>% 
  kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapNFHZAPopConcern8, fig.align = "center", fig.cap="Children Under 5 Living within Flood Zones or River Corridors."}
# create pop subset
flood_sub <- vt_flood_vulnerable %>% 
  filter(Group == "Under 5")

# Map over-represented pops and flood zones
m <- tm_layout(bg.color = "#e6f3f7") +
  tm_shape(vt_blkgrp_2852, 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(vt_awater_sf_simple) + tm_fill(col = "#e6f3f7") +
  tm_shape(flood_sub) +
  tm_dots(col = "darkgoldenrod3",
          labels = "1 dot = 10 people",
          alpha = 0.8) +
  # tm_facets(by = "Group",
  #           ncol = 2,
  #           free.coords = FALSE) +
  tm_shape(vt_nfhza_2852_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(rc_simple) + tm_fill(col = "forestgreen", alpha = 0.5) +
  tm_shape(ne_states_sf_cb) + tm_borders(lwd = 0.2) +
  tm_shape(ny_state_sf) + tm_borders(lwd = 0.2) +
  tm_shape(vt_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(vt_highways2nd) + 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(I91roadSegment2) +
  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(vt_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(0.6,0.005)) +
  tm_add_legend(type = "fill", col = "darkgoldenrod3",
                alpha = 0.6,
                border.col = "white", border.alpha = 0,
                labels = "1 dot = 10 people") +
  tm_add_legend(type = "fill", col = "gray",
                alpha = 0.4,
                border.col = "white", border.alpha = 0,
                labels = "No FEMA Data",
                title = "FEMA Coverage") +
  tm_add_legend(type = "fill", col = "forestgreen",
                alpha = 0.4,
                border.col = "white", border.alpha = 0,
                labels = "River Corridor",
                title = "VT Agency of Natural Resources") +
  tm_layout(title = "Children\nUnder 5\nLiving within\nFlood Zones or\nRiver Corridors",
            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/VT/VT_flood_vulnerable8.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/FEMA/VT/VT_flood_vulnerable8.png")
```

```{r tabNFHZAU5, fig.align = "center", fig.cap="Children under 5 living within FEMA Flood Zones or River Corridors."}
vt_blkgrps_flood %>% 
  st_centroid(.) %>% 
  st_intersection(vt_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`)) %>% 
  kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Children Under 5 Living in FEMA Flood Zones or River Corridors", align = "r") %>% 
  kable_styling(latex_options = c("repeat_header","striped"))
```

```{r tabNFHZALEH, fig.align = "center", fig.cap="Limited English speaking households living within FEMA Flood Zones or River Corridors."}
vt_blkgrps_flood %>% 
  st_centroid(.) %>% 
  st_intersection(vt_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`)) %>% 
  kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Limited English Speaking Households Living in FEMA Flood Zones or River Corridors", align = "r") %>% 
  kable_styling(latex_options = c("repeat_header","striped")) 
```

No map is provided for Limited English speaking households because the numbers are too small to represent in any individual Census block group.