Evacuation_NewEngland.R

# Evacuation risks analysis. Identify the number of people within areas that are at risk of flooding and thus potentially need evacuation in the event of an emergency. 

library(tidyverse)
library(sf)
library(tmap)
library(tmaptools)
library(lwgeom)
library(tigris)
options(tigris_use_cache = TRUE, tigris_class = "sf")
# library(raster)
# library(fasterize)

load("DATA/ne_layers.rds")



# Download Census TIGERLine hydrography for RI
## First, extract list of county names to use with tigris::water
ri_counties <- counties("RI") %>% 
  pull(NAME)
# Next, download water features for each county and rbind to one layer
ri_awater_sf <- rbind_tigris(
  lapply(
    ri_counties, function(x) area_water(state = "RI", county = x)
  )
) %>% 
  st_union() %>% 
  st_as_sf() %>% 
  st_transform(., crs = 2840)

# Download linear water features for each county and rbind to one layer
ri_lwater_sf <- rbind_tigris(
  lapply(
    ri_counties, function(x) linear_water(state = "RI", county = x)
  )
) %>% 
  st_union() %>% 
  st_as_sf()

# Read in conservation areas to eliminate unpopulated from areal interpolation. See http://www.rigis.org/datasets/local-conservation-areas
ri_localCon_sf <- st_read(dsn = "DATA/FEMA/RI", 
                          layer = "Local_Conservation_Areas") %>% 
  st_transform(., crs = 2840) %>% 
  st_union() %>% 
  st_as_sf() %>% 
  st_make_valid()

ri_stateCon_sf <- st_read(dsn = "DATA/FEMA/RI", 
                          layer = "State_Conservation_Areas") %>% 
  st_transform(., crs = 2840) %>% 
  st_union() %>% 
  st_as_sf() %>% 
  st_make_valid()


# 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.
ri_towns_sf_pts <- county_subdivisions(state = "RI", cb = TRUE) %>% 
  filter(NAME %in% c("Providence", 
                     "Woonsocket", 
                     "Pawtucket", 
                     "Warwick", 
                     "Bristol", 
                     "Portsmouth", 
                     "Newport", 
                     "Narragansett", 
                     "North Kingstown", 
                     "Charlestown", 
                     "Situate",
                     "Glocester")) %>% 
  st_transform(., crs = 2840) %>% 
  st_centroid(of_largest_polygon = TRUE)


# Create road layer for context
ri_highways <- primary_roads() %>% 
  filter(FULLNAME %in% c("I- 95","I- 195","I- 295", "US Hwy 6")) %>%
  # group_by(FULLNAME) %>%
  # summarize(RTTYP = unique(RTTYP),
  #           MTFCC = unique(MTFCC)) %>%
  tmaptools::crop_shape(., ne_states_sf_cb) %>% 
  st_transform(., crs = 2840)

# Extract highway segments for labeling
I95roadSegment <- ri_highways %>% 
  filter(LINEARID == "110468245978")

I95roadSegment2 <- ri_highways %>% 
  filter(LINEARID == "1107052605232")

I295roadSegment <- ri_highways %>% 
  filter(LINEARID == "1104755623349")

I195roadSegment <- ri_highways %>% 
  filter(LINEARID == "110448466166")


# Create custom icons of highway shields
I95 <- tmap_icons(file = "https://upload.wikimedia.org/wikipedia/commons/thumb/6/61/I-95.svg/200px-I-95.svg.png")
I195 <- tmap_icons("https://upload.wikimedia.org/wikipedia/commons/thumb/f/f7/I-195.svg/200px-I-195.svg.png")
I295 <- tmap_icons("https://upload.wikimedia.org/wikipedia/commons/thumb/4/41/I-295.svg/200px-I-295.svg.png")

# Read in NFHL for RI. Data comes from FEMA. 
# List available layers in geodatabase
# st_layers("DATA/FEMA/RI/NFHL_44_20181118.gdb")
# Read in flood hazard areas
ri_fhza_2840 <- st_read(dsn = "DATA/FEMA/RI/NFHL_44_20181118.gdb", 
                        layer = "S_Fld_Haz_Ar") %>% 
  filter(FLD_ZONE != "OPEN WATER" & 
           !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 = 2840) %>% 
  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"))

# crop flood zones to land areas only
ri_state_sf <- ne_states_sf_cb %>% 
  filter(NAME == "Rhode Island")
ri_fhza_2840_land <- ri_fhza_2840 %>% 
  crop_shape(., ri_state_sf, polygon = TRUE) %>% 
  st_difference(., ri_awater_sf) %>% 
  mutate(Area = st_area(.))

# Total land area of RI within flood zones
ri_fhza_2840_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))

# Percentage of RI land within flood zones
round(as.numeric(sum(ri_fhza_2840_land$Area)/ri_state_sf$ALAND)*100,1)

# map out flood hazard areas
tmap_mode("view")
tm_shape(ri_awater_sf) + tm_fill(col = "skyblue", alpha = 0.4) +
tm_shape(ri_fhza_2840_land) + tm_fill(col = "FLD_ZONE", alpha = 0.7) +
  tm_shape(ri_state_sf) + tm_borders()

# read in hurricane evacuation zone layer
ri_hea_sf <- st_read(dsn = "DATA/FEMA/RI", 
                     layer = "Hurricane_Evacuation_Areas") %>% 
  mutate(EVAC = as.character(EVAC)) %>% 
  filter(EVAC %in% c("A","B","C")) %>% 
  st_transform(., crs = 2840) %>% 
  st_make_valid() %>% 
  mutate(Area = st_area(.))

# Total land area with hurricane evacuation zones
ri_hea_sf %>% 
  as.data.frame() %>% 
  group_by(EVAC) %>% 
  summarize(SqKm = round(as.numeric(sum(Area)/10^6),1), 
            SqMi = round(as.numeric(SqKm/2.59),1))

# Percentage of area within evac zones
round(as.numeric(sum(ri_hea_sf$Area)/ri_state_sf$ALAND)*100,1)

# map out hurricane evacuation zones
tm_shape(ri_hea_sf) + tm_fill(col = "EVAC", alpha = 0.7)

# # Map out high and moderate risk flood zones
# RI_FHZA %>% 
#   filter(FLD_ZONE != "OPEN WATER" & 
#            !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
#   tm_shape(.) + tm_fill("FLD_ZONE") +
#   tm_shape(ne_states_sf_cb) + tm_borders()
# 
# # Map out moderate risk flood zones
# RI_FHZA %>% 
#   filter(ZONE_SUBTY %in% c("0.2 PCT ANNUAL CHANCE FLOOD HAZARD","FLOODWAY")) %>%
#   mutate(FLD_ZONE = as.character(FLD_ZONE)) %>% # omit unused factor levels
#   tm_shape(.) + tm_fill("FLD_ZONE", palette = c(AE = "blue", X = "red")) +
#   tm_shape(ne_states_sf_cb) + tm_borders()
# 
# # Map out minimal risk flood zones
# RI_FHZA %>% 
#   filter(ZONE_SUBTY == "AREA OF MINIMAL FLOOD HAZARD") %>%
#   mutate(FLD_ZONE = as.character(FLD_ZONE)) %>% # omit unused factor levels
#   tm_shape(.) + tm_fill("FLD_ZONE") +
#   tm_shape(ne_states_sf_cb) + tm_borders()


# Use areal interpolation to calculate populations within flood zones. Approach follows method used by Qiang (2019) to eliminate unpopulated areas and then downscale populations to raster pixels. 
# Convert to projected local CRS EPSG:2840: NAD83(HARN) / Rhode Island
ri_blkgrp_2840 <- ne_blkgrp_sf %>% 
  filter(STATE == "Rhode Island") %>% 
  st_transform(., crs = 2840)

ri_tracts_2840 <- ne_tracts_sf %>% 
  filter(STATE == "Rhode Island") %>% 
  st_transform(., crs = 2840)

# Get rid of empty geometries
empty_geo <- st_is_empty(ri_fhza_2840)
ri_fhza_2840 <- ri_fhza_2840[!empty_geo,]
empty_geo <- st_is_empty(ri_blkgrp_2840)
ri_blkgrp_2840 <- ri_blkgrp_2840[!empty_geo,]
empty_geo <- st_is_empty(ri_tracts_2840)
ri_tracts_2840 <- ri_tracts_2840[!empty_geo,]

# Filter out block group areas overlapping with water or conservation areas so that we do not count exposed populations in areas where there are unlikely to be residents
start_time <- Sys.time()
ri_blkrps_fix <- st_difference(ri_blkgrp_2840,ri_awater_sf) %>% 
  st_difference(., ri_localCon_sf) %>% 
  st_difference(., ri_stateCon_sf)
end_time <- Sys.time()
end_time - start_time # takes about 7.5 - 9.5 mins

# repeat for tracts
start_time <- Sys.time()
ri_tracts_fix <- st_difference(ri_tracts_2840,ri_awater_sf) %>% 
  st_difference(., ri_localCon_sf) %>% 
  st_difference(., ri_stateCon_sf)
end_time <- Sys.time()
end_time - start_time # takes about 7.75 - 8.5 mins
# # Save layers since it takes a long time to run
# save(ri_blkrps_fix, ri_tracts_fix, file = "DATA/FEMA/RI/ri_fix.rds")

# read in National Land Cover Database raster data
nlcd <- raster("DATA/NLCD_2016_Land_Cover_L48_20190424.img") 

# create polygon layer in same crs as raster to crop
ri_polys <- st_transform(ri_blkgrp_2840, crs = projection(nlcd))

# create a function to reclassify NLCD pixels that are not 'developed' as NA
rc <- function(x) {
  ifelse(x >= 21 & x <= 24, 1, NA)
}

# Define the Proj.4 spatial reference 
# http://spatialreference.org/ref/epsg/2840/proj4/
newCRS <- "+proj=tmerc +lat_0=41.08333333333334 +lon_0=-71.5 +k=0.99999375 +x_0=100000 +y_0=0 +ellps=GRS80 +units=m +no_defs" 

# crop raster to region and reproject to RI state plane
start_time <- Sys.time()
ri_nlcd_2840 <- crop(nlcd, ri_polys) %>% 
  calc(., fun=rc) %>% 
  projectRaster(., crs = newCRS)
end_time <- Sys.time()
end_time - start_time 

# count number of developed NLCD cells per block group and calc pop per cell. Once we know the pop per pixel we can count the number of pixels falling within the flood zone to get the exposed population. 
start_time <- Sys.time()
ri_blkgrp_cnt <- ri_blkrps_fix %>% 
  mutate(RI_LOWINC = if_else(RI_INCOME == "I",totalpopE,0)) %>%
  mutate(RI_LOWINC = replace_na(RI_LOWINC,0)) %>%
  mutate(RI_MINORITIES = if_else(RI_MINORITY == "M", totalpopE,0)) %>%
  mutate(RI_MINORITIES = replace_na(RI_MINORITIES,0)) %>%
  dplyr::select(GEOID,
                totalpopE,
                minorityE,
                under5E,
                over64E,
                under18E,
                eng_limitE,
                num2povE,
                lthsE,
                RI_LOWINC,
                RI_MINORITIES) %>% 
  as_Spatial() %>% 
  raster::extract(ri_nlcd_2840, ., 
          fun = function(x, ...) length(na.omit(x)), 
          sp=TRUE, na.rm=TRUE, small=TRUE) %>% 
  st_as_sf() %>% 
  rename(cells = layer) %>% 
  mutate(popCell = totalpopE/cells,
         minorityCell = minorityE/cells,
         under5Cell = under5E/cells,
         over64Cell = over64E/cells,
         under18Cell = under18E/cells,
         engLimitCell = eng_limitE/cells,
         num2povCell = num2povE/cells,
         lthsCell = lthsE/cells,
         RIlowIncCell = RI_LOWINC/cells,
         RIminorityCell = RI_MINORITIES/cells)
end_time <- Sys.time()
end_time - start_time # takes about 2.5 - 3.4 mins

# repeat for tracts
start_time <- Sys.time()
ri_tract_cnt <- ri_tracts_fix %>% 
  dplyr::select(GEOID,
                disabledOver18E,
                HHnoCarE) %>% 
  as_Spatial() %>% 
  raster::extract(ri_nlcd_2840, ., 
                  fun = function(x, ...) length(na.omit(x)), 
                  sp=TRUE, na.rm=TRUE, small=TRUE) %>% 
  st_as_sf() %>% 
  rename(cells = layer) %>% 
  mutate(disabledCell = disabledOver18E/cells,
         HHnoCarCell = HHnoCarE/cells)
end_time <- Sys.time()
end_time - start_time # takes about 3 mins

# # convert block group polygons to raster and mask out NA values
# start_time <- Sys.time()
# ri_blkgrp_ras <- fasterize(sf = ri_blkgrp_cnt, raster = ri_nlcd_2840, 
#                            field = "popCell", fun = "first") %>% 
#   mask(., ri_nlcd_2840)
# end_time <- Sys.time()
# end_time - start_time

# Convert block groups to raster with pop per pixel as value
# iterate through variables in ri_blkgrp_cnt and create raster for each
# create a vector of variable names over which to iterate
fields <- names(as.data.frame(ri_blkgrp_cnt)[14:23])
fields_tracts <- names(as.data.frame(ri_tract_cnt)[6:7])

ri_blkgrp_ras_list <- lapply(fields, FUN = function(x) {
  fasterize(sf = ri_blkgrp_cnt, raster = ri_nlcd_2840, 
            field = x, fun = "first") %>% 
    mask(., ri_nlcd_2840)
})

ri_tracts_ras_list <- lapply(fields_tracts, FUN = function(x) {
  fasterize(sf = ri_tract_cnt, raster = ri_nlcd_2840, 
            field = x, fun = "first") %>% 
    mask(., ri_nlcd_2840)
})

# assign names to list elements
names(ri_blkgrp_ras_list) <- fields
names(ri_tracts_ras_list) <- fields_tracts

# map it out
tmap_mode("view")
ri_blkgrp_ras_list[[1]] %>% 
  tm_shape(.) + tm_raster(alpha = 0.7) +
  tm_shape(RI_FHZA_2840) + tm_borders(col = "blue")

tmap_mode("view")
ri_tracts_ras_list[[1]] %>% 
  tm_shape(.) + tm_raster(alpha = 0.7) +
  tm_shape(RI_FHZA_2840) + tm_borders(col = "blue")

# # remove last item in list because it's just the geometry column
# ri_blkgrp_ras_list[length(ri_blkgrp_ras_list)] <- NULL

# # extract sum of population within flood zones
# start_time <- Sys.time()
# ri_fhza_pop <- RI_FHZA_2840 %>% 
#   as_Spatial() %>% 
#   extract(ri_blkgrp_ras, ., fun = sum,
#           sp = TRUE, na.rm = TRUE, small = TRUE) %>% 
#   st_as_sf() %>% 
#   mutate(pop = layer)
# end_time <- Sys.time()
# end_time - start_time

# iterate through list and extract sum of pops within flood zone
start_time <- Sys.time()
ri_blkgrp_flood_list <- lapply(ri_blkgrp_ras_list, function(x) {
  RI_FHZA_2840 %>% 
    as_Spatial() %>% 
    extract(x, ., fun = sum,
            sp = TRUE, na.rm = TRUE, small = TRUE) %>% 
    st_as_sf()
})
end_time <- Sys.time()
end_time - start_time # takes about 46 - 51 mins

# iterate through list and extract sum of pops within hurricane evacuation zones
start_time <- Sys.time()
ri_blkgrp_hevac_list <- lapply(ri_blkgrp_ras_list, function(x) {
  ri_hea_sf %>% 
    as_Spatial() %>% 
    raster::extract(x, ., fun = sum,
            sp = TRUE, na.rm = TRUE, small = TRUE) %>% 
    st_as_sf()
})
end_time <- Sys.time()
end_time - start_time # takes about 22  mins

# repeat for tracts
start_time <- Sys.time()
ri_tracts_flood_list <- lapply(ri_tracts_ras_list, function(x) {
  RI_FHZA_2840 %>% 
    as_Spatial() %>% 
    raster::extract(x, ., fun = sum,
            sp = TRUE, na.rm = TRUE, small = TRUE) %>% 
    st_as_sf()
})
end_time <- Sys.time()
end_time - start_time # takes about 10 mins

# repeat for tracts for hurricane evacuation zones
start_time <- Sys.time()
ri_tracts_hevac_list <- lapply(ri_tracts_ras_list, function(x) {
  ri_hea_sf %>% 
    as_Spatial() %>% 
    raster::extract(x, ., fun = sum,
                    sp = TRUE, na.rm = TRUE, small = TRUE) %>% 
    st_as_sf()
})
end_time <- Sys.time()
end_time - start_time # takes about 4.5 mins

# save a copy jic
# save(ri_blkgrp_flood_list,ri_tracts_flood_list, ri_tracts_hevac_list, ri_blkgrp_hevac_list, file = "DATA/FEMA/RI/ri_flood_lists.rds")

# class(ri_blkgrp_flood_list[[1]])
# names(ri_blkgrp_flood_list)
# glimpse(as.data.frame(ri_blkgrp_flood_list[[1]]))
# # Sum of population within flood zones
# as.data.frame(ri_blkgrp_flood_list[[2]]) %>% 
#   summarize(sum(layer))
# # map it out
# tmap_mode("view")
# ri_blkgrp_flood_list[[1]] %>% 
#   tm_shape(.) + tm_fill("layer", alpha = 0.5)

# create a df of sum of flooded pops for each group
# first, convert each list item from sf to df
ri_flooded_df <- lapply(ri_blkgrp_flood_list, FUN = as.data.frame) %>% 
  # next, extract layer column from each list df. 
  lapply(., "[", "layer") %>% 
  # finally, append columns to one df
  do.call(cbind, .)
# name the columns and replace Cell with BG
names(ri_flooded_df) <- str_replace(fields, "Cell", "BG")
# Add the flood zone column
ri_flooded_df <- ri_blkgrp_flood_list[[1]] %>% 
  as.data.frame() %>% 
  dplyr::select(FLD_ZONE) %>% 
  cbind(., ri_flooded_df)

# do the same for hurricane evac data for block groups
# first, convert each list item from sf to df
ri_blkgrp_hevac_df <- lapply(ri_blkgrp_hevac_list, FUN = as.data.frame) %>% 
  # next, extract layer column from each list df. 
  lapply(., "[", "layer") %>% 
  # finally, append columns to one df
  do.call(cbind, .)
# name the columns and replace Cell with BG
names(ri_blkgrp_hevac_df) <- str_replace(fields, "Cell", "BG")
# Add the hurricane evac zone column
ri_blkgrp_hevac_df <- ri_blkgrp_hevac_list[[1]] %>% 
  as.data.frame() %>% 
  dplyr::select(EVAC) %>% 
  cbind(., ri_blkgrp_hevac_df)

# repeat for tracts
ri_floodedTracts_df <- lapply(ri_tracts_flood_list, FUN = as.data.frame) %>% 
  lapply(., "[", "layer") %>% 
  do.call(cbind, .)
# name the columns and replace Cell with BG
names(ri_floodedTracts_df) <- str_replace(fields_tracts, "Cell", "CT")
# Add the flood zone column
ri_floodedTracts_df <- ri_tracts_flood_list[[1]] %>% 
  as.data.frame() %>% 
  dplyr::select(FLD_ZONE) %>% 
  cbind(., ri_floodedTracts_df)

# do the same for hurricane evac data for tracts
ri_tracts_hevac_df <- lapply(ri_tracts_hevac_list, FUN = as.data.frame) %>% 
  lapply(., "[", "layer") %>% 
  do.call(cbind, .)
# name the columns and replace Cell with BG
names(ri_tracts_hevac_df) <- str_replace(fields_tracts, "Cell", "CT")
# Add the hurricane evac zone column
ri_tracts_hevac_df <- ri_tracts_hevac_list[[1]] %>% 
  as.data.frame() %>% 
  dplyr::select(EVAC) %>% 
  cbind(., ri_tracts_hevac_df)

# Sum total block group populations within flood zones
ri_bgsum_flooded <- ri_flooded_df %>% 
  dplyr::select(-FLD_ZONE) %>% 
  rename(`Total Pop` = popBG,
         Minority = minorityBG,
         `Under 5` = under5BG,
         `Over 64` = over64BG,
         `Under 18` = under18BG,
         `Limited English HH` = engLimitBG,
         `Low Income` = num2povBG,
         `No HS Dip` = lthsBG,
         `RI Low Income` = RIlowIncBG,
         `RI Minority` = RIminorityBG) %>%
  gather(key = "Group", value = "FloodPop") %>% 
  group_by(Group) %>% 
  summarize(FloodPop = sum(FloodPop, na.rm = TRUE))

# Sum block group populations of concern for the state
ri_blkgrp_sum <- ri_blkgrp_2840 %>% 
  as.data.frame() %>% 
  mutate(RI_LOWINC = if_else(RI_INCOME == "I",totalpopE,0)) %>%
  mutate(RI_LOWINC = replace_na(RI_LOWINC,0)) %>%
  mutate(RI_MINORITIES = if_else(RI_MINORITY == "M", totalpopE,0)) %>%
  mutate(RI_MINORITIES = replace_na(RI_MINORITIES,0)) %>%
  transmute(`Total Pop` = totalpopE,
            Minority = minorityE,
            `Under 5` = under5E,
            `Over 64` = over64E,
            `Under 18` = under18E,
            `Limited English HH` = eng_limitE,
            `Low Income` = num2povE,
            `No HS Dip` = lthsE,
            `RI Low Income` = RI_LOWINC,
            `RI Minority` = RI_MINORITIES) %>% 
  gather(., key = "Group", value = "RIPop") %>% 
  group_by(Group) %>% 
  summarize(RIPop = sum(RIPop, na.rm = TRUE))

# Sum total block group populations within hurricane evacuation zones
ri_bgsum_hevac <- ri_blkgrp_hevac_df %>% 
  dplyr::select(-EVAC) %>% 
  rename(`Total Pop` = popBG,
         Minority = minorityBG,
         `Under 5` = under5BG,
         `Over 64` = over64BG,
         `Under 18` = under18BG,
         `Limited English HH` = engLimitBG,
         `Low Income` = num2povBG,
         `No HS Dip` = lthsBG,
         `RI Low Income` = RIlowIncBG,
         `RI Minority` = RIminorityBG) %>%
  gather(key = "Group", value = "HevacPop") %>% 
  group_by(Group) %>% 
  summarize(HevacPop = sum(HevacPop, na.rm = TRUE))

# repeat for tracts
# Sum total tract populations within flood zones
ri_ctsum_flooded <- ri_floodedTracts_df %>% 
  dplyr::select(-FLD_ZONE) %>% 
  rename(Disabled = disabledCT,
         `No Car HH` = HHnoCarCT) %>%
  gather(key = "Group", value = "FloodPop") %>% 
  group_by(Group) %>% 
  summarize(FloodPop = sum(FloodPop, na.rm = TRUE))

# Sum tract populations of concern for the state
ri_tracts_sum <- ri_tracts_2840 %>% 
  as.data.frame() %>% 
  transmute(Disabled = disabledOver18E,
            `No Car HH` = HHnoCarE) %>% 
  gather(., key = "Group", value = "RIPop") %>% 
  group_by(Group) %>% 
  summarize(RIPop = sum(RIPop, na.rm = TRUE))

# Sum total tract populations within hurricane evacuation zones
ri_ctsum_hevac <- ri_tracts_hevac_df %>% 
  dplyr::select(-EVAC) %>% 
  rename(Disabled = disabledCT,
         `No Car HH` = HHnoCarCT) %>%
  gather(key = "Group", value = "HevacPop") %>% 
  group_by(Group) %>% 
  summarize(HevacPop = sum(HevacPop, na.rm = TRUE))

# join tract flooded and state pops together
ri_tract_join <- left_join(ri_tracts_sum,ri_ctsum_flooded,by="Group") %>% 
  left_join(., ri_ctsum_hevac,by="Group")

# Bring df to together and calculate percentge of population groups within flood zones and hurricane evacuation zones
ri_flooded <- left_join(ri_blkgrp_sum, ri_bgsum_flooded, by = "Group") %>% 
  left_join(., ri_bgsum_hevac, by = "Group") %>% 
  rbind(., ri_tract_join) %>% 
  mutate(PctFlooded = FloodPop/RIPop*100,
         PctHevac = HevacPop/RIPop*100)

# save for later
save(ri_flooded, file = "DATA/FEMA/RI/ri_flooded.rds")

# Create table in descending order by flooded pop
arrange(ri_flooded,desc(FloodPop))

# Create lollipop plot of flood exposed populations
ri_flooded %>% 
  ggplot(aes(x = reorder(Group,-PctFlooded), 
             y = PctFlooded)) + 
  geom_segment(aes(x = reorder(Group,-PctFlooded), 
                   xend = reorder(Group,-PctFlooded),
                   y = pull(ri_flooded[8,5]), yend = PctFlooded), 
               color = "skyblue") +
  geom_point(color = "blue", size = 4, alpha = 0.8) +
  coord_flip() + xlab("") + ylab("") + ggtitle("Rhode Island Populations within Flood Zones") + theme_light() +
  theme(panel.grid.major.y = element_blank(),
        panel.border = element_blank(),
        axis.ticks.y = element_blank()) +
  geom_text(aes(x = Group, y = PctFlooded + 0.2 * sign(PctFlooded), 
                label = paste0(round(PctFlooded,1),"%")), 
            hjust = 0.1, vjust = -0.5, size = 3,
            color=rgb(100,100,100, maxColorValue=255)) +
  scale_y_continuous(labels = function(x) paste0(x, "%")) +
  geom_hline(yintercept = pull(ri_flooded[8,5]), linetype = "dashed")

# Create lollipop plot of populations within hurricane evacuation zones
ri_flooded %>% 
  ggplot(aes(x = reorder(Group,-PctHevac), 
             y = PctHevac)) + 
  geom_segment(aes(x = reorder(Group,-PctHevac), 
                   xend = reorder(Group,-PctHevac),
                   y = pull(ri_flooded[8,6]), yend = PctHevac), 
               color = "skyblue") +
  geom_point(color = "blue", size = 4, alpha = 0.8) +
  coord_flip() + xlab("") + ylab("") + ggtitle("Rhode Island Populations within Hurricane Evacuation Zones") + theme_light() +
  theme(panel.grid.major.y = element_blank(),
        panel.border = element_blank(),
        axis.ticks.y = element_blank()) +
  geom_text(aes(x = Group, y = PctHevac + 0.2 * sign(PctHevac), 
                label = paste0(round(PctHevac,1),"%")), 
            hjust = 0.1, vjust = -0.5, size = 3,
            color=rgb(100,100,100, maxColorValue=255)) +
  scale_y_continuous(labels = function(x) paste0(x, "%")) +
  geom_hline(yintercept = pull(ri_flooded[8,6]), linetype = "dashed") +
  geom_text(aes(x = "No Car HH", y = 11.5, label = "Below state avg"),
            color = "gray48") +
  geom_segment(aes(x = "Total Pop", xend = "Total Pop", y = 13, yend = 11), 
             arrow = arrow(length = unit(0.3,"cm"))) +
  geom_text(aes(x = "Low Income", y = 14.7, label = "Above state avg"),
            color = "gray48") +
  geom_segment(aes(x = "Under 18", xend = "Under 18", y = 13.4, yend = 15), 
               arrow = arrow(length = unit(0.3,"cm")))

# Create map of pop over 64 exposed to flood hazard
# create raster of over 64 clipped to flood zones
floodOver64 <- ri_blkgrp_ras_list$over64Cell %>% 
  mask(., RI_FHZA_2840)
# map it out
tm_shape(RI_FHZA_2840) + 
  tm_fill(col = "blue", alpha = 0.3) + 
  tm_shape(floodOver64) + tm_raster(palette = "Reds",alpha = 0.7)

# Create map with blue ocean background color
tmap_mode("plot")
tm_layout(bg.color = "skyblue", saturation = .8) +
  tm_shape(ri_blkgrp_2840) + tm_fill(col = "white") +
  tm_shape(ne_states_sf_cb) + tm_fill(col="white") +
  tm_shape(ri_awater_sf) + tm_fill(col = "skyblue", alpha = 0.3) + 
  tm_shape(ri_lwater_sf) + tm_lines(col = "skyblue", alpha = 0.3) +
  tm_shape(floodOver64) + 
  tm_raster(palette = "Reds") + 
  tm_shape(ne_states_sf_cb) + tm_borders() +
  tm_shape(ri_highways) + tm_lines(col = "seashell4", lwd = 2) +
  tm_shape(I95roadSegment) + 
  tm_symbols(shape = I95, border.lwd = NA, size = 0.25) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE)

# # export block groups to shapefile
# ri_blkgrp_fix <- st_make_valid(ri_blkrps_fix)
# class(ri_blkgrp_fix)
# ri_blkgrp_fix2 <- st_cast(ri_blkgrp_fix,"MULTIPOLYGON")
# ri_blkgrp_fix2 %>% 
#   dplyr::select(-starts_with("Shape")) %>% 
#   st_write(., dsn = "DATA/FEMA/RI/ri_blkgrp_fix09.shp",
#            driver = "ESRI Shapefile", delete_layer = TRUE)
# # this works
# st_write(ri_localCon_sf, dsn = "DATA/FEMA/RI/ri_localCon_sf2.shp",
#          driver = "ESRI Shapefile")
# 
# ne_blkgrp_sf %>% 
#   dplyr::select(-starts_with("Shape")) %>% 
#   st_write(., dsn = "DATA/FEMA/RI/ne_blkgrp_sf4.shp", driver = "ESRI Shapefile", delete_layer = TRUE)
# 
# library(maptools)
# ri_blkgrp_fix2_sp <- as_Spatial(ri_blkgrp_fix2)
# rgdal::writeOGR(ri_blkgrp_fix2_sp, dsn = "DATA/FEMA/RI", 
#                 layer = "riTestsf3.shp", driver = "ESRI Shapefile",
#                 GDAL1_integer64_policy=TRUE)


###### REANALYSIS USING NLCD AS POLYGON LAYER
load("DATA/FEMA/RI/ri_fix.rds")

# Write out fixed block groups for processing in ArcGIS
ri_blkrps_fix %>%
  dplyr::select(GEOID) %>%
  st_write(., "DATA/FEMA/RI/ri_blkgrps_fix.shp", delete_layer = TRUE)

# repeat for tracts
ri_tracts_fix %>% 
  dplyr::select(GEOID) %>%
  st_write(., "DATA/FEMA/RI/ri_tracts_fix.shp", delete_layer = TRUE)


# Perform NLCD raster-to-vector conversion, vector erase/difference, and vector intersections in ArcMap because it takes too long in R. 
# Convert NLCD raster to shapefile. Isolate undeveloped areas. 
# Erase areas of ri_blkgrp_fixed and ri_tracts_fixed that overlap with undeveloped areas in NLCD shapefiles. Compute OldArea of erased polygons in sqm to identify area of developed polygons remaining. 
# Intersect ri_blkgrps and ri_tracts with NFHZA and Hurricane evacuation zones. Read back into R. 

# read in processed ri_blkgrps and ri_tracts
st_layers(dsn = "DATA/FEMA/RI")

ri_blkgrps_nfhza <- st_read(dsn = "DATA/FEMA/RI", 
                           layer = "ri_blkgrps_nfhza") %>% 
  left_join(., as.data.frame(ri_blkrps_fix), by = "GEOID") %>% 
  st_transform(., crs = 2840) %>% 
  mutate(NewArea = st_area(.)) %>%
  st_make_valid()

ri_blkgrps_hevac <- st_read(dsn = "DATA/FEMA/RI", 
                            layer = "ri_blkgrps_hevac") %>% 
  left_join(., as.data.frame(ri_blkrps_fix), by = "GEOID") %>%
  st_transform(., crs = 2840) %>% 
  mutate(NewArea = st_area(.)) %>%
  st_make_valid()

ri_tracts_nfhza <- st_read(dsn = "DATA/FEMA/RI", 
                           layer = "ri_tracts_nfhza") %>% 
  left_join(., as.data.frame(ri_tracts_fix), by = "GEOID") %>% 
  st_transform(., crs = 2840) %>% 
  mutate(NewArea = st_area(.)) %>%
  st_make_valid()

ri_tracts_hevac <- st_read(dsn = "DATA/FEMA/RI", 
                           layer = "ri_tracts_hevac") %>% 
  left_join(., as.data.frame(ri_tracts_fix), by = "GEOID") %>%
  st_transform(., crs = 2840) %>% 
  mutate(NewArea = st_area(.)) %>% 
  st_make_valid()

# Apportion populations based on geographic proportion of intersect
ri_blkgrps_nfhza <- ri_blkgrps_nfhza %>%
  mutate(RI_LOWINC = if_else(RI_INCOME == "I",totalpopE,0)) %>%
  mutate(RI_LOWINC = replace_na(RI_LOWINC,0)) %>%
  mutate(RI_MINORITIES = if_else(RI_MINORITY == "M", totalpopE,0)) %>%
  mutate(RI_MINORITIES = replace_na(RI_MINORITIES,0)) %>%
  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,
         NewRI_LOWINC = RI_LOWINC*Proportion,
         NewRI_MINORITIES = RI_MINORITIES*Proportion)

ri_blkgrps_hevac <- ri_blkgrps_hevac %>%
  mutate(RI_LOWINC = if_else(RI_INCOME == "I",totalpopE,0)) %>%
  mutate(RI_LOWINC = replace_na(RI_LOWINC,0)) %>%
  mutate(RI_MINORITIES = if_else(RI_MINORITY == "M", totalpopE,0)) %>%
  mutate(RI_MINORITIES = replace_na(RI_MINORITIES,0)) %>%
  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,
         NewRI_LOWINC = RI_LOWINC*Proportion,
         NewRI_MINORITIES = RI_MINORITIES*Proportion)

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

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

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

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

# Compute total tract populations within flood zones
ri_flood_tracts_df <- ri_tracts_nfhza %>%
  as.data.frame() %>%
  summarize(`Disabled` = as.integer(sum(NewDisabled)),
            `No Car HH` = as.integer(sum(NewNoCar))) %>%
  gather(key = Group, value = FloodPop)

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

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

ri_tract_hevac_pops_df <- ri_tracts_2840 %>%
  as.data.frame() %>%
  summarize(`Disabled` = sum(disabledOver18E),
            `No Car HH` = sum(HHnoCarE)) %>%
  gather(key = Group, value = RIPop) %>%
  left_join(.,ri_hevac_tracts_df, by = "Group")

# Compute populations for state,and join with flood pops
ri_FloodPops_df <- ri_blkgrp_2840 %>%
  as.data.frame() %>%
  mutate(RI_LOWINC = if_else(RI_INCOME == "I",totalpopE,0)) %>%
  mutate(RI_LOWINC = replace_na(RI_LOWINC,0)) %>%
  mutate(RI_MINORITIES = if_else(RI_MINORITY == "M",totalpopE,0)) %>%
  mutate(RI_MINORITIES = replace_na(RI_MINORITIES,0)) %>%
  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),
            `RI Low Income` = sum(RI_LOWINC, na.rm = TRUE),
            `RI Minority` = sum(RI_MINORITIES, na.rm = TRUE)) %>%
  gather(key = Group, value = RIPop) %>%
  left_join(., ri_flood_blkgrp_df, by = "Group") %>%
  rbind(.,ri_tract_flood_pops_df) %>%
  mutate(PctFlood = FloodPop/RIPop*100)

# Compute populations for state,  and join with hurricane evac pops
ri_HevacPops_df <- ri_blkgrp_2840 %>%
  as.data.frame() %>%
  mutate(RI_LOWINC = if_else(RI_INCOME == "I",totalpopE,0)) %>%
  mutate(RI_LOWINC = replace_na(RI_LOWINC,0)) %>%
  mutate(RI_MINORITIES = if_else(RI_MINORITY == "M",totalpopE,0)) %>%
  mutate(RI_MINORITIES = replace_na(RI_MINORITIES,0)) %>%
  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),
            `RI Low Income` = sum(RI_LOWINC, na.rm = TRUE),
            `RI Minority` = sum(RI_MINORITIES, na.rm = TRUE)) %>%
  gather(key = Group, value = RIPop) %>%
  left_join(., ri_hevac_blkgrp_df, by = "Group") %>%
  rbind(.,ri_tract_hevac_pops_df) %>%
  mutate(PctHevac = HevacPop/RIPop*100)

# Show table of pops within flood zones
ri_FloodPops_df %>% 
  arrange(-FloodPop)

# Create lollipop plot of pops within flood zones
ri_FloodPops_df %>% 
  ggplot(aes(x = reorder(Group,-PctFlood), 
             y = PctFlood)) +
  geom_segment(aes(x = reorder(Group,-PctFlood), 
                   xend = reorder(Group,-PctFlood),
                   y = ri_FloodPops_df[1,4], yend = PctFlood), 
               color = "skyblue") +
  geom_point(color = "blue", size = 4, alpha = 0.8) +
  coord_flip() + xlab("") + ylab("") + 
  ggtitle("Rhode Island Populations within Flood Zones") + theme_light() +
  theme(panel.grid.major.y = element_blank(),
        panel.border = element_blank(),
        axis.ticks.y = element_blank()) +
  geom_text(aes(x = Group, y = PctFlood + 0.2 * sign(PctFlood), 
                label = paste0(round(PctFlood,1),"%")), 
            hjust = 0.1, vjust = -0.5, size = 3,
            color=rgb(100,100,100, maxColorValue=255)) +
  scale_y_continuous(labels = function(x) paste0(x, "%")) +
  geom_hline(yintercept = ri_FloodPops_df[1,4], linetype = "dashed") +
  geom_text(aes(x = "Disabled", y = 7, label = "Below state avg"),
            color = "gray48") +
  geom_segment(aes(x = "No Car HH", xend = "No Car HH", y = 7.8, yend = 6), 
               arrow = arrow(length = unit(0.3,"cm"))) +
  geom_text(aes(x = "Low Income", y = 8.8, label = "Above state avg"),
            color = "gray48") +
  geom_segment(aes(x = "Under 18", xend = "Under 18", y = 8, yend = 9.9), 
               arrow = arrow(length = unit(0.3,"cm")))

# Map elderly populations within floodzones
# tmap_mode("view")
# ri_blkgrps_nfhza %>% 
#   filter(NewOver64 > 0) %>% 
#   tm_shape(.) + 
#   tm_fill(col = "NewOver64", palette = "Reds", alpha = 0.8)

NewOver64_sf <- ri_blkgrps_nfhza %>% 
  filter(NewOver64 >= 1)
NewOver64_sf %>% 
  tm_shape(.) + tm_borders(col = "skyblue", lwd = 0.1, alpha = 0.3) +
  tm_shape(NewOver64_sf) + 
  tm_fill("NewOver64", breaks = c(1,60,120,180,240,303), 
          style = "fixed", palette = "Reds", title = "Over 64") +
  tm_shape(ne_states_sf_cb) + tm_borders() +
  tm_shape(ri_highways) + tm_lines(col = "seashell4", lwd = 2) +
  tm_shape(I95roadSegment) + 
  tm_symbols(shape = I95, border.lwd = NA, size = 0.25) +
  tm_shape(I95roadSegment2) + 
  tm_symbols(shape = I95, border.lwd = NA, size = 0.25) +
  tm_shape(I195roadSegment) + 
  tm_symbols(shape = I195, border.lwd = NA, size = 0.25) +
  tm_shape(I295roadSegment) + 
  tm_symbols(shape = I295, border.lwd = NA, size = 0.25) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE)

# Create map with blue ocean background color
tm_layout(bg.color = "skyblue", saturation = .8) +
  tm_shape(ri_blkgrp_2840) + tm_fill(col = "white") +
  tm_shape(ne_states_sf_cb) + tm_fill(col="white") +
  tm_shape(ri_awater_sf) + tm_fill(col = "skyblue", alpha = 0.3) + 
  tm_shape(NewOver64_sf) + 
  tm_borders(col = "skyblue", lwd = 0.1, alpha = 0.3) +
  tm_shape(NewOver64_sf) + 
  tm_fill("NewOver64", breaks = c(1,60,120,180,240,303), 
          style = "fixed", palette = "Reds", title = "Over 64") + 
  tm_shape(ne_states_sf_cb) + tm_borders() +
  tm_shape(ri_highways) + tm_lines(col = "seashell4", lwd = 2) +
  tm_shape(I95roadSegment) + 
  tm_symbols(shape = I95, border.lwd = NA, size = 0.25) +
  tm_shape(I95roadSegment2) + 
  tm_symbols(shape = I95, border.lwd = NA, size = 0.25) +
  tm_shape(I195roadSegment) + 
  tm_symbols(shape = I195, border.lwd = NA, size = 0.25) +
  tm_shape(I295roadSegment) + 
  tm_symbols(shape = I295, border.lwd = NA, size = 0.25) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE)

# Create a dot density map of total populations and overlay on flood zones
# Create random points, with 1 point for every 100 people
ri_totalpop_pts <- ri_blkgrp_2840 %>% 
  select(totalpopE) %>% 
  filter(totalpopE >= 100) %>% 
  st_sample(., size = round(.$totalpopE/100)) %>% # create 1 random point for every 100 people
  st_sf(.) %>% 
  mutate(Group = "Total Pop")

# Map totalpop and flood zones
tm_layout(bg.color = "#e6f3f7") +
  tm_shape(ri_blkgrp_2840, unit = "mi") + tm_fill(col = "white") +
  tm_shape(ne_states_sf_cb) + tm_fill(col="white") +
  tm_shape(ri_awater_sf) + tm_fill(col = "#e6f3f7") + 
  tm_shape(ne_states_sf_cb) + tm_borders() +
  tm_shape(ri_totalpop_pts) + tm_dots(col = "forestgreen",
                                      labels = "1 dot = 100 persons") +
  tm_shape(ri_fhza_2840_land) + 
  tm_fill(col = "Interval", 
          palette = c("gold", "goldenrod3"), 
          labels = c("1% AEP (100-year)", "0.2% AEP (500-year)"), 
          title = "FEMA Flood Zones",
          alpha = 0.6,
          border.alpha = 0) +
  tm_shape(ri_fhza_2840_land) +
  tm_borders(col = "goldenrod1",
             lwd = 0.5,
             alpha = 0.6) +
  tm_shape(ne_states_sf_cb) + tm_borders() +
  tm_shape(ri_highways) + tm_lines(col = "seashell4", lwd = 2) +
  tm_shape(I95roadSegment) + 
  tm_symbols(shape = I95, border.lwd = NA, size = 0.25) +
  tm_shape(I95roadSegment2) + 
  tm_symbols(shape = I95, border.lwd = NA, size = 0.25) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.25) +
  tm_shape(I295roadSegment) + 
  tm_symbols(shape = I295, border.lwd = NA, size = 0.25) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_add_legend(type = "fill", col = "forestgreen",
                border.col = "white", border.alpha = 0,
                labels = "1 dot = 100 persons",
                title = "Total Population") +
  tm_layout(title = "Population Distribution \nand Flood Zones", 
            frame = TRUE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

# Create a dot density map of transit-dependent populations in flood zone
# Create random points, with 1 point for every 5 people
Over64_nfhza_pts <- ri_blkgrps_nfhza %>% 
  dplyr::select(NewOver64) %>% 
  filter(NewOver64 >= 5) %>% 
  st_sample(., size = round(.$NewOver64/5)) %>% # create 1 random point for every 5 people
  st_sf(.) %>% 
  mutate(Group = "Over 64")

Disabled_nfhza_pts <- ri_tracts_nfhza %>% 
  dplyr::select(NewDisabled) %>% 
  filter(NewDisabled >= 5) %>% 
  st_sample(., size = round(.$NewDisabled/5)) %>% 
  st_sf(.) %>% 
  mutate(Group = "Disabled")

NoCarHH_nfhza_pts <- ri_tracts_nfhza %>% 
  dplyr::select(NewNoCar) %>% 
  filter(NewNoCar >= 5) %>% 
  st_sample(., size = round(.$NewNoCar/5)) %>% 
  st_sf(.) %>% 
  mutate(Group = "No Car HH")

# Bring them together
ri_nfhza_vulnerable <- rbind(Over64_nfhza_pts,Disabled_nfhza_pts,
                             NoCarHH_nfhza_pts) %>% 
  slice(sample(1:n())) # randomise order to avoid bias in plotting order

# Map transit-dependent pops and flood zones
tm_layout(bg.color = "#e6f3f7") +
  tm_shape(ri_blkgrp_2840, unit = "mi") + tm_fill(col = "white") +
  tm_shape(ne_states_sf_cb) + tm_fill(col="white") +
  tm_shape(ri_awater_sf) + tm_fill(col = "#e6f3f7") + 
  tm_shape(ri_fhza_2840_land) + 
  tm_fill(col = "Interval", 
          palette = c("gold", "goldenrod3"), 
          labels = c("1% AEP (100-year)", "0.2% AEP (500-year)"), 
          title = "FEMA Flood Zones",
          border.alpha = 0) +
  tm_shape(ne_states_sf_cb) + tm_borders() +
  tm_shape(ri_nfhza_vulnerable) + tm_dots(col = "Group", 
                                 palette = c("green","red","blue"),
                                 legend.show = FALSE) +
  tm_shape(ne_states_sf_cb) + tm_borders() +
  tm_shape(ri_highways) + tm_lines(col = "seashell4", lwd = 2) +
  tm_shape(I95roadSegment) + 
  tm_symbols(shape = I95, border.lwd = NA, size = 0.25) +
  tm_shape(I95roadSegment2) + 
  tm_symbols(shape = I95, border.lwd = NA, size = 0.25) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.25) +
  tm_shape(I295roadSegment) + 
  tm_symbols(shape = I295, border.lwd = NA, size = 0.25) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_add_legend(type = "fill", col = c("green","red","blue"),
                border.col = "white", border.alpha = 0,
                labels = c("Disabled", "No Car HH", "Over 64"),
                title = "Population\nGroup*") +
  tm_layout(title = "Transit Dependent Populations \nwithin Flood Zones", 
            frame = TRUE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))



# Show table of pops within hurricane evacuation zones
ri_HevacPops_df %>% 
  arrange(-HevacPop)

# Create lollipop plot of pops within hurricane evac zones
ri_HevacPops_df %>% 
  ggplot(aes(x = reorder(Group,-PctHevac), 
             y = PctHevac)) +
  geom_segment(aes(x = reorder(Group,-PctHevac), 
                   xend = reorder(Group,-PctHevac),
                   y = ri_HevacPops_df[1,4], yend = PctHevac), 
               color = "skyblue") +
  geom_point(color = "blue", size = 4, alpha = 0.8) +
  coord_flip() + xlab("") + ylab("") + 
  ggtitle("Rhode Island Populations within Hurricane Evacuation Zones") +
  theme_light() +
  theme(panel.grid.major.y = element_blank(),
        panel.border = element_blank(),
        axis.ticks.y = element_blank()) +
  geom_text(aes(x = Group, y = PctHevac + 0.2 * sign(PctHevac), 
                label = paste0(round(PctHevac,1),"%")), 
            hjust = 0.1, vjust = -0.5, size = 3,
            color=rgb(100,100,100, maxColorValue=255)) +
  scale_y_continuous(labels = function(x) paste0(x, "%")) +
  geom_hline(yintercept = ri_HevacPops_df[1,4], linetype = "dashed") +
  geom_text(aes(x = "Disabled", y = 11.6, label = "Below state avg"),
            color = "gray48") +
  geom_segment(aes(x = "No Car HH", xend = "No Car HH", y = 12.8, yend = 10.8), 
               arrow = arrow(length = unit(0.3,"cm"))) +
  geom_text(aes(x = "Low Income", y = 14.5, label = "Above state avg"),
            color = "gray48") +
  geom_segment(aes(x = "Under 18", xend = "Under 18", y = 13.3, yend = 15.3), 
               arrow = arrow(length = unit(0.3,"cm")))

# create simplified version of hurricane evac polygons
ri_hea_sf_agg <- 
  ri_hea_sf %>% 
  group_by(EVAC) %>% 
  summarize()

r_hea_cata <- ri_hea_sf_agg %>% 
  filter(EVAC == "A")

r_hea_catb <- ri_hea_sf_agg %>% 
  filter(EVAC == "B")

# Map totalpop and hurricane evacuation zones
tm_layout(bg.color = "#e6f3f7") +
  tm_shape(ri_blkgrp_2840, unit = "mi") + tm_fill(col = "white") +
  tm_shape(ne_states_sf_cb) + tm_fill(col="white") +
  tm_shape(ri_awater_sf) + tm_fill(col = "#e6f3f7") + 
  tm_shape(ne_states_sf_cb) + tm_borders() +
  tm_shape(ri_totalpop_pts) + tm_dots(col = "forestgreen",
                                      labels = "1 dot = 100 persons") +
  tm_shape(ri_hea_sf_agg) + 
  tm_fill(col = "EVAC",
          palette = c("darkkhaki", "rosybrown1", "khaki"),
          labels = c("A: Catgeory 1 - 2",
                     "B: Category 3 - 4",
                     "C: Category 5"),
          title = "Evacuation Zone and\nHurricane Category",
          alpha = 0.6) +
  tm_shape(r_hea_cata) +
  tm_borders(col = "darkgoldenrod",
             lwd = 0.5,
             alpha = 0.6) +
  tm_shape(r_hea_catb) +
  tm_borders(col = "rosybrown3",
             lwd = 0.5,
             alpha = 0.6) +
  tm_shape(ne_states_sf_cb) + tm_borders() +
  tm_shape(ri_highways) + tm_lines(col = "seashell4", lwd = 2) +
  tm_shape(I95roadSegment) + 
  tm_symbols(shape = I95, border.lwd = NA, size = 0.25) +
  tm_shape(I95roadSegment2) + 
  tm_symbols(shape = I95, border.lwd = NA, size = 0.25) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.25) +
  tm_shape(I295roadSegment) + 
  tm_symbols(shape = I295, border.lwd = NA, size = 0.25) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_add_legend(type = "fill", col = "forestgreen",
                border.col = "white", border.alpha = 0,
                labels = "1 dot = 100 persons",
                title = "Total Population") +
  tm_layout(title = "Population Distribution \nand Hurricane\nEvacuation Zones", 
            frame = TRUE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

# Create a dot density of transit-dependent populations within hurricane evacuation zones
# Create random points, with 1 point for every 5 people
Over64_hevac_pts <- ri_blkgrps_hevac %>% 
  dplyr::select(NewOver64) %>% 
  filter(NewOver64 >= 5) %>% 
  st_sample(., size = round(.$NewOver64/5)) %>% # create 1 random point for every 5 people
  st_sf(.) %>% 
  mutate(Group = "Over 64")

Disabled_hevac_pts <- ri_tracts_hevac %>% 
  dplyr::select(NewDisabled) %>% 
  filter(NewDisabled >= 5) %>% 
  st_sample(., size = round(.$NewDisabled/5)) %>% 
  st_sf(.) %>% 
  mutate(Group = "Disabled")

NoCarHH_hevac_pts <- ri_tracts_hevac %>% 
  dplyr::select(NewNoCar) %>% 
  filter(NewNoCar >= 5) %>% 
  st_sample(., size = round(.$NewNoCar/5)) %>% 
  st_sf(.) %>% 
  mutate(Group = "No Car HH")

# Bring them together
ri_hevac_vulnerable <- rbind(Over64_hevac_pts,Disabled_hevac_pts,
                             NoCarHH_hevac_pts) %>% 
  slice(sample(1:n())) # randomise order to avoid bias in plotting order

# Map them out
tm_layout(bg.color = "#e6f3f7") +
  tm_shape(ri_blkgrp_2840, unit = "mi") + tm_fill(col = "white") +
  tm_shape(ne_states_sf_cb) + tm_fill(col="white") +
  tm_shape(ri_awater_sf) + tm_fill(col = "#e6f3f7") + 
  tm_shape(ri_hea_sf) + 
  tm_fill(col = "EVAC",
          palette = c("darkkhaki", "rosybrown1", "khaki"),
          labels = c("A: Catgeory 1 - 2",
                     "B: Category 3 - 4",
                     "C: Category 5"),
          title = "Evacuation Zone and\nHurricane Category") +
  tm_shape(ne_states_sf_cb) + tm_borders() +
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ri_hevac_vulnerable) + tm_dots(col = "Group", 
                                 palette = c("green","red","blue"),
                                 legend.show = FALSE) +
  tm_shape(ne_states_sf_cb) + tm_borders() +
  tm_shape(ri_highways) + tm_lines(col = "seashell4", lwd = 2) +
  tm_shape(I95roadSegment) + 
  tm_symbols(shape = I95, border.lwd = NA, size = 0.25) +
  tm_shape(I95roadSegment2) + 
  tm_symbols(shape = I95, border.lwd = NA, size = 0.25) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.25) +
  tm_shape(I295roadSegment) + 
  tm_symbols(shape = I295, border.lwd = NA, size = 0.25) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_add_legend(type = "fill", col = c("green","red","blue"),
                border.col = "white", border.alpha = 0,
                labels = c("Disabled", "No Car HH", "Over 64"),
                title = "Population\nGroup*") +
  tm_layout(title = "Transit Dependent Populations \nwithin Hurricane Evacuation\nZones", 
            frame = TRUE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))


# Intersect block groups with flood zones and allocate populations by proportion of area intersected
# ri_blkgrps_fhza_sf <- ri_blkgrps_nlcd %>%
#   mutate(OldArea = st_area(.)) %>%
#   st_intersection(RI_FHZA_2840,.) %>%
#   mutate(NewArea = st_area(.),
#          Proportion = NewArea/OldArea,
#          NewPop = as.integer(totlppE*Proportion),
#          NewMinority = as.integer(minrtyE*Proportion),
#          NewUnder5 = as.integer(under5E*Proportion),
#          NewOver64 = as.integer(over64E*Proportion),
#          NewUnder18 = as.integer(undr18E*Proportion),
#          NewEng_limit = as.integer(eng_lmE*Proportion),
#          NewPov = as.integer(num2pvE*Proportion),
#          NewLths = as.integer(lthsE*Proportion),
#          NewRI_LOWINC = as.integer(RI_LOWI*Proportion),
#          NewRI_MINORITIES = as.integer(RI_MINO*Proportion))
# 
# # Repeat for tracts
# ri_flood_tracts_sf <- ri_tracts_2840 %>% 
#   select(GEOID,
#          disabledOver18E,
#          HHnoCarE) %>% 
#   mutate(OldArea = st_area(.)) %>% 
#   st_intersection(RI_FHZA_2840,.) %>% 
#   mutate(NewArea = st_area(.),
#          Proportion = NewArea/OldArea,
#          NewDisabled = as.integer(disabledOver18E*Proportion),
#          NewNoCar = as.integer(HHnoCarE*Proportion))
# 
# # Compute total block group populations within flood zones
# ri_flood_blkgrp_df <- ri_flood_sf %>% 
#   as.data.frame() %>% 
#   summarize(`Total Pop` = sum(NewPop),
#             Minority = sum(NewMinority),
#             `Under 5` = sum(NewUnder5),
#             `Over 64` = sum(NewOver64),
#             `Under 18` = sum(NewUnder18),
#             `Limited English HH` = sum(NewEng_limit),
#             `Low Income` = sum(NewPov),
#             `No HS Dip` = sum(NewLths),
#             `RI Low Income` = sum(NewRI_LOWINC),
#             `RI Minority` = sum(NewRI_MINORITIES)) %>% 
#   gather(key = Group, value = FloodPop)
# 
# # Compute total tract populations within flood zones
# ri_flood_tracts_df <- ri_flood_tracts_sf %>% 
#   as.data.frame() %>% 
#   summarize(`Disabled` = sum(NewDisabled),
#             `No Car HH` = sum(NewNoCar)) %>% 
#   gather(key = Group, value = FloodPop)
# # Compute total tract populations within the state for same groups
# ri_tract_pops_df <- ri_tracts_2840 %>% 
#   as.data.frame() %>% 
#   summarize(`Disabled` = sum(disabledOver18E),
#             `No Car HH` = sum(HHnoCarE)) %>% 
#   gather(key = Group, value = RIPop) %>% 
#   left_join(.,ri_flood_tracts_df, by = "Group")
# 
# # Compute populations for state, join with flood pops, and create plot
# ri_FloodPops_df <- ri_blkgrp_2840 %>% 
#   as.data.frame() %>% 
#   mutate(RI_LOWINC = if_else(RI_INCOME == "I",totalpopE,0)) %>% 
#   mutate(RI_LOWINC = replace_na(RI_LOWINC,0)) %>% 
#   mutate(RI_MINORITIES = if_else(RI_MINORITY == "M",totalpopE,0)) %>% 
#   mutate(RI_MINORITIES = replace_na(RI_MINORITIES,0)) %>% 
#   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),
#             `RI Low Income` = sum(RI_LOWINC, na.rm = TRUE),
#             `RI Minority` = sum(RI_MINORITIES, na.rm = TRUE)) %>% 
#   gather(key = Group, value = RIPop) %>% 
#   left_join(., ri_flood_blkgrp_df, by = "Group") %>% 
#   rbind(.,ri_tract_pops_df) %>% 
#   mutate(PctFlood = FloodPop/RIPop*100)
# 
# # glimpse(ri_FloodPops_df)
# ri_FloodPops_df

# Create lollipop plot
ri_FloodPops_df %>% 
  ggplot(aes(x = reorder(Group,-PctFlood), 
             y = PctFlood)) +
  geom_segment(aes(x = reorder(Group,-PctFlood), 
                   xend = reorder(Group,-PctFlood),
                   y = ri_FloodPops_df[1,4], yend = PctFlood), 
               color = "skyblue") +
  geom_point(color = "blue", size = 4, alpha = 0.8) +
  coord_flip() + xlab("") + ylab("") + ggtitle("Rhode Island Populations within Flood Zones") + theme_light() +
  theme(panel.grid.major.y = element_blank(),
        panel.border = element_blank(),
        axis.ticks.y = element_blank()) +
  geom_text(aes(x = Group, y = PctFlood + 0.2 * sign(PctFlood), 
                label = paste0(round(PctFlood,0),"%")), 
            hjust = 0.1, vjust = -0.5, size = 3,
            color=rgb(100,100,100, maxColorValue=255)) +
  scale_y_continuous(labels = function(x) paste0(x, "%")) 

# Map Populations within flood zones
ri_flood_sf <- ri_flood_sf %>% 
  mutate(PopPerKm2 = NewPop/(as.numeric(NewArea)/10^6))

ri_flood_sf %>% 
  tm_shape(.) + tm_borders(col = "skyblue", lwd = 0.1, alpha = 0.3) +
  tm_shape(ri_flood_sf) + tm_fill("PopPerKm2", style = "quantile") +
  tm_shape(ne_states_sf_cb) + tm_borders() +
  tm_shape(ri_highways) + tm_lines(col = "yellow", lwd = 2) +
  tm_shape(I95roadSegment) + 
  tm_symbols(shape = I95, border.lwd = NA, size = 0.25) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE)

# Map Disabled Populations within flood zones
ri_flood_tracts_sf <- ri_flood_tracts_sf %>% 
  mutate(DisabledPerKm2 = NewDisabled/(as.numeric(NewArea)/10^6),
         NoCarPerKm2 = NewNoCar/(as.numeric(NewArea)/10^6))
# Create map with blue ocean background color
tm_layout(bg.color = "skyblue", saturation = .8) +
  tm_shape(ri_blkgrp_2840) + tm_fill(col = "white") +
  tm_shape(ne_states_sf_cb) + tm_fill(col="white") +
  tm_shape(ri_awater_sf) + tm_fill(col = "skyblue", alpha = 0.3) + 
  tm_shape(ri_lwater_sf) + tm_lines(col = "skyblue", alpha = 0.3) +
  tm_shape(ri_flood_tracts_sf) + 
  tm_borders(col = "skyblue", lwd = 0.1, alpha = 0.3) +
  tm_shape(ri_flood_tracts_sf) + 
  tm_fill("DisabledPerKm2", style = "quantile") + 
  tm_shape(ne_states_sf_cb) + tm_borders() +
  tm_shape(ri_highways) + tm_lines(col = "seashell4", lwd = 2) +
  tm_shape(I95roadSegment) + 
  tm_symbols(shape = I95, border.lwd = NA, size = 0.25) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE)

# # Extract bbox of RI for WFS request from FEMA
# ne_states_sf_cb %>% 
#   filter(NAME == "Rhode Island") %>% 
#   st_bbox()
# 
# # Extract smaller bbox
# ne_towns_sf %>% 
#   filter(STATE == "Rhode Island" & 
#            NAME == "Providence city, Providence County, Rhode Island") %>% 
#   st_bbox()
# 
# # Download WFS of flood hazard regions for RI
# RI_fhaz <- st_read("https://hazards.fema.gov/gis/nfhl/services/public/NFHL/MapServer/WFSServer?SERVICE=WFS&VERSION=1.1.0&REQUEST=GetFeature&TYPENAME=public_NFHL:Flood_Hazard_Zones&BBOX=41.77238,-71.36969,41.86182,-71.47267,urn:ogc:def:crs:EPSG:6.9:4269&MAXFEATURES=1000")
# 
# RI_fhaz <- st_read("https://hazards.fema.gov/gis/nfhl/services/public/NFHL/MapServer/WFSServer?SERVICE=WFS&REQUEST=GetFeature&TYPENAME=public_NFHL:Flood_Hazard_Zones&BBOX=41,-121.5,41.125,-121.375,urn:ogc:def:crs:EPSG::4269&MAXFEATURES=1000")
# 
# RI_fhaz <- st_read("http://ovc.catastro.meh.es/INSPIRE/wfsCP.aspx?service=WFS&request=GetFeature&Typenames=cp.cadastralparcel&SRSname=EPSG::25830&bbox=233673,4015968,233761,4016008&version=2.0.0&")