LST_NewEnglandPDF.Rmd

---
title: "Heat Risk in New England"
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)
```
\pagebreak
# Heat risk in New England

Heat, or hot weather, is the leading weather-related cause of death in the U.S.[^NWShazards] Between 1999 and 2010, the U.S. Centers for Disease Control recorded over 8,000 heat-related deaths in the U.S.[^CDCpicture] Heat-related hospitalizations or emergency department visits are estimated to be at least 10 times the death rate.[^Samuelson2020] Heat exhaustion and heat stroke are the most serious heat-related illnesses. Exposure to excessive heat can directly or indirectly cause some illnesses and can exacerbate many preexisting conditions, such as heart and lung disease. People at greatest risk for heat-related illness include children under 5 years of age, people 65 years of age and older, people who are overweight or have existing medical conditions, such as diabetes and heart disease, people who are socially isolated, and low income individuals.[^CDCpicture] 

Risk of exposure to excessive heat tends to be higher for people who work outside (e.g., agriculture, construction, and landscaping), and for those living in densely developed urban areas where there is a dearth of vegetation and an abundance of dense materials such as asphalt and concrete that absorb heat and release it more slowly (i.e., urban heat island effect). However, risk of heat-related illness or death varies within urban environments due to a variety of mitigating factors, especially age, race, and wealth.[^Williams2020] Exposure to higher temperatures than a population is accustomed to can also make people more vulnerable to heat-related illnesses and death. Although New England experiences a generally cool climate, research has shown that people in this region are actually more sensitive to elevated temperatures than those living in warmer regions of the country.[^Hondula2015]

Climate change is increasing average global temperatures, as well as the frequency, duration, and severity of extreme heat events.[^IPCCsummary] In the contiguous U.S., annual average temperatures have increased by 1.2°F (0.7°C) over the last few decades and by 1.8°F (1.0°C) relative to the beginning of the 20th century. These changes are happening more quickly in the Northeast. Average annual temperatures have increased by about 3°F (1.7°C) or more in New England since 1901. By 2050, average annual temperatures in the Northeast are projected to increase by up to 5.1°F (2.8°C) relative to the period 1975–2005, with several more days of extreme heat occurring throughout the region each year.[^4thNatlAssessment] The combination of rapidly increasing temperatures and higher sensitivity puts people in New England at especially high risk of heat-related illness and death. 

The heat risk analysis presented here is based on Land Surface Temperatures (LST) derived from NASA's Moderate-resolution Imaging Spectroradiometer (MODIS) satellite sensor.[^MODIS11] Unlike ambient air temperature measured by ground-based weather stations, LST is a measure of the radiant energy (i.e., emissivity) of the ground or surface. LST values tend to be slightly higher than ambient air temperatures, but are highly correlated with air temperatures. Abundant research has shown that satellite-derived LST is an acceptable proxy for air temperature.[^Kloog2014] Moreover, while weather stations are sparsely and unevenly distributed, satellite-derived LST has the advantage of providing an unbroken and continuous snapshot of temperature at high resolution across the region at any given moment. 

The MODIS data used here covers average day and night LST for an 8-day period from July 28 to August 4, 2019 for the entire region. That week constituted the conclusion of a historically warm July for New England This data originates at 1km spatial resolution and was extracted to average temperatures at each Census Block Group separately for day-night averages, daytime, nighttime, and urban heat island effects. LST data was validated against ground weather station data for the same time periods (see Appendix A).

```{r dataLST, include=FALSE}
library(tmap)
library(tidyverse)
library(sf)
library(tmaptools)
library(maptools)
library(ggcorrplot)
library(spdep)
library(rgdal)
library(gdalUtils)
library(raster)
library(tigris)
options(tigris_use_cache = TRUE, tigris_class = "sf")
library(lwgeom)
library(OpenStreetMap)

# load basic layers
load("DATA/ne_layers.rds")

# # create OSM basemap. see help on openmap for type options.
# ne_bmap <- ne_states_sf_cb %>% 
#   st_transform(., crs = 4326) %>% 
#   st_bbox(.) %>% 
#   read_osm(., type = "stamen-watercolor")
# 
# # map basemap
# tm_shape(ne_bmap) + tm_rgb() + tm_shape(ne_states_sf_cb) + tm_borders()

# # STEPS FOR PROCESSING LST DATA
# # MODIS data acquired through USGS EarthExplorer portal. https://earthexplorer.usgs.gov/
# # read in HDF (Hierarchical Data Format) scientific data sets (SDSs) for MODIS from Terra satellite
# sdsJuly282019h12Terra <- get_subdatasets("DATA/LST/MODISJulyAug2019/MOD11A2.A2019209.h12v04.006.2019218045630.hdf")
# sdsJuly282019h13Terra <- get_subdatasets("DATA/LST/MODISJulyAug2019/MOD11A2.A2019209.h13v04.006.2019218045706.hdf")
# 
# # sdsJuly282019h12 <- get_subdatasets("LST/MODIS8dayLST1km/2345441438/MOD11A2_A2019209_h12v04_006_2019218045630_HEGOUT.hdf")
# # sdsJuly282019h13 <- get_subdatasets("LST/MODIS8dayLST1km/2345441418/MOD11A2_A2019209_h13v04_006_2019218045706_HEGOUT.hdf")
# 
# # read in HDF (Hierarchical Data Format) scientific data sets (SDSs) for MODIS  from Aqua satellite
# sdsJuly282019h12Aqua <- get_subdatasets("DATA/LST/MODISJulyAug2019/MYD11A2.A2019209.h12v04.006.2019218030935.hdf")
# sdsJuly282019h13Aqua <- get_subdatasets("DATA/LST/MODISJulyAug2019/MYD11A2.A2019209.h13v04.006.2019218033834.hdf")
# 
# # retrieve desired data from MODIS Terra (late morning 10:06am - 12:18pm and early evening 8:42pm - 11pm)
# LSTdayJuly282019h12Terra <- readGDAL(sdsJuly282019h12Terra[1])
# daytimeH12Terra <- readGDAL(sdsJuly282019h12Terra[3])
# LSTnightJuly282019h12Terra <- readGDAL(sdsJuly282019h12Terra[5])
# nighttimeH12Terra <- readGDAL(sdsJuly282019h12Terra[7])
# LSTdayJuly282019h13Terra <- readGDAL(sdsJuly282019h13Terra[1])
# daytimeH13Terra <- readGDAL(sdsJuly282019h13Terra[3])
# LSTnightJuly282019h13Terra <- readGDAL(sdsJuly282019h13Terra[5])
# nighttimeH13Terra <- readGDAL(sdsJuly282019h13Terra[7])
# 
# # retrieve desired data sets from MODIS Aqua (early afternoon 11:48am - 2pm and later evening 12am - 3:06am)
# LSTdayJuly282019h12Aqua <- readGDAL(sdsJuly282019h12Aqua[1])
# daytimeH12Aqua <- readGDAL(sdsJuly282019h12Aqua[3])
# LSTnightJuly282019h12Aqua <- readGDAL(sdsJuly282019h12Aqua[5])
# nighttimeH12Aqua <- readGDAL(sdsJuly282019h12Aqua[7])
# LSTdayJuly282019h13Aqua <- readGDAL(sdsJuly282019h13Aqua[1])
# daytimeH13Aqua <- readGDAL(sdsJuly282019h13Aqua[3])
# LSTnightJuly282019h13Aqua <- readGDAL(sdsJuly282019h13Aqua[5])
# nighttimeH13Aqua <- readGDAL(sdsJuly282019h13Aqua[7])
# 
# # assign to raster image to manipulate data for visualization and analysis and convert from Kelvin to Fahrenheit
# LSTdayJuly282019h12r <- raster(LSTdayJuly282019h12Aqua)*(9/5) - 459.67
# LSTnightJuly282019h12r <- raster(LSTnightJuly282019h12Aqua)*(9/5) - 459.67
# LSTdayJuly282019h13r <- raster(LSTdayJuly282019h13Aqua)*(9/5) - 459.67
# LSTnightJuly282019h13r <- raster(LSTnightJuly282019h13Aqua)*(9/5) - 459.67
# 
# # # If origins are different, need to adjust tolerance to mosaic
# # origin(LSTdayJuly282019h12r)
# # origin(LSTdayJuly282019h13r)
# 
# # mosic rasters
# # create LST day mosaic raster
# LSTdayJuly2019r <- mosaic(LSTdayJuly282019h12r,LSTdayJuly282019h13r,
#                           fun=mean, tolerance = 1)
# 
# # create LST night mosaic raster
# LSTnightJuly2019r <- mosaic(LSTnightJuly282019h12r,LSTnightJuly282019h13r,
#                           fun=mean, tolerance = 1)
# 
# # create average of day and night
# LSTavgJuly2019r <- (LSTdayJuly2019r + LSTnightJuly2019r)/2
# 
# # create average day - night difference
# LSTDayNighDiffJuly2019r <- LSTdayJuly2019r - LSTnightJuly2019r
# 
# # clean up
# rm(list = ls(pattern = "h12|h13"))
# 
# # Crop LST to New England
# # define new projection
# newproj <- st_crs(ne_blkgrp_sf)[[2]]
# 
# LSTavgJuly2019r_ne <- LSTavgJuly2019r %>%
#   projectRaster(., crs = newproj) %>%
#   crop(., ne_states_sf_cb) %>%
#   mask(., ne_states_sf_cb)
# 
# LSTdayJuly2019r_ne <- LSTdayJuly2019r %>%
#   projectRaster(., crs = newproj) %>%
#   crop(., ne_states_sf_cb) %>%
#   mask(., ne_states_sf_cb)
# 
# LSTnightJuly2019r_ne <- LSTnightJuly2019r %>%
#   projectRaster(., crs = newproj) %>%
#   crop(., ne_states_sf_cb) %>%
#   mask(., ne_states_sf_cb)
# 
# LSTDayNighDiffJuly2019r_ne <- LSTDayNighDiffJuly2019r %>%
#   projectRaster(., crs = newproj) %>%
#   crop(., ne_states_sf_cb) %>%
#   mask(., ne_states_sf_cb)
# 
# # download Census-designated urban areas to contrast with rural areas
# us_urban_sf <- urban_areas(cb = TRUE)
# 
# # create urban area polygon for New England
# ne_urban_sf <- us_urban_sf %>%
#   st_centroid(.) %>%
#   st_join(., ne_states_sf_cb, join = st_within, left = FALSE) %>%
#   as.data.frame() %>%
#   dplyr::select(GEOID10) %>%
#   inner_join(us_urban_sf,.,by="GEOID10") %>%
#   group_by() %>%
#   summarize()
# 
# # create rural area polygon
# ne_rural_sf <- ne_states_sf_cb %>%
#   group_by() %>%
#   summarize() %>%
#   st_difference(., ne_urban_sf) %>%
#   st_make_valid()
# 
# # # select blkgrps that do not intersect wtih urban polygons - rural blkgrps
# # ne_blkgrpLST_rural <- st_disjoint(ne_blkgrp_sf,
# #                                   ne_urban_sf,sparse = FALSE) %>%
# #   apply(., 1, all) %>% # convert dense matrix to logical vector
# #   ne_blkgrpLST_sf[.,] # subset blkgrps where there is NO intersection
# #
# # # select blkgrps that intersect with urban polygons - urban blkgrps
# # ne_blkgrpLST_urban <- st_intersects(ne_blkgrp_sf,
# #                                     ne_urban_sf,sparse = FALSE) %>%
# #   apply(., 1, any) %>% # convert dense matrix to logical vector
# #   ne_blkgrpLST_sf[.,] # subset blkgrps where there is ANY intersection
# #
# # ne_blkgrp_urban_simple <- ne_blkgrpLST_urban %>%
# #   group_by() %>%
# #   summarize()
# #
# # ne_blkgrp_rural_simple <- ne_blkgrpLST_rural %>%
# #   group_by() %>%
# #   summarize()
# 
# # crop LST to urban
# LSTavgJuly2019r_neUrban <- LSTavgJuly2019r %>%
#   projectRaster(., crs = newproj) %>%
#   crop(., ne_urban_sf) %>%
#   mask(., ne_urban_sf)
# 
# LSTdayJuly2019r_neUrban <- LSTdayJuly2019r %>%
#   projectRaster(., crs = newproj) %>%
#   crop(., ne_urban_sf) %>%
#   mask(., ne_urban_sf)
# 
# LSTnightJuly2019r_neUrban <- LSTnightJuly2019r %>%
#   projectRaster(., crs = newproj) %>%
#   crop(., ne_urban_sf) %>%
#   mask(., ne_urban_sf)
# 
# LSTDayNighDiffJuly2019r_neUrban <- LSTDayNighDiffJuly2019r %>%
#   projectRaster(., crs = newproj) %>%
#   crop(., ne_urban_sf) %>%
#   mask(., ne_urban_sf)
# 
# # crop LST to rural
# LSTavgJuly2019r_neRural <- LSTavgJuly2019r %>%
#   projectRaster(., crs = newproj) %>%
#   crop(., ne_rural_sf) %>%
#   mask(., ne_rural_sf)
# 
# LSTdayJuly2019r_neRural <- LSTdayJuly2019r %>%
#   projectRaster(., crs = newproj) %>%
#   crop(., ne_rural_sf) %>%
#   mask(., ne_rural_sf)
# 
# LSTnightJuly2019r_neRural <- LSTnightJuly2019r %>%
#   projectRaster(., crs = newproj) %>%
#   crop(., ne_rural_sf) %>%
#   mask(., ne_rural_sf)
# 
# LSTDayNighDiffJuly2019r_neRural <- LSTDayNighDiffJuly2019r %>%
#   projectRaster(., crs = newproj) %>%
#   crop(., ne_rural_sf) %>%
#   mask(., ne_rural_sf)
# 
# # Identify stats to compute thresholds
# LSTavgJuly2019mean <- cellStats(x = LSTavgJuly2019r_ne, stat = "mean")
# LSTavgJuly2019Urbanmean <- cellStats(x = LSTavgJuly2019r_neUrban, stat = "mean")
# LSTavgJuly2019Ruralmean <- cellStats(x = LSTavgJuly2019r_neRural, stat = "mean")
# LSTdayJuly2019mean <- cellStats(x = LSTdayJuly2019r_ne, stat = "mean")
# LSTdayJuly2019Urbanmean <- cellStats(x = LSTdayJuly2019r_neUrban, stat = "mean")
# LSTdayJuly2019Ruralmean <- cellStats(x = LSTdayJuly2019r_neRural, stat = "mean")
# LSTnightJuly2019mean <- cellStats(x = LSTnightJuly2019r_ne, stat = "mean")
# LSTnightJuly2019Urbanmean <- cellStats(x = LSTnightJuly2019r_neUrban, stat = "mean")
# LSTnightJuly2019Ruralmean <- cellStats(x = LSTnightJuly2019r_neRural, stat = "mean")
# LSTavgJuly2019sd <- cellStats(x = LSTavgJuly2019r_ne, stat = "sd")
# LSTavgJuly2019Urbansd <- cellStats(x = LSTavgJuly2019r_neUrban, stat = "sd")
# LSTavgJuly2019Ruralsd <- cellStats(x = LSTavgJuly2019r_neRural, stat = "sd")
# LSTdayJuly2019sd <- cellStats(x = LSTdayJuly2019r_ne, stat = "sd")
# LSTdayJuly2019Urbansd <- cellStats(x = LSTdayJuly2019r_neUrban, stat = "sd")
# LSTdayJuly2019Ruralsd <- cellStats(x = LSTdayJuly2019r_neRural, stat = "sd")
# LSTnightJuly2019sd <- cellStats(x = LSTnightJuly2019r_ne, stat = "sd")
# LSTnightJuly2019Urbansd <- cellStats(x = LSTnightJuly2019r_neUrban, stat = "sd")
# LSTnightJuly2019Ruralsd <- cellStats(x = LSTnightJuly2019r_neRural, stat = "sd")
# 
# # Create rasters of cells exceeding threshold
# LSTavgJulyHI <- LSTavgJuly2019r_ne >
#   (LSTavgJuly2019mean + LSTavgJuly2019sd * 3)
# LSTdayJulyHI <- LSTdayJuly2019r_ne >
#   (LSTdayJuly2019mean + LSTdayJuly2019sd * 3)
# LSTnightJulyHI <- LSTnightJuly2019r_ne >
#   (LSTnightJuly2019mean + LSTnightJuly2019sd * 3)
# 
# # Extract raster values to block groups
# # check for empty geometries
# any(is.na(st_dimension(ne_blkgrp_sf)))
# # identiy empty geometries
# empty_geo <- st_is_empty(ne_blkgrp_sf)
# # filter out empty geometries
# ne_blkgrp_sf <- ne_blkgrp_sf[!empty_geo,]
# # clean up
# rm(empty_geo)
# 
# # Extract mean LST values within each block group to a df
# meanLSTavg_df <- ne_blkgrp_sf %>%
#   dplyr::select(GEOID) %>%
#   as_Spatial() %>%
#   extract(LSTavgJuly2019r_ne, .,
#           fun=mean, sp=TRUE, na.rm=TRUE, small=TRUE) %>%
#   as.data.frame() %>%
#   rename(meanAvgLST = layer)
# 
# meanLSTday_df <- ne_blkgrp_sf %>%
#   dplyr::select(GEOID) %>%
#   as_Spatial() %>%
#   extract(LSTdayJuly2019r_ne, .,
#           fun=mean, sp=TRUE, na.rm=TRUE, small=TRUE) %>%
#   as.data.frame() %>%
#   rename(meanDayLST = layer)
# 
# meanLSTnight_df <- ne_blkgrp_sf %>%
#   dplyr::select(GEOID) %>%
#   as_Spatial() %>%
#   extract(LSTnightJuly2019r_ne, .,
#           fun=mean, sp=TRUE, na.rm=TRUE, small=TRUE) %>%
#   as.data.frame() %>%
#   rename(meanNightLST = layer)
# 
# meanLSTDayNightDiff_df <- ne_blkgrp_sf %>%
#   dplyr::select(GEOID) %>%
#   as_Spatial() %>%
#   extract(LSTDayNighDiffJuly2019r_ne,.,
#           fun=mean, sp=TRUE, na.rm=TRUE, small=TRUE) %>%
#   as.data.frame() %>%
#   rename(meanDayNightDiffLST = layer)
# 
# # # Extract mean LST value for rural areas in order to compute urban heat island
# # meanLSTavgRural_df <- ne_blkgrp_rural_simple %>%
# #   as_Spatial() %>%
# #   extract(LSTavgJuly2019r_ne, .,
# #           fun=mean, sp=TRUE, na.rm=TRUE, small=TRUE) %>%
# #   as.data.frame() %>%
# #   rename(meanAvgLSTrural = layer)
# 
# # Join LST statistics to block groups
# ne_blkgrpLST_sf <- ne_blkgrp_sf %>%
#   left_join(., meanLSTavg_df, by = "GEOID") %>%
#   left_join(., meanLSTday_df, by = "GEOID") %>%
#   left_join(., meanLSTnight_df, by = "GEOID") %>%
#   left_join(., meanLSTDayNightDiff_df, by = "GEOID") %>%
#   mutate(UHI24avg = meanAvgLST - LSTavgJuly2019Ruralmean,
#          UHIDay = meanDayLST - LSTdayJuly2019Ruralmean,
#          UHINight = meanNightLST - LSTnightJuly2019Ruralmean)
# 
# save(ne_blkgrpLST_sf, ne_urban_sf, ne_rural_sf, file = "DATA/LST/ne_blkgrpLST_sf.rds")
# 
# load previously processed LST data
load("DATA/LST/ne_blkgrpLST_sf.rds")

# Calculate percentiles for pops of concern and temperature and create ranks for comparison
ne_blkgrpLST_sf <- ne_blkgrpLST_sf %>% 
  mutate(pctile24hrLST = percent_rank(meanAvgLST),
         rank24hrLST = case_when(
           pctile24hrLST < 0.2 ~ 1,
           pctile24hrLST >= 0.2 & pctile24hrLST < 0.4 ~ 2,
           pctile24hrLST >= 0.4 & pctile24hrLST < 0.6 ~ 3,
           pctile24hrLST >= 0.6 & pctile24hrLST < 0.8 ~ 4,
           pctile24hrLST >= 0.8 ~ 5),
         pctileDayLST = percent_rank(meanDayLST),
         rankDayLST = case_when(
           pctileDayLST < 0.2 ~ 1,
           pctileDayLST >= 0.2 & pctileDayLST < 0.4 ~ 2,
           pctileDayLST >= 0.4 & pctileDayLST < 0.6 ~ 3,
           pctileDayLST >= 0.6 & pctileDayLST < 0.8 ~ 4,
           pctileDayLST >= 0.8 ~ 5),
         pctileNightLST = percent_rank(meanNightLST),
         rankNightLST = case_when(
           pctileNightLST < 0.2 ~ 1,
           pctileNightLST >= 0.2 & pctileNightLST < 0.4 ~ 2,
           pctileNightLST >= 0.4 & pctileNightLST < 0.6 ~ 3,
           pctileNightLST >= 0.6 & pctileNightLST < 0.8 ~ 4,
           pctileNightLST >= 0.8 ~ 5),
         pctileMinority = percent_rank(minority_pctE),
         rankMinority = case_when(
           pctileMinority < 0.2 ~ 1,
           pctileMinority >= 0.2 & pctileMinority < 0.4 ~ 2,
           pctileMinority >= 0.4 & pctileMinority < 0.6 ~ 3,
           pctileMinority >= 0.6 & pctileMinority < 0.8 ~ 4,
           pctileMinority >= 0.8 ~ 5),
         LST24MinorScore = (rank24hrLST + rankMinority)/2,
         LSTDayMinorScore = (rankDayLST + rankMinority)/2,
         LSTNightMinorScore = (rankNightLST + rankMinority)/2,
         pctileEngLang = percent_rank(eng_limit_pctE),
         rankEngLang = case_when(
           pctileEngLang < 0.2 ~ 1,
           pctileEngLang >= 0.2 & pctileEngLang < 0.4 ~ 2,
           pctileEngLang >= 0.4 & pctileEngLang < 0.6 ~ 3,
           pctileEngLang >= 0.6 & pctileEngLang < 0.8 ~ 4,
           pctileEngLang >= 0.8 ~ 5),
         LST24EngLangScore = (rank24hrLST + rankEngLang)/2,
         LSTDayEngLangScore = (rankDayLST + rankEngLang)/2,
         LSTNightEngLangScore = (rankNightLST + rankEngLang)/2,
         pctileLowInc = percent_rank(pct2povE),
         rankLowInc = case_when(
           pctileLowInc < 0.2 ~ 1,
           pctileLowInc >= 0.2 & pctileLowInc < 0.4 ~ 2,
           pctileLowInc >= 0.4 & pctileLowInc < 0.6 ~ 3,
           pctileLowInc >= 0.6 & pctileLowInc < 0.8 ~ 4,
           pctileLowInc >= 0.8 ~ 5),
         LST24LowIncScore = (rank24hrLST + rankLowInc)/2,
         LSTDayLowIncScore = (rankDayLST + rankLowInc)/2,
         LSTNightLowIncScore = (rankNightLST + rankLowInc)/2,
         pctileNoHSDip = percent_rank(pct_lthsE),
         rankNoHSDip = case_when(
           pctileNoHSDip < 0.2 ~ 1,
           pctileNoHSDip >= 0.2 & pctileNoHSDip < 0.4 ~ 2,
           pctileNoHSDip >= 0.4 & pctileNoHSDip < 0.6 ~ 3,
           pctileNoHSDip >= 0.6 & pctileNoHSDip < 0.8 ~ 4,
           pctileNoHSDip >= 0.8 ~ 5),
         LST24NoHSDipScore = (rank24hrLST + rankNoHSDip)/2,
         LSTDayNoHSDipScore = (rankDayLST + rankNoHSDip)/2,
         LSTNightNoHSDipScore = (rankNightLST + rankNoHSDip)/2,
         pctileUnder5 = percent_rank(pct_under5E),
         rankUnder5 = case_when(
           pctileUnder5 < 0.2 ~ 1,
           pctileUnder5 >= 0.2 & pctileUnder5 < 0.4 ~ 2,
           pctileUnder5 >= 0.4 & pctileUnder5 < 0.6 ~ 3,
           pctileUnder5 >= 0.6 & pctileUnder5 < 0.8 ~ 4,
           pctileUnder5 >= 0.8 ~ 5),
         LST24Under5Score = (rank24hrLST + rankUnder5)/2,
         LSTDayUnder5Score = (rankDayLST + rankUnder5)/2,
         LSTNightUnder5Score = (rankNightLST + rankUnder5)/2,
         pctileOver64 = percent_rank(pct_over64E),
         rankOver64 = case_when(
           pctileOver64 < 0.2 ~ 1,
           pctileOver64 >= 0.2 & pctileOver64 < 0.4 ~ 2,
           pctileOver64 >= 0.4 & pctileOver64 < 0.6 ~ 3,
           pctileOver64 >= 0.6 & pctileOver64 < 0.8 ~ 4,
           pctileOver64 >= 0.8 ~ 5),
         LST24Over64Score = (rank24hrLST + rankOver64)/2,
         LSTDayOver64Score = (rankDayLST + rankOver64)/2,
         LSTNightOver64Score = (rankNightLST + rankOver64)/2
         )

# Calculate totals of populations of concern for use in computing percentages in later tables
statepopsdf <- ne_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  group_by(STATE) %>% 
  summarize(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))
```


## Day-Night Average Land Surface Temperatures

Day-Night Average Land Surface Temperatures (LST) represent a simple average of LST values collected during the day (11:48am - 2pm) and during the night (12am - 3:06am) over eight days, from July 28 to August 4, 2019. The day-night average for the entire region that week was `r round(mean(ne_blkgrpLST_sf$meanAvgLST,na.rm=TRUE),1)`°F. These day-night average temperatures varied significantly across New England, with highest temperatures apparent in the most urbanized areas of the region. This is particularly apparent along the I-95 corridor from Stamford, Connecticut to New Haven, and along the I-91 corridor from New Haven through Hartford to Springfield, Massachusetts. Higher temperatures are also apparent around greater Providence, Rhode Island, and the greater Boston region, reaching into southern New Hampshire (see Figure \@ref(fig:mapLSTavgNE)). 

```{r mapLSTavgNE, fig.align = "center", fig.cap="Map of July-Aug 2019 Day-Night Average Land Surface Temperatures (LST) across New England at Census Block Group level.", strip.white=TRUE}
# Map of average LST across New England
m <- tm_shape(ne_blkgrpLST_sf, unit = "mi") + 
  tm_fill("meanAvgLST", style = "quantile", palette = "Oranges",
          title = expression("Temperature " ( degree*F)),
          legend.hist = FALSE,
          colorNA = NULL,
          textNA = NULL,
          legend.format=list(list(digits=1)),
          legend.is.portrait = TRUE) +
  tm_shape(ne_states_sf_cb) + tm_borders(alpha = 0.4, lwd = 0.5) + 
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ne_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", xmod = 0.7, ymod = 0.2) +
  tm_scale_bar(breaks = c(0, 50, 100), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "Day-Night Average\nLand Surface\nTemperatures\nJuly-Aug 2019", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"),
            legend.hist.width = 0.9)

tmap_save(m, "DATA/LST/ne_24hrLST.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/ne_24hrLST.png")
```

```{r LST24outliers, include=FALSE}
# download municipal polygons
nh_towns_sf <- county_subdivisions("NH") %>%
  dplyr::select(NAME) %>%
  st_transform(., crs = st_crs(ne_blkgrpLST_sf))

me_towns_sf <- county_subdivisions("ME") %>%
  dplyr::select(NAME) %>%
  st_transform(., crs = st_crs(ne_blkgrpLST_sf))

# Identify municipalities with high heat outliers in New Hampshire and Maine
nh_LST24outliers <- ne_blkgrpLST_sf %>% 
  filter(STATE == "New Hampshire") %>% 
  filter(meanAvgLST == max(meanAvgLST, na.rm = TRUE)) %>% 
  st_centroid(.) %>% 
  st_intersection(nh_towns_sf,.) %>% 
  filter(row_number() == 1) %>% 
  pull(NAME)

me_LST24outliers <- ne_blkgrpLST_sf %>% 
  filter(STATE == "Maine") %>% 
  filter(meanAvgLST == max(meanAvgLST, na.rm = TRUE)) %>% 
  st_centroid(.) %>% 
  st_intersection(me_towns_sf,.) %>% 
  filter(row_number() == 1) %>% 
  pull(NAME)
```

Rhode Island, Massachusetts, and Connecticut are the warmest states in New England and exhibit significantly different median day-night average temperatures, although the highest temperatures are very similar (see Figure \@ref(fig:boxplotLSTavgNE)). By contrast, Vermont, Maine, and New Hampshire are significantly cooler overall, with Vermont exhibiting the coolest day-night average temperatures in New England. Although New Hampshire and Maine are cooler on average, some Block Groups show much higher temperatures, comparable to the highest temperatures in southern New England. This is the case in `r nh_LST24outliers`, New Hampshire and `r me_LST24outliers`, Maine. 

```{r boxplotLSTavgNE, fig.align = "center", fig.cap="Boxplot of July-Aug 2019 day-night average Land Surface Temperatures (LST) across New England at Census Block Group level."}
# boxplot of PM2.5 by state for 2016
ne_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  drop_na(meanAvgLST) %>% 
  ggplot(aes(x = reorder(STATE, meanAvgLST, FUN = "median", .desc = TRUE), 
             y = meanAvgLST, fill = STATE)) + 
  geom_boxplot(notch = TRUE) + 
  ggtitle("Day-Night Average Land Surface Temperatures July-Aug 2019") +
  theme_minimal() +
  scale_y_continuous(labels = function(x) paste0(x, "°")) +
  theme(legend.position = "none", axis.text=element_text(size=8)) + 
  xlab(NULL) + 
  ylab("Temperature (°F)") +
  coord_flip()
```


### Day-Night Average Urban Heat Island Effect

The urban heat island (UHI) effect is a phenomenon in which temperatures in urban areas tend to be higher than surrounding non-urban or rural areas. These elevated urban air temperatures are a consequence of the high density of buildings, roads, and other impervious infrastructure in urban areas that absorb heat and release it more slowly. The UHI effect is compounded by the relative absence of vegetation which would otherwise moderate temperatures through evaporative cooling.[^Zhou2019] Higher summer temperatures due to UHI increase the risk of heat-related illnesses, result in increased energy use for air conditioning, and exacerbate air pollution, especially ground-level ozone. As the climate warms, UHI are likely to exacerbate both temperature and air pollution risks. 

The UHI effect analyzed here is defined as the difference in temperature between urbanized areas and rural areas in New England.[^urban] Specifically, the average temperature of the rural areas of New England is subtracted from the average temperatures of each Census Block Group. The resulting value shows how much warmer or cooler each Census Block Group is when compared to the rural background reference. Positive values indicate that a Census Block Group is warmer than the rural background average and may be experiencing the UHI effect. The larger the value, the stronger the effect. 

The UHI effect for day-night average temperatures in New England are most pronounced in the inner core of urbanized areas (see Figure \@ref(fig:mapUHIavgNE)), particularly around Hartford, Connecticut, Providence, Rhode Island, and Boston, Massachusetts. Temperatures in these areas are up to `r round(max(ne_blkgrpLST_sf$UHI24avg, na.rm = TRUE),1)`° warmer than the rural background average. By contrast, most of Vermont, northern New Hampshire, and northwestern Maine are cooler than rural background average. 

```{r mapUHIavgNE, fig.align = "center", fig.cap="Map of July-Aug 2019 day-night average Urban Heat Island (UHI) effect across New England at Census Block Group level. UHI is the difference between Block Group average temperature and average temperature of rural areas.", strip.white=TRUE}
# Map of average LST UHI across New England
m <- tm_shape(ne_blkgrpLST_sf, unit = "mi") + 
  tm_fill("UHI24avg", palette = "-RdYlGn", midpoint = 0,
          legend.reverse = TRUE,
          title = expression("Temperature " ( degree*F)),
          legend.hist = FALSE,
          colorNA = NULL,
          textNA = NULL,
          legend.format=list(list(digits=1)),
          legend.is.portrait = TRUE) +
  tm_shape(ne_states_sf_cb) + tm_borders(alpha = 0.4, lwd = 0.5) + 
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ne_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", xmod = 0.7, ymod = 0.2) +
  tm_scale_bar(breaks = c(0, 50, 100), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "Day-Night Average\nUrban Heat Island\nJuly-Aug 2019", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"),
            legend.hist.width = 0.9)

tmap_save(m, "DATA/LST/ne_24hrUHI.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/ne_24hrUHI.png")
```

Rhode Island, Massachusetts, and Connecticut exhibit median UHI effects that are significantly different from one another, although the highest values are similar with respect to the UHI effect (see Figure \@ref(fig:boxplotUHIavgNE)). In these states, median UHI effects range from `r ne_blkgrpLST_sf %>% as.data.frame() %>% filter(STATE == "Connecticut") %>% summarize(round(median(UHI24avg, na.rm = TRUE),1)) %>% pull()`° to `r ne_blkgrpLST_sf %>% as.data.frame() %>% filter(STATE == "Rhode Island") %>% summarize(round(median(UHI24avg, na.rm = TRUE),1)) %>% pull()`° above the rural background. Northern New England states show significantly lower UHI effect values. Vermont's median value is only slightly above 0°, indicating little difference from the rural background temperature. 

```{r boxplotUHIavgNE, fig.align = "center", fig.cap="Boxplot of July-Aug 2019 day-night average Urban Heat Island (UHI) across New England at Census Block Group level. UHI is the difference between Block Group average temperature and average temperature of rural areas."}
# boxplot of PM2.5 by state for 2016
ne_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  drop_na(UHI24avg) %>% 
  ggplot(aes(x = reorder(STATE, UHI24avg, FUN = "median", .desc = TRUE), 
             y = UHI24avg, fill = STATE)) + 
  geom_boxplot(notch = TRUE) + 
  ggtitle("Day-Night Average UHI July-Aug 2019") +
  theme_minimal() +
  scale_y_continuous(labels = function(x) paste0(x, "°")) +
  theme(legend.position = "none", axis.text=element_text(size=8)) + 
  xlab(NULL) + 
  ylab("Temperature (°F)") +
  coord_flip()
```



### Day-Night Average Land Surface Temperatures and Populations of Concern

```{r LST24hrTempPop, include=FALSE}
# Pop Weighted avg of 24hr LST exposure for all Groups in New England relative to NE average
LST24hrPop <- ne_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  # filter(STATE == "Massachusetts") %>%
  dplyr::select(totalpopE,
                nhwhitepopE,
                minorityE,
                nhblackpopE,
                nhamerindpopE,
                nhasianpopE,
                nhnativhpopE,
                nhotherpopE,
                nh2morepopE,
                hisppopE,
                povknownE,
                num2povE, 
                eng_hhE,
                eng_limitE,
                age25upE,
                lthsE, 
                allAgesE, 
                under5E, 
                over64E, 
                meanAvgLST) %>% 
  gather(key = Group, value = Pop, totalpopE:over64E) %>% 
  group_by(Group) %>% 
  summarize(LST24hrwMean = weighted.mean(x = meanAvgLST, w = Pop, na.rm = TRUE),
            LST24hrMean = mean(meanAvgLST, na.rm = TRUE)) %>% 
  filter(Group %in% c("minorityE","eng_limitE","num2povE","lthsE","under5E","over64E")) %>% 
  mutate(Group = case_when(
    Group == "minorityE" ~ "Minority",
    Group == "eng_limitE" ~ "Limited English HH",
    Group == "num2povE" ~ "Low Income",
    Group == "lthsE" ~ "No HS Dip",
    Group == "under5E" ~ "Under 5",
    Group == "over64E" ~ "Over 64"))


# Pop Weighted avg of LST24hr for all Groups in New England relative to NE average
LST24hrPopTab <- ne_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  # filter(STATE == "Massachusetts") %>%
  dplyr::select(totalpopE,
                nhwhitepopE,
                minorityE,
                nhblackpopE,
                nhamerindpopE,
                nhasianpopE,
                nhnativhpopE,
                nhotherpopE,
                nh2morepopE,
                hisppopE,
                povknownE,
                num2povE, 
                eng_hhE,
                eng_limitE,
                age25upE,
                lthsE, 
                allAgesE, 
                under5E, 
                over64E, 
                meanAvgLST) %>% 
  gather(key = Group, value = Pop, totalpopE:over64E) %>% 
  group_by(Group) %>% 
  summarize(LST24hrwMean = weighted.mean(x = meanAvgLST, w = Pop, na.rm = TRUE),
            LST24hrMean = mean(meanAvgLST, na.rm = TRUE)) %>% 
  spread(key = Group, value = LST24hrwMean) %>% 
  transmute(Minority = (minorityE/LST24hrMean - 1)*100,
         #Minority_NHW = (minorityE/nhwhitepopE - 1)*100,
         `Limited English HH` = (eng_limitE/LST24hrMean - 1)*100,
         `Low Income` = (num2povE/LST24hrMean - 1)*100,
         `No HS Dip` = (lthsE/LST24hrMean - 1)*100,
         `Under 5` = (under5E/LST24hrMean - 1)*100,
         `Over 64` = (over64E/LST24hrMean - 1)*100) %>%
  gather(key = Group, value = Pct) %>% 
  left_join(x = LST24hrPop, y = ., by = "Group") %>% 
  mutate(Difference = LST24hrwMean - LST24hrMean) %>% 
  transmute(Group = Group,
            LST24hrwMean = paste0(round(LST24hrwMean,2),"°"),
            Difference = if_else(Difference > 0,
                                 paste0("+",round(Difference,2),"°"),
                                 paste0(round(Difference,2),"°")),
            Pct = if_else(Pct > 0, paste0("+",round(Pct,2),"%"),
                       paste0(round(Pct,2),"%")))
```

In addition to variations in the general geography of Land Surface Temperatures (LST), temperature exposure also varies demographically. Figure \@ref(fig:plotLST24hrTempPopAvgNE) and Table \@ref(tab:tabLST24hrPopAvgNE) show population-weighted exposures for populations of concern relative to the day-night average LST for the region. For example, limited English speaking households in New England were exposed to a population-weighted day-night average LST of `r round(max(LST24hrPop$LST24hrwMean),1)`°, approximately `r round(max(LST24hrPop$LST24hrwMean) - max(LST24hrPop$LST24hrMean),1)`° (`r LST24hrPopTab[1,4]`) above the regional day-night average LST. Similarly, minorities were exposed to a population-weighted day-night average LST of `r LST24hrPop %>% filter(Group == "Minority") %>% summarize(round(max(LST24hrwMean),1)) %>% pull()`°, over `r round(LST24hrPop %>% filter(Group == "Minority") %>% summarize(max(LST24hrwMean)) %>% pull() - max(LST24hrPop$LST24hrMean),1)`° (`r LST24hrPopTab[3,4]`) above the regional day-night average LST. By contrast, persons over age 64 were, on average, exposed to population-weighted day-night average LSTs of `r LST24hrPopTab[5,3]` below the regional average.

```{r plotLST24hrTempPopAvgNE, echo=FALSE, message=FALSE, warning=FALSE, fig.align = "center", fig.cap="Population-weighted average exposures to day-night average temperatures for populations of concern in New England relative to the regional average."}
LST24hrPop %>% 
  ggplot(aes(x = reorder(Group,LST24hrwMean), 
             y = LST24hrwMean)) +
  geom_segment(aes(x = reorder(Group,LST24hrwMean), 
                   xend = reorder(Group,LST24hrwMean),
                   y = LST24hrMean, yend = LST24hrwMean), 
               color = "firebrick1") +
  geom_point(color = "firebrick4", size = 4, alpha = 0.8) +
  coord_flip() + 
  labs(x = "", y = "", title = "Population-Weighted Temperature Exposure to\nDay-Night Average LST (°F)") + 
  theme_light() +
  theme(panel.grid.major.y = element_blank(),
        panel.border = element_blank(),
        axis.ticks.y = element_blank()) +
  geom_text(aes(x = Group, y = LST24hrwMean + 0.2 * sign(LST24hrwMean), 
                label = paste0(round(LST24hrwMean,1),"°")), 
            hjust = 1.6, vjust = -.8, size = 3,
            color=rgb(100,100,100, maxColorValue=255)) +
  scale_y_continuous(labels = function(x) paste0(x, "°")) +
  geom_hline(yintercept = LST24hrPop$LST24hrMean, linetype = "dashed") +
  geom_text(aes(x = "Low Income", y = 80.5, label = "Above regional\naverage"),
            color = "gray48") +
  geom_text(aes(x = "Low Income", y = 77, label = "Below regional\naverage"),
            color = "gray48") +
  # geom_segment(aes(x = "Under 5", xend = "Under 5", y = 79.5, yend = 81), 
  #              arrow = arrow(length = unit(0.3,"cm"))) +
  expand_limits(y = c(76,83))
```


```{r tabLST24hrPopAvgNE, fig.align = "center", fig.cap="Population-weighted average exposures to day-night average LST for populations of concern in New England relative to the regional average."}
LST24hrPopTab %>% 
  arrange(Group) %>% 
  kableExtra::kable(longtable = T, booktabs = T, digits = 2, 
                    col.names = c("Group","Avg 24hr LST (°F)",
                                  "Difference from Regional Avg* (°F)",
                                  "Pct Above/Below Regional Avg*"), 
                    caption = "Population-Weighted Temperature Exposure - Day-Night Average LST", align = "r") %>% 
  kableExtra::column_spec(2:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header")) %>% 
  kableExtra::footnote(., symbol = paste0("Regional average day-night LST is ", round(mean(ne_blkgrpLST_sf$meanAvgLST,na.rm=TRUE),2), "°"))
```

```{r UHI24hrTempPop, include=FALSE}
# Pop Weighted avg of 24hr LST exposure for all Groups in New England relative to NE average
UHI24hrPop <- ne_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  # filter(STATE == "Massachusetts") %>%
  dplyr::select(totalpopE,
                nhwhitepopE,
                minorityE,
                nhblackpopE,
                nhamerindpopE,
                nhasianpopE,
                nhnativhpopE,
                nhotherpopE,
                nh2morepopE,
                hisppopE,
                povknownE,
                num2povE, 
                eng_hhE,
                eng_limitE,
                age25upE,
                lthsE, 
                allAgesE, 
                under5E, 
                over64E, 
                UHI24avg) %>% 
  gather(key = Group, value = Pop, totalpopE:over64E) %>% 
  group_by(Group) %>% 
  summarize(UHI24hrwMean = weighted.mean(x = UHI24avg, w = Pop, na.rm = TRUE),
            UHI24hrMean = mean(UHI24avg, na.rm = TRUE)) %>% 
  filter(Group %in% c("minorityE","eng_limitE","num2povE","lthsE","under5E","over64E")) %>% 
  mutate(Group = case_when(
    Group == "minorityE" ~ "Minority",
    Group == "eng_limitE" ~ "Limited English HH",
    Group == "num2povE" ~ "Low Income",
    Group == "lthsE" ~ "No HS Dip",
    Group == "under5E" ~ "Under 5",
    Group == "over64E" ~ "Over 64"))

# Pop Weighted avg of PM2.5 for all Groups in New England relative to NE average
UHI24hrPopTab <- ne_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  # filter(STATE == "Massachusetts") %>%
  dplyr::select(totalpopE,
                nhwhitepopE,
                minorityE,
                nhblackpopE,
                nhamerindpopE,
                nhasianpopE,
                nhnativhpopE,
                nhotherpopE,
                nh2morepopE,
                hisppopE,
                povknownE,
                num2povE, 
                eng_hhE,
                eng_limitE,
                age25upE,
                lthsE, 
                allAgesE, 
                under5E, 
                over64E, 
                UHI24avg) %>% 
  gather(key = Group, value = Pop, totalpopE:over64E) %>% 
  group_by(Group) %>% 
  summarize(UHI24hrwMean = weighted.mean(x = UHI24avg, w = Pop, na.rm = TRUE),
            UHI24hrMean = mean(UHI24avg, na.rm = TRUE)) %>% 
  spread(key = Group, value = UHI24hrwMean) %>% 
  transmute(Minority = (minorityE/UHI24hrMean - 1)*100,
         #Minority_NHW = (minorityE/nhwhitepopE - 1)*100,
         `Limited English HH` = (eng_limitE/UHI24hrMean - 1)*100,
         `Low Income` = (num2povE/UHI24hrMean - 1)*100,
         `No HS Dip` = (lthsE/UHI24hrMean - 1)*100,
         `Under 5` = (under5E/UHI24hrMean - 1)*100,
         `Over 64` = (over64E/UHI24hrMean - 1)*100) %>%
  gather(key = Group, value = Pct) %>% 
  left_join(x = UHI24hrPop, y = ., by = "Group") %>% 
  mutate(Difference = UHI24hrwMean - UHI24hrMean) %>% 
  transmute(Group = Group,
            UHI24hrwMean = paste0(round(UHI24hrwMean,2),"°"),
            Difference = if_else(Difference > 0,
                                 paste0("+",round(Difference,2),"°"),
                                 paste0(round(Difference,2),"°")),
            Pct = if_else(Pct > 0, paste0("+",round(Pct,2),"%"),
                       paste0(round(Pct,2),"%")))
```

Figure \@ref(fig:plotUHI24hrTempPopAvgNE) and Table \@ref(tab:tabUHI24hrPopAvgNE) show population-weighted exposures for populations of concern relative to the day-night average Urban Heat Island (UHI) effect for New England. All populations of concern are exposed to the UHI effect on average. Indeed, the average UHI effect across New England is approximately `r round(mean(ne_blkgrpLST_sf$UHI24avg,na.rm=TRUE),1)`°. However, some populations experience significantly higher UHI effects. For example, limited English speaking households in New England were exposed to a population-weighted day-night average UHI effect of `r round(max(UHI24hrPop$UHI24hrwMean),1)`° above the rural average. This value is more than `r round(max(UHI24hrPop$UHI24hrwMean) - max(UHI24hrPop$UHI24hrMean),1)`° (`r UHI24hrPopTab[1,4]`) above the regional day-night average UHI. Similarly, minorities were exposed to a population-weighted day-night average UHI effect of `r UHI24hrPop %>% filter(Group == "Minority") %>% summarize(round(max(UHI24hrwMean),1)) %>% pull()`°, approximately `r round(UHI24hrPop %>% filter(Group == "Minority") %>% summarize(max(UHI24hrwMean)) %>% pull() - max(UHI24hrPop$UHI24hrMean),1)`° (`r UHI24hrPopTab[3,4]`) above the regional day-night average UHI. By contrast, persons over age 64 were, on average, were exposed to a population-weighted day-night average UHI of `r UHI24hrPopTab[5,3]` (`r UHI24hrPopTab[5,4]`) below the regional average.

```{r plotUHI24hrTempPopAvgNE, echo=FALSE, message=FALSE, warning=FALSE, fig.align = "center", fig.cap="Population-weighted average exposures to day-night average Urban Heat Island effect for populations of concern in New England."}
UHI24hrPop %>% 
  ggplot(aes(x = reorder(Group,UHI24hrwMean), 
             y = UHI24hrwMean)) +
  geom_segment(aes(x = reorder(Group,UHI24hrwMean), 
                   xend = reorder(Group,UHI24hrwMean),
                   y = 0, yend = UHI24hrwMean), 
               color = "firebrick1") +
  geom_point(color = "firebrick4", size = 4, alpha = 0.8) +
  coord_flip() + 
  labs(x = "", y = "", title = "Population-Weighted Temperature Exposure to\nDay-Night Average UHI (°F)") + 
  theme_light() +
  theme(panel.grid.major.y = element_blank(),
        panel.border = element_blank(),
        axis.ticks.y = element_blank()) +
  geom_text(aes(x = Group, y = UHI24hrwMean + 0.2 * sign(UHI24hrwMean), 
                label = paste0(round(UHI24hrwMean,1),"°")), 
            hjust = 1.6, vjust = -.8, size = 3,
            color=rgb(100,100,100, maxColorValue=255)) +
  scale_y_continuous(labels = function(x) paste0(x, "°")) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  # geom_text(aes(x = "Under 5", y = 12, label = "Above\nrural average"),
  #           color = "gray48") +
  # geom_segment(aes(x = "Over 64", xend = "Over 64", y = 9, yend = 12), 
  #              arrow = arrow(length = unit(0.3,"cm"))) +
  expand_limits(y = c(-2,15))
```


```{r tabUHI24hrPopAvgNE, fig.align = "center", fig.cap="Population-weighted average exposures to day-night average Urban Heat Island effect for populations of concern in New England relative to the rural average."}
UHI24hrPopTab %>% 
  arrange(Group) %>% 
  kableExtra::kable(longtable = T, booktabs = T, digits = 2, 
                    col.names = c("Group","Day-Night Avg UHI (°F)",
                                  "Difference from Regional UHI Avg* (°F)",
                                  "Pct Above/Below Regional UHI Avg*"), 
                    caption = "Population-Weighted Temperature Exposure - Day-Night Avg UHI", align = "r") %>% 
  kableExtra::column_spec(2:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header")) %>% 
  kableExtra::footnote(., symbol = paste0("Regional average day-night UHI effect is ", round(mean(ne_blkgrpLST_sf$UHI24avg,na.rm=TRUE),2), "°"))
```

Because the risk of heat-related illness is a combination of both exposure and vulnerability, it is important to consider where these factors come together. Figures \@ref(fig:mapLST24hrHighRiskLEH) and \@ref(fig:mapLST24hrHighRiskMIN) highlight Census Block Groups across New England which exhibited day-night average LSTs in the top quintile (i.e., 80th percentile) and also contained high percentages of limited English speaking households or minorities in the top quintile for the region.[^LST24quintile] These high risk locations were present across four New England states in the late summer of 2019. Tables \@ref(tab:tabLST24hrHighRiskLEH) and \@ref(tab:tabLST24hrHighRiskMIN) show breakdowns of how many households these represent for each state, and the day-night average LST values for those same Census Block Groups. Note that Massachusetts contains the largest number of these at-risk households, over 67,000, as well as at-risk minority persons, over 650,000, while Rhode Island shows the highest day-night average LSTs for these at-risk housholds, on average over 86.7°, and over 87° for at-risk minorities. Maps and tables for other populations of concern can be found in Figures \@ref(fig:mapLST24hrHighRiskLI) to \@ref(fig:mapLST24hrHighRiskU5) and Tables \@ref(tab:tabLST24hrHighRiskLI) to \@ref(tab:tabLST24hrHighRiskU5) in Appendix B.

```{r mapLST24hrHighRiskLEH, fig.align = "center", fig.cap="Map of Block Groups in the highest quintile (80th percentile) of day-night average LST and 80th percentile of percentage of Limited English Speaking Households by Census Block Group, July-Aug 2019."}
LST24EngLang5 <- ne_blkgrpLST_sf %>%
  filter(LST24EngLangScore == 5)

bbox_new <- st_bbox(LST24EngLang5) # current bounding box

xrange <- bbox_new$xmax - bbox_new$xmin # range of x values
yrange <- bbox_new$ymax - bbox_new$ymin # range of y values

# bbox_new[1] <- bbox_new[1] - (0.5 * xrange) # xmin - left
bbox_new[3] <- bbox_new[3] + (0.2 * xrange) # xmax - right
# bbox_new[2] <- bbox_new[2] - (0.5 * yrange) # ymin - bottom
bbox_new[4] <- bbox_new[4] + (0.2 * yrange) # ymax - top

bbox_new <- bbox_new %>%  # take the bounding box ...
  st_as_sfc() # ... and make it a sf polygon

m <- tm_shape(LST24EngLang5, unit = "mi" , bbox = bbox_new) + tm_fill(col = "red", title = "80th Percentile", legend.is.portrait = TRUE) +
  tm_shape(ne_states_sf_cb) + tm_polygons(alpha = 0.4, lwd = 0.2) +
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ne_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", xmod = 0.7, ymod = 0.2) +
  tm_scale_bar(breaks = c(0, 50, 100), text.size = 0.5, 
               position = c(0.4,0.005)) +
  tm_add_legend(type = "fill", 
                labels = "80th Percentile Block Groups",
                col = "red",
                border.col = "white", border.alpha = 0) +
  tm_layout(title = "Highest Exposure\nand Vulnerability\nto Day-Night Avg\nTemps for Limited\nEnglish Speaking\nHouseholds\nJuly-Aug 2019", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"),
            legend.hist.width = 0.9)

tmap_save(m, "DATA/LST/ne_LST24EngLang5.png",
          height = 5, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/ne_LST24EngLang5.png")
```


```{r tabLST24hrHighRiskLEH, fig.align = "center", fig.cap="Limited English Housheolds in Census Block Groups at or above 80th percentile of day-night average Land Surface Temperature."}
ne_blkgrpLST_sf %>%
  as.data.frame() %>%
  left_join(., statepopsdf,by="STATE") %>%
  filter(LST24EngLangScore == 5) %>%
  group_by(STATE) %>%
  summarize(`Limited English Households` = sum(eng_limitE),
            `Pct of Limited English Households` = round(sum(eng_limitE)/max(TotalLEH)*100,1),
            Min = min(meanAvgLST,na.rm=TRUE),
            Mean = mean(meanAvgLST,na.rm=TRUE),
            Max = max(meanAvgLST,na.rm=TRUE)) %>% 
  mutate(`Pct of Limited English Households` = paste0(`Pct of Limited English Households`,"%"),
         Min = paste0(round(Min,2),"°"),
         Mean = paste0(round(Mean,2),"°"),
         Max = paste0(round(Max,2),"°")) %>% 
  kableExtra::kable(longtable = F, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Limited English Speaking Households in Highest Quintile of Day-Night Avg LST", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header")) %>% 
  kableExtra::column_spec(2:3, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 3, "Day-Night LST (°F)" = 3))
```

```{r mapLST24hrHighRiskMIN, fig.align = "center", fig.cap="Map of Block Groups in the highest quintile (80th percentile) of day-night average LST and 80th percentile of percentage of Minorities by Census Block Group, July-Aug 2019."}
LST24Minority5 <- ne_blkgrpLST_sf %>%
  filter(LST24MinorScore == 5)

bbox_new <- st_bbox(LST24EngLang5) # current bounding box

xrange <- bbox_new$xmax - bbox_new$xmin # range of x values
yrange <- bbox_new$ymax - bbox_new$ymin # range of y values

# bbox_new[1] <- bbox_new[1] - (0.5 * xrange) # xmin - left
bbox_new[3] <- bbox_new[3] + (0.2 * xrange) # xmax - right
# bbox_new[2] <- bbox_new[2] - (0.5 * yrange) # ymin - bottom
bbox_new[4] <- bbox_new[4] + (0.2 * yrange) # ymax - top

bbox_new <- bbox_new %>%  # take the bounding box ...
  st_as_sfc() # ... and make it a sf polygon

m <- tm_shape(LST24Minority5, unit = "mi", bbox = bbox_new) + tm_fill(col = "red", title = "80th Percentile", legend.is.portrait = TRUE) +
  tm_shape(ne_states_sf_cb) + tm_polygons(alpha = 0.4, lwd = 0.2) +
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ne_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", xmod = 0.7, ymod = 0.2) +
  tm_scale_bar(breaks = c(0, 50, 100), text.size = 0.5, 
               position = c(0.4,0.005)) +
  tm_add_legend(type = "fill", 
                labels = "80th Percentile Block Groups",
                col = "red",
                border.col = "white", border.alpha = 0) +
  tm_layout(title = "Highest Exposure\nand Vulnerability\nto Day-Night Avg\nTemps for\nMinorities\nJuly-Aug 2019", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"),
            legend.hist.width = 0.9)

tmap_save(m, "DATA/LST/ne_LST24Minority5.png",
          height = 5, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/ne_LST24Minority5.png")
```

\pagebreak
```{r tabLST24hrHighRiskMIN, fig.align = "center", fig.cap="Minorities in Census Block Groups at or above 80th percentile of day-night average Land Surface Temperature."}
ne_blkgrpLST_sf %>%
  as.data.frame() %>%
  left_join(., statepopsdf,by="STATE") %>%
  filter(LST24MinorScore == 5) %>%
  group_by(STATE) %>%
  summarize(Minorities = sum(minorityE),
            `Pct of Minorities` = round(sum(minorityE)/max(TotalMin)*100,1),
            Min = min(meanAvgLST,na.rm=TRUE),
            Mean = mean(meanAvgLST,na.rm=TRUE),
            Max = max(meanAvgLST,na.rm=TRUE)) %>% 
  mutate(`Pct of Minorities` = paste0(`Pct of Minorities`,"%"),
         Min = paste0(round(Min,2),"°"),
         Mean = paste0(round(Mean,2),"°"),
         Max = paste0(round(Max,2),"°")) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Minorities in Highest Quintile of Day-Night Avg LST", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header")) %>% 
  kableExtra::column_spec(2:3, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 3, "Day-Night LST (°F)" = 3))
```

\pagebreak
## Daytime Average Land Surface Temperatures

Daytime Average Land Surface Temperatures (LST) represent LST values collected during the day (11:48am - 2pm) over eight days, from July 28 to August 4, 2019. The average daytime LST for the region was `r round(mean(ne_blkgrpLST_sf$meanDayLST,na.rm=TRUE),1)`°F. Daytime Land Surface Temperatures (LST) varied significantly across New England, with highest temperatures apparent in the most urbanized areas of the region, particularly along the I-95 corridor from Stamford, Connecticut to New Haven, and along the I-91 corridor from New Haven through Hartford to Springfield, Massachusetts. This is also apparent around greater Providence, Rhode Island, and the greater Boston region out to I-495, and reaching into southern New Hampshire. Note that pockets of high daytime LST are also present in pockets across Maine (see Figure \@ref(fig:mapLSTdayNE) below). Across New England, approximately `r ne_blkgrpLST_sf %>% as.data.frame() %>% filter(pctileDayLST >= 0.8) %>% group_by() %>% summarize(sum(totalpopE, na.rm = TRUE)/10^6) %>% pull() %>% round(.,2)` million people live in Block Groups that exceeded the 80th percentile (temperatures over `r round(quantile(ne_blkgrpLST_sf$meanDayLST,0.8,na.rm=TRUE),1)`°) for daytime average temperatures.  

```{r mapLSTdayNE, fig.align = "center", fig.cap="Map of July-August 2019 average daytime Land Surface Temperatures (LST) across New England at Census Block Group level.", strip.white=TRUE}
# Map of average LST across New England
m <- tm_shape(ne_blkgrpLST_sf, unit = "mi") + 
  tm_fill("meanDayLST", style = "quantile", palette = "Reds",
          title = expression("Temperature " ( degree*F)),
          legend.hist = FALSE,
          colorNA = NULL,
          textNA = NULL,
          legend.format=list(list(digits=1)),
          legend.is.portrait = TRUE) +
  tm_shape(ne_states_sf_cb) + tm_borders(alpha = 0.4, lwd = 0.5) + 
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ne_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", xmod = 0.7, ymod = 0.2) +
  tm_scale_bar(breaks = c(0, 50, 100), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "Average Day\nLand Surface\nTemperatures\nJuly-Aug 2019", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"),
            legend.hist.width = 0.9)

tmap_save(m, "DATA/LST/ne_DayLST.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/ne_DayLST.png")
```

```{r LSTDayoutliers, include=FALSE}
# # download municipal polygons
# nh_towns_sf <- county_subdivisions("NH") %>%
#   dplyr::select(NAME) %>%
#   st_transform(., crs = st_crs(ne_blkgrpLST_sf))
# 
# me_towns_sf <- county_subdivisions("ME") %>%
#   dplyr::select(NAME) %>%
#   st_transform(., crs = st_crs(ne_blkgrpLST_sf))

# Identify municipalities with high heat outliers in New Hampshire and Maine
nh_LSTDayoutliers <- ne_blkgrpLST_sf %>% 
  filter(STATE == "New Hampshire") %>% 
  filter(meanDayLST == max(meanDayLST, na.rm = TRUE)) %>% 
  st_centroid(.) %>% 
  st_intersection(nh_towns_sf,.) %>% 
  filter(row_number() == 1) %>% 
  pull(NAME)

me_LSTDayoutliers <- ne_blkgrpLST_sf %>% 
  filter(STATE == "Maine") %>% 
  filter(meanDayLST == max(meanDayLST, na.rm = TRUE)) %>% 
  st_centroid(.) %>% 
  st_intersection(me_towns_sf,.) %>% 
  filter(row_number() == 1) %>% 
  pull(NAME)
```

Rhode Island, Massachusetts, and Connecticut are the warmest states in New England, and their average daytime temperatures are statistically significantly different (see Figure \@ref(fig:boxplotLSTdayNE)). Vermont, Maine, and New Hampshire are significantly cooler overall, with Vermont exhibiting the coolest daytime  temperatures in New England. Although New Hampshire and Maine are cooler on average, some Block Groups show much higher temperatures, comparable to southern New England. This is the case in `r nh_LSTDayoutliers`, New Hampshire and `r me_LSTDayoutliers`, Maine.

```{r boxplotLSTdayNE, fig.align = "center", fig.cap="Boxplot of July-Aug 2019 daytime average Land Surface Temperatures (LST) across New England at Census Block Group level."}
# # Test for statistically significant differences
# # First, set up a one-way ANOVA
# amod <- ne_blkgrpLST_sf %>% 
#   mutate(STATE = as.factor(STATE)) %>% 
#   aov(meanDayLST ~ STATE, data = .)
# library(multcomp)
# # Next, set up all-pair comparisons for factor `STATE`. See https://stats.stackexchange.com/questions/341463/test-for-significant-differences-between-groups-in-r for example
# g <- glht(amod, linfct = mcp(STATE = "Tukey"))
# summary(g)

# boxplot of PM2.5 by state for 2016
ne_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  drop_na(meanDayLST) %>% 
  ggplot(aes(x = reorder(STATE, meanDayLST, FUN = "median", .desc = TRUE), 
             y = meanDayLST, fill = STATE)) + 
  geom_boxplot(notch = TRUE) + 
  ggtitle("Average Day Land Surface Temperatures July-Aug 2019") +
  theme_minimal() +
  scale_y_continuous(labels = function(x) paste0(x, "°")) +
  theme(legend.position = "none", axis.text=element_text(size=8)) + 
  xlab(NULL) + 
  ylab("Temperature (°F)") +
  coord_flip()
```


### Daytime Urban Heat Island Effect

The UHI effect for daytime average temperatures in New England are most pronounced in the more densely urbanized areas (see Figure \@ref(fig:mapUHIdayNE)), particularly around Bridgeport and Hartford, Connecticut, around Boston, Massachusetts, and around Providence, Rhode Island. Pockets of high UHI are also evident around cities in southern New Hampshire. Temperatures in these areas are up to `r round(max(ne_blkgrpLST_sf$UHIDay,na.rm = TRUE),1)`° warmer than the rural background average. By contrast, at least half of Vermont, northern New Hampshire, and northwestern Maine are cooler than rural background average. 

```{r mapUHIdayNE, fig.align = "center", fig.cap="Map of July-Aug 2019 daytime Urban Heat Island (UHI) effect across New England at Census Block Group level. UHI is the difference between Block Group average temperature and average temperature of rural areas.", strip.white=TRUE}
# Map of average LST UHI across New England
m <- tm_shape(ne_blkgrpLST_sf, unit = "mi") + 
  tm_fill("UHIDay", palette = "-RdYlGn", midpoint = 0,
          legend.reverse = TRUE,
          title = expression("Temperature " ( degree*F)),
          legend.hist = FALSE,
          colorNA = NULL,
          textNA = NULL,
          legend.format=list(list(digits=1)),
          legend.is.portrait = TRUE) +
  tm_shape(ne_states_sf_cb) + tm_borders(alpha = 0.4, lwd = 0.5) + 
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  # tm_shape(ne_urban_sf) + tm_fill(alpha = 0.5) +
  tm_shape(ne_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", xmod = 0.7, ymod = 0.2) +
  tm_scale_bar(breaks = c(0, 50, 100), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "Daytime\nUrban Heat Island\nJuly-Aug 2019", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"),
            legend.hist.width = 0.9)

tmap_save(m, "DATA/LST/ne_DayUHI.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/ne_DayUHI.png")
```

Rhode Island, Massachusetts, and Connecticut exhibit the highest UHI effects, and each is statistically significantly different from the others (see Figure \@ref(fig:boxplotUHIDayNE)). In these states, median UHI effect ranged from `r ne_blkgrpLST_sf %>% as.data.frame() %>% filter(STATE == "Connecticut") %>% summarize(round(median(UHIDay,na.rm = TRUE),1)) %>% pull()`° in Connecticut to `r ne_blkgrpLST_sf %>% as.data.frame() %>% filter(STATE == "Rhode Island") %>% summarize(round(median(UHIDay,na.rm = TRUE),1)) %>% pull()`° in Rhode Island above the rural background. Northern New England states show significantly lower UHI effect values. Vermont's median value is only slightly above 0°, indicating minor difference from the rural background temperature. 

```{r boxplotUHIDayNE, fig.align = "center", fig.cap="Boxplot of July-Aug 2019 daytime Urban Heat Island (UHI) across New England at Census Block Group level. UHI is the difference between Block Group average temperature and average temperature of rural areas."}
# Test for statistically significant differences
# First, set up a one-way ANOVA
# amod <- ne_blkgrpLST_sf %>% 
#   mutate(STATE = as.factor(STATE)) %>% 
#   aov(UHIDay ~ STATE, data = .)
# library(multcomp)
# # Next, set up all-pair comparisons for factor `STATE`. See https://stats.stackexchange.com/questions/341463/test-for-significant-differences-between-groups-in-r for example
# g <- glht(amod, linfct = mcp(STATE = "Tukey"))
# summary(g)

# boxplot of PM2.5 by state for 2016
ne_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  drop_na(UHIDay) %>% 
  ggplot(aes(x = reorder(STATE, UHIDay, FUN = "median", .desc = TRUE), 
             y = UHIDay, fill = STATE)) + 
  geom_boxplot(notch = TRUE) + 
  ggtitle("Daytime UHI July-Aug 2019") +
  theme_minimal() +
  scale_y_continuous(labels = function(x) paste0(x, "°")) +
  theme(legend.position = "none", axis.text=element_text(size=8)) + 
  xlab(NULL) + 
  ylab("Temperature (°F)") +
  coord_flip()
```



### Daytime Average Land Surface Temperatures and Populations of Concern

```{r LSTDayTempPop, include=FALSE}
# Pop Weighted avg of 24hr LST exposure for all Groups in New England relative to NE average
LSTDayPop <- ne_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  # filter(STATE == "Massachusetts") %>%
  dplyr::select(totalpopE,
                nhwhitepopE,
                minorityE,
                nhblackpopE,
                nhamerindpopE,
                nhasianpopE,
                nhnativhpopE,
                nhotherpopE,
                nh2morepopE,
                hisppopE,
                povknownE,
                num2povE, 
                eng_hhE,
                eng_limitE,
                age25upE,
                lthsE, 
                allAgesE, 
                under5E, 
                over64E, 
                meanDayLST) %>% 
  gather(key = Group, value = Pop, totalpopE:over64E) %>% 
  group_by(Group) %>% 
  summarize(LSTDaywMean = weighted.mean(x = meanDayLST, w = Pop, na.rm = TRUE),
            LSTDayMean = mean(meanDayLST, na.rm = TRUE)) %>% 
  filter(Group %in% c("minorityE","eng_limitE","num2povE","lthsE","under5E","over64E")) %>% 
  mutate(Group = case_when(
    Group == "minorityE" ~ "Minority",
    Group == "eng_limitE" ~ "Limited English HH",
    Group == "num2povE" ~ "Low Income",
    Group == "lthsE" ~ "No HS Dip",
    Group == "under5E" ~ "Under 5",
    Group == "over64E" ~ "Over 64"))


# Pop Weighted avg of LST day for all Groups in New England relative to NE average
LSTDayPopTab <- ne_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  # filter(STATE == "Massachusetts") %>%
  dplyr::select(totalpopE,
                nhwhitepopE,
                minorityE,
                nhblackpopE,
                nhamerindpopE,
                nhasianpopE,
                nhnativhpopE,
                nhotherpopE,
                nh2morepopE,
                hisppopE,
                povknownE,
                num2povE, 
                eng_hhE,
                eng_limitE,
                age25upE,
                lthsE, 
                allAgesE, 
                under5E, 
                over64E, 
                meanDayLST) %>% 
  gather(key = Group, value = Pop, totalpopE:over64E) %>% 
  group_by(Group) %>% 
  summarize(LSTDaywMean = weighted.mean(x = meanDayLST, w = Pop, na.rm = TRUE),
            LSTDayMean = mean(meanDayLST, na.rm = TRUE)) %>% 
  spread(key = Group, value = LSTDaywMean) %>% 
  transmute(Minority = (minorityE/LSTDayMean - 1)*100,
         #Minority_NHW = (minorityE/nhwhitepopE - 1)*100,
         `Limited English HH` = (eng_limitE/LSTDayMean - 1)*100,
         `Low Income` = (num2povE/LSTDayMean - 1)*100,
         `No HS Dip` = (lthsE/LSTDayMean - 1)*100,
         `Under 5` = (under5E/LSTDayMean - 1)*100,
         `Over 64` = (over64E/LSTDayMean - 1)*100) %>%
  gather(key = Group, value = Pct) %>% 
  left_join(x = LSTDayPop, y = ., by = "Group") %>% 
  mutate(Difference = LSTDaywMean - LSTDayMean) %>% 
  transmute(Group = Group,
            LSTDaywMean = paste0(round(LSTDaywMean,2),"°"),
            Difference = if_else(Difference > 0,
                                 paste0("+",round(Difference,2),"°"),
                                 paste0(round(Difference,2),"°")),
            Pct = if_else(Pct > 0, paste0("+",round(Pct,2),"%"),
                       paste0(round(Pct,2),"%")))
```

In addition to variations in the general geography of Land Surface Temperatures (LST), temperature exposure also varies demographically. Figure \@ref(fig:plotLSTDayTempPopAvgNE) and Table \@ref(tab:tabLSTDayPopAvgNE) show population-weighted exposures for populations of concern relative to the average daytime average LST for the region of `r round(mean(ne_blkgrpLST_sf$meanDayLST,na.rm=TRUE),1)`°F. For example, limited English speaking households in New England were exposed to a population-weighted daytime average LST of `r round(max(LSTDayPop$LSTDaywMean),1)`°, approximately `r round(max(LSTDayPop$LSTDaywMean) - max(LSTDayPop$LSTDayMean),1)`° (`r LSTDayPopTab[1,4]`) above the regional daytime average LST. Similarly, minorities were exposed to a population-weighted daytime average LST of approximately `r LSTDayPop %>% filter(Group == "Minority") %>% summarize(round(max(LSTDaywMean),1)) %>% pull()`°, `r round(LSTDayPop %>% filter(Group == "Minority") %>% summarize(max(LSTDaywMean)) %>% pull() - LSTDayPop %>% summarize(max(LSTDayMean)) %>% pull(),1)`° (`r LSTDayPopTab[3,4]`) above the regional daytime average LST. By contrast, persons over age 64 were, on average, exposed to population-weighted daytime average LSTs of `r LSTDayPop %>% filter(Group == "Over 64") %>% summarize(round(max(LSTDaywMean),1)) %>% pull()`°, approximately `r LSTDayPopTab[5,3]` below the regional average.

```{r plotLSTDayTempPopAvgNE, fig.align = "center", fig.cap="Population-weighted average exposures to average daytime Land Surface Temperatures for populations of concern in New England relative to the regional average."}
LSTDayPop %>% 
  ggplot(aes(x = reorder(Group,LSTDaywMean), 
             y = LSTDaywMean)) +
  geom_segment(aes(x = reorder(Group,LSTDaywMean), 
                   xend = reorder(Group,LSTDaywMean),
                   y = LSTDayMean, yend = LSTDaywMean), 
               color = "firebrick1") +
  geom_point(color = "firebrick4", size = 4, alpha = 0.8) +
  coord_flip() + 
  labs(x = "", y = "", title = "Population-Weighted Temperature Exposure to\nDaytime Average LST (°F)") + 
  theme_light() +
  theme(panel.grid.major.y = element_blank(),
        panel.border = element_blank(),
        axis.ticks.y = element_blank()) +
  geom_text(aes(x = Group, y = LSTDaywMean + 0.2 * sign(LSTDaywMean), 
                label = paste0(round(LSTDaywMean,1),"°")), 
            hjust = 1.6, vjust = -.8, size = 3,
            color=rgb(100,100,100, maxColorValue=255)) +
  scale_y_continuous(labels = function(x) paste0(x, "°")) +
  geom_hline(yintercept = LSTDayPop$LSTDayMean, linetype = "dashed") +
  geom_text(aes(x = "Low Income", y = 94.5, label = "Above regional\naverage"),
            color = "gray48") +
  geom_text(aes(x = "Low Income", y = 88, label = "Below regional\naverage"),
            color = "gray48") +
  # geom_segment(aes(x = "Under 5", xend = "Under 5", y = 92, yend = 95), 
  #              arrow = arrow(length = unit(0.3,"cm"))) +
  expand_limits(y = c(87,94))
```


```{r tabLSTDayPopAvgNE, fig.align = "center", fig.cap="Population-weighted average exposures to Daytime LST for populations of concern in New England relative to the regional average."}
LSTDayPopTab %>% 
  arrange(Group) %>% 
  kableExtra::kable(longtable = T, booktabs = T, digits = 2, 
                    col.names = c("Group","Avg Day LST (°F)",
                                  "Difference from Regional Avg* (°F)",
                                  "Pct Above/Below Regional Avg*"), 
                    caption = "Population-Weighted Temperature Exposure - Daytime Avg LST", align = "r") %>% 
  kableExtra::column_spec(2:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header")) %>% 
  kableExtra::footnote(., symbol = paste0("Regional average daytime LST is ", round(mean(ne_blkgrpLST_sf$meanDayLST,na.rm=TRUE),2), "°"))
```

```{r UHIDayPop, include=FALSE}
# Pop Weighted avg of 24hr LST exposure for all Groups in New England relative to NE average
UHIDayPop <- ne_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  # filter(STATE == "Massachusetts") %>%
  dplyr::select(totalpopE,
                nhwhitepopE,
                minorityE,
                nhblackpopE,
                nhamerindpopE,
                nhasianpopE,
                nhnativhpopE,
                nhotherpopE,
                nh2morepopE,
                hisppopE,
                povknownE,
                num2povE, 
                eng_hhE,
                eng_limitE,
                age25upE,
                lthsE, 
                allAgesE, 
                under5E, 
                over64E, 
                UHIDay) %>% 
  gather(key = Group, value = Pop, totalpopE:over64E) %>% 
  group_by(Group) %>% 
  summarize(UHIDaywMean = weighted.mean(x = UHIDay, w = Pop, na.rm = TRUE),
            UHIDayMean = mean(UHIDay, na.rm = TRUE)) %>% 
  filter(Group %in% c("minorityE","eng_limitE","num2povE","lthsE","under5E","over64E")) %>% 
  mutate(Group = case_when(
    Group == "minorityE" ~ "Minority",
    Group == "eng_limitE" ~ "Limited English HH",
    Group == "num2povE" ~ "Low Income",
    Group == "lthsE" ~ "No HS Dip",
    Group == "under5E" ~ "Under 5",
    Group == "over64E" ~ "Over 64"))


# Pop Weighted avg of PM2.5 for all Groups in New England relative to NE average
UHIDayPopTab <- ne_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  # filter(STATE == "Massachusetts") %>%
  dplyr::select(totalpopE,
                nhwhitepopE,
                minorityE,
                nhblackpopE,
                nhamerindpopE,
                nhasianpopE,
                nhnativhpopE,
                nhotherpopE,
                nh2morepopE,
                hisppopE,
                povknownE,
                num2povE, 
                eng_hhE,
                eng_limitE,
                age25upE,
                lthsE, 
                allAgesE, 
                under5E, 
                over64E, 
                UHIDay) %>% 
  gather(key = Group, value = Pop, totalpopE:over64E) %>% 
  group_by(Group) %>% 
  summarize(UHIDaywMean = weighted.mean(x = UHIDay, w = Pop, na.rm = TRUE),
            UHIDayMean = mean(UHIDay, na.rm = TRUE)) %>% 
  spread(key = Group, value = UHIDaywMean) %>% 
  transmute(Minority = (minorityE/UHIDayMean - 1)*100,
         #Minority_NHW = (minorityE/nhwhitepopE - 1)*100,
         `Limited English HH` = (eng_limitE/UHIDayMean - 1)*100,
         `Low Income` = (num2povE/UHIDayMean - 1)*100,
         `No HS Dip` = (lthsE/UHIDayMean - 1)*100,
         `Under 5` = (under5E/UHIDayMean - 1)*100,
         `Over 64` = (over64E/UHIDayMean - 1)*100) %>%
  gather(key = Group, value = Pct) %>% 
  left_join(x = UHIDayPop, y = ., by = "Group") %>% 
  mutate(Difference = UHIDaywMean - UHIDayMean) %>% 
  transmute(Group = Group,
            UHIDaywMean = paste0(round(UHIDaywMean,2),"°"),
            Difference = if_else(Difference > 0,
                                 paste0("+",round(Difference,2),"°"),
                                 paste0(round(Difference,2),"°")),
            Pct = if_else(Pct > 0, paste0("+",round(Pct,2),"%"),
                       paste0(round(Pct,2),"%")))
```

Figure \@ref(fig:plotUHIDayTempPopAvgNE) and Table \@ref(tab:tabUHIDayPopAvgNE) show population-weighted exposures for populations of concern relative to the average daytime average Urban Heat Island (UHI) effect for the region. All populations of concern are exposed to the UHI effect on average. Indeed, the average UHI effect across the region is approximately `r round(mean(ne_blkgrpLST_sf$UHIDay,na.rm=TRUE),1)`°. However, some populations experience significantly higher UHI effects. For example, limited English speaking households in New England were exposed to a population-weighted daytime average UHI effect of `r round(max(UHIDayPop$UHIDaywMean),1)`° above the rural average. This value is more than `r round(max(UHIDayPop$UHIDaywMean) - max(UHIDayPop$UHIDayMean),1)`° (`r UHIDayPopTab[1,4]`) above the regional daytime average UHI. Similarly, minorities were exposed to a population-weighted daytime average UHI effect of `r UHIDayPop %>% filter(Group == "Minority") %>% summarize(round(max(UHIDaywMean),1)) %>% pull()`°, approximately `r round(UHIDayPop %>% filter(Group == "Minority") %>% summarize(max(UHIDaywMean)) %>% pull() - max(UHIDayPop$UHIDayMean),1)`° (`r UHIDayPopTab[3,4]`) above the regional daytime average UHI. By contrast, persons over age 64 were, on average, exposed to population-weighted daytime average UHI of `r UHIDayPop %>% filter(Group == "Over 64") %>% summarize(round(max(UHIDaywMean),1)) %>% pull()`° (`r UHIDayPopTab[5,4]`) below the regional average.

```{r plotUHIDayTempPopAvgNE, echo=FALSE, message=FALSE, warning=FALSE, fig.align = "center", fig.cap="Population-weighted average exposures to daytime Urban Heat Island effect for populations of concern in New England."}
UHIDayPop %>% 
  ggplot(aes(x = reorder(Group,UHIDaywMean), 
             y = UHIDaywMean)) +
  geom_segment(aes(x = reorder(Group,UHIDaywMean), 
                   xend = reorder(Group,UHIDaywMean),
                   y = 0, yend = UHIDaywMean), 
               color = "firebrick1") +
  geom_point(color = "firebrick4", size = 4, alpha = 0.8) +
  coord_flip() + 
  labs(x = "", y = "", title = "Population-Weighted Temperature Exposure to\nDaytime Average UHI (°F)") + 
  theme_light() +
  theme(panel.grid.major.y = element_blank(),
        panel.border = element_blank(),
        axis.ticks.y = element_blank()) +
  geom_text(aes(x = Group, y = UHIDaywMean + 0.2 * sign(UHIDaywMean), 
                label = paste0(round(UHIDaywMean,1),"°")), 
            hjust = 1.6, vjust = -.8, size = 3,
            color=rgb(100,100,100, maxColorValue=255)) +
  scale_y_continuous(labels = function(x) paste0(x, "°")) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  # geom_text(aes(x = "Under 5", y = 16, label = "Above\nrural average"),
  #           color = "gray48") +
  # geom_segment(aes(x = "Over 64", xend = "Over 64", y = 12.5, yend = 16), 
  #              arrow = arrow(length = unit(0.3,"cm"))) +
  expand_limits(y = c(-1,17))
```


```{r tabUHIDayPopAvgNE, fig.align = "center", fig.cap="Population-weighted average exposures to daytime Urban Heat Island effect for populations of concern in New England relative to the rural average."}
UHIDayPopTab %>% 
  arrange(Group) %>% 
  kableExtra::kable(longtable = T, booktabs = T, digits = 2, 
                    col.names = c("Group","Avg Daytime UHI (°F)",
                                  "Difference from Regional UHI Avg* (°F)",
                                  "Pct Above/Below Regional UHI Avg*"), 
                    caption = "Population-Weighted Temperature Exposure - Daytime Avg UHI", align = "r") %>% 
  kableExtra::column_spec(2:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header")) %>% 
  kableExtra::footnote(., symbol = paste0("Regional average daytime UHI effect is ", round(mean(ne_blkgrpLST_sf$UHIDay,na.rm=TRUE),2), "°"))
```

Because the risk of heat-related illness is a combination of both exposure and vulnerability, it is important to consider where these factors come together. Figures \@ref(fig:mapLSTDayHighRiskLEH) and \@ref(fig:mapLSTDayHighRiskMIN) highlight Census Block Groups across New England which exhibited daytime average LSTs in the top quintile (i.e., 80th percentile) and also contained high percentages of limited English households and minorities in the top quintile for the region.[^Dayquintile] These high risk locations were present across five New England states in the late summer of 2019. Tables \@ref(tab:tabLSTDayHighRiskLEH) and \@ref(tab:tabLSTDayHighRiskMIN) show breakdowns of how many households these represent for each state, and the average daytime LST values for those same Census Block Groups. Note that Massachusetts contains the largest number of these at-risk households, over 73,000, and at-risk minority persons, over 690,000, while Rhode Island shows the highest average daytime LST for at-risk housholds and minorities, at 104°. Maps and tables for other populations of concern can be found in Figures \@ref(fig:mapLSTDayHighRiskLI) to \@ref(fig:mapLSTDayHighRiskU5) and Tables \@ref(tab:tabLSTDayHighRiskLI) to \@ref(tab:tabLSTDayHighRiskU5) in Appendix B.

```{r mapLSTDayHighRiskLEH, fig.align = "center", fig.cap="Map of Block Groups in the highest quintile (80th percentile) of daytime average LST and 80th percentile of percentage of Limited English Speaking Households by Census Block Group, July-Aug 2019."}
LSTDayEngLang5 <- ne_blkgrpLST_sf %>%
  filter(LSTDayEngLangScore == 5)

bbox_new <- st_bbox(LSTDayEngLang5) # current bounding box

xrange <- bbox_new$xmax - bbox_new$xmin # range of x values
yrange <- bbox_new$ymax - bbox_new$ymin # range of y values

# bbox_new[1] <- bbox_new[1] - (0.5 * xrange) # xmin - left
bbox_new[3] <- bbox_new[3] + (0.2 * xrange) # xmax - right
# bbox_new[2] <- bbox_new[2] - (0.5 * yrange) # ymin - bottom
# bbox_new[4] <- bbox_new[4] + (0.2 * yrange) # ymax - top

bbox_new <- bbox_new %>%  # take the bounding box ...
  st_as_sfc() # ... and make it a sf polygon

 m <- LSTDayEngLang5 %>%
  tm_shape(., unit = "mi", bbox = bbox_new) + tm_fill(col = "red", 
                        title = "80th Percentile Temp and Pop",
                        legend.is.portrait = TRUE) +
  tm_shape(ne_states_sf_cb) + tm_polygons(alpha = 0.4, lwd = 0.2) + 
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ne_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", xmod = 0.7, ymod = 0.2) +
  tm_scale_bar(breaks = c(0, 50, 100), text.size = 0.5, 
               position = c(0.4,0.005)) +
   tm_add_legend(type = "fill", 
                labels = "80th Percentile Block Groups",
                col = "red",
                border.col = "white", border.alpha = 0) +
  tm_layout(title = "Highest Exposure\nand Vulnerability\nto Daytime Avg\nTemps for Limited\nEnglish Speaking\nHouseholds\nJuly-Aug 2019", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"),
            legend.hist.width = 0.9)
 
 tmap_save(m, "DATA/LST/ne_LSTDayEngLang5.png",
          height = 5, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/ne_LSTDayEngLang5.png")
```


```{r tabLSTDayHighRiskLEH, fig.align = "center", fig.cap="Limited English Housheolds in Census Block Groups at or above 80th percentile of daytime average Land Surface Temperature."}
ne_blkgrpLST_sf %>%
  as.data.frame() %>%
  left_join(., statepopsdf,by="STATE") %>%
  filter(LSTDayEngLangScore == 5) %>%
  group_by(STATE) %>%
  summarize(`Limited English Households` = sum(eng_limitE),
            `Pct of Limited English Households` = round(sum(eng_limitE)/max(TotalLEH)*100,1),
            Min = min(meanDayLST,na.rm=TRUE),
            Mean = mean(meanDayLST,na.rm=TRUE),
            Max = max(meanDayLST,na.rm=TRUE)) %>% 
  mutate(`Pct of Limited English Households` = paste0(`Pct of Limited English Households`,"%"),
         Min = paste0(round(Min,2),"°"),
         Mean = paste0(round(Mean,2),"°"),
         Max = paste0(round(Max,2),"°")) %>% 
  kableExtra::kable(longtable = F, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Limited English Households in Highest Quintile of Daytime Avg LST", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header")) %>% 
  kableExtra::column_spec(2:3, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 3, "Daytime LST (°F)" = 3))
```

```{r mapLSTDayHighRiskMIN, fig.align = "center", fig.cap="Map of Block Groups in the highest quintile (80th percentile) of daytime average LST and 80th percentile of percentage of Minorities by Census Block Group, July-Aug 2019."}
LSTMinority5 <- ne_blkgrpLST_sf %>%
  filter(LSTDayMinorScore == 5)

bbox_new <- st_bbox(LSTMinority5) # current bounding box

xrange <- bbox_new$xmax - bbox_new$xmin # range of x values
yrange <- bbox_new$ymax - bbox_new$ymin # range of y values

# bbox_new[1] <- bbox_new[1] - (0.5 * xrange) # xmin - left
bbox_new[3] <- bbox_new[3] + (0.5 * xrange) # xmax - right
# bbox_new[2] <- bbox_new[2] - (0.5 * yrange) # ymin - bottom
bbox_new[4] <- bbox_new[4] + (0.25 * yrange) # ymax - top

bbox_new <- bbox_new %>%  # take the bounding box ...
  st_as_sfc() # ... and make it a sf polygon

 m <- LSTMinority5 %>%
  tm_shape(., unit = "mi", bbox = bbox_new) + tm_fill(col = "red", 
                        title = "80th Percentile Temp and Pop",
                        legend.is.portrait = TRUE) +
  tm_shape(ne_states_sf_cb) + tm_polygons(alpha = 0.4, lwd = 0.2) +
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ne_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", xmod = 0.7, ymod = 0.2) +
  tm_scale_bar(breaks = c(0, 50, 100), text.size = 0.5, 
               position = c(0.4,0.005)) +
   tm_add_legend(type = "fill", 
                labels = "80th Percentile Block Groups",
                col = "red",
                border.col = "white", border.alpha = 0) + 
  tm_layout(title = "Highest Exposure\nand Vulnerability\nto Daytime Avg\nTemps for\nMinorities\nJuly-Aug 2019", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"),
            legend.hist.width = 0.9)
 
 tmap_save(m, "DATA/LST/ne_LSTDayMinority5.png",
          height = 5, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/ne_LSTDayMinority5.png")
```

\pagebreak
```{r tabLSTDayHighRiskMIN, fig.align = "center", fig.cap="Minorities in Census Block Groups at or above 80th percentile of daytime average Land Surface Temperature."}
ne_blkgrpLST_sf %>%
  as.data.frame() %>%
  left_join(., statepopsdf,by="STATE") %>%
  filter(LSTDayMinorScore == 5) %>%
  group_by(STATE) %>%
  summarize(Minorities = sum(minorityE),
            `Pct of Minorities` = round(sum(minorityE)/max(TotalMin)*100,1),
            Min = min(meanDayLST,na.rm=TRUE),
            Mean = mean(meanDayLST,na.rm=TRUE),
            Max = max(meanDayLST,na.rm=TRUE)) %>% 
  mutate(`Pct of Minorities` = paste0(`Pct of Minorities`,"%"),
         Min = paste0(round(Min,2),"°"),
         Mean = paste0(round(Mean,2),"°"),
         Max = paste0(round(Max,2),"°")) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Minorities in Highest Quintile of Daytime Avg LST", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header")) %>% 
  kableExtra::column_spec(2:3, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 3, "Daytime LST (°F)" = 3))
```

\pagebreak
## Nighttime Average Land Surface Temperatures

High nighttime temperatures carry added significance for concerns about heat risk. Research has repeatedly shown that exposure to high nighttime temperatures over several days increases the probability of death during a heat wave in urban areas, especially for the elderly.[^Gabriel2011]

Nighttime Land Surface Temperatures (LST) represent LST values collected during the night (12am - 3:06am) over eight days, from July 28 to August 4, 2019.  Nighttime Land Surface Temperatures (LST) varied significantly across New England, with highest temperatures apparent in and around the most urbanized areas of the region, particularly along the I-95 corridor from Stamford to New Haven,Connecticut, and along the I-91 corridor from New Haven through Hartford and Springfield, Massachusetts. Other high temperature areas are apparent around greater Providence, Rhode Island, and the greater Boston region out to I-495, and reaching into southern New Hampshire (see Figure \@ref(fig:mapLSTnightNE) below). In contrast to daytime LSTs, nighttime values are high in pockets outside of the urban cores, including areas such as the Champlain Valley in northwestern Vermont and Cape Cod in Massachusetts.

```{r mapLSTnightNE, fig.align = "center", fig.cap="Map of July-August 2019 average nighttime Land Surface Temperatures (LST) across New England at Census Block Group level.", strip.white=TRUE}
# Map of average LST across New England
m <- tm_shape(ne_blkgrpLST_sf, unit = "mi") + 
  tm_fill("meanNightLST", style = "quantile", palette = "PuRd",
          title = expression("Temperature " ( degree*F)),
          legend.hist = FALSE,
          colorNA = NULL,
          textNA = NULL,
          legend.format=list(list(digits=1)),
          legend.is.portrait = TRUE) +
  tm_shape(ne_states_sf_cb) + tm_borders(alpha = 0.4, lwd = 0.5) + 
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ne_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", xmod = 0.7, ymod = 0.2) +
  tm_scale_bar(breaks = c(0, 50, 100), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "Average Night\nLand Surface\nTemperatures\nJuly-Aug 2019", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"),
            legend.hist.width = 0.9)

tmap_save(m, "DATA/LST/ne_NightLST.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/ne_NightLST.png")
```

```{r LSTNightoutliers, include=FALSE}
# download municipal polygons
# nh_towns_sf <- county_subdivisions("NH") %>%
#   dplyr::select(NAME) %>%
#   st_transform(., crs = st_crs(ne_blkgrpLST_sf))
# 
# me_towns_sf <- county_subdivisions("ME") %>%
#   dplyr::select(NAME) %>%
#   st_transform(., crs = st_crs(ne_blkgrpLST_sf))

ri_towns_sf <- county_subdivisions("RI") %>%
  dplyr::select(NAME) %>%
  st_transform(., crs = st_crs(ne_blkgrpLST_sf))

vt_towns_sf <- county_subdivisions("VT") %>%
  dplyr::select(NAME) %>%
  st_transform(., crs = st_crs(ne_blkgrpLST_sf))

# Identify municipalities with high heat outliers in New Hampshire and Maine
nh_LSTNightoutliers <- ne_blkgrpLST_sf %>% 
  filter(STATE == "New Hampshire") %>% 
  filter(meanNightLST == max(meanNightLST, na.rm = TRUE)) %>% 
  st_centroid(.) %>% 
  st_intersection(nh_towns_sf,.) %>% 
  filter(row_number() == 1) %>% 
  pull(NAME)

me_LSTNightoutliers <- ne_blkgrpLST_sf %>% 
  filter(STATE == "Maine") %>% 
  filter(meanNightLST == max(meanNightLST, na.rm = TRUE)) %>% 
  st_centroid(.) %>% 
  st_intersection(me_towns_sf,.) %>% 
  filter(row_number() == 1) %>% 
  pull(NAME)

ri_LSTNightoutliers <- ne_blkgrpLST_sf %>% 
  filter(STATE == "Rhode Island") %>% 
  filter(meanNightLST == max(meanNightLST, na.rm = TRUE)) %>% 
  st_centroid(.) %>% 
  st_intersection(ri_towns_sf,.) %>% 
  filter(row_number() == 1) %>% 
  pull(NAME)

vt_LSTNightoutliers <- ne_blkgrpLST_sf %>% 
  filter(STATE == "Vermont") %>% 
  filter(meanNightLST == max(meanNightLST, na.rm = TRUE)) %>% 
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  filter(row_number() == 1) %>% 
  pull(NAME)
```

Rhode Island, Connecticut, and Massachusetts are the warmest states in New England, with Rhode Island standing out as statistically significantly warmer (see Figure \@ref(fig:boxplotLSTnightNE)). The warmest nighttime temperatures for New England occurred in `r ri_LSTNightoutliers`, Rhode Island. By contrast, Vermont, Maine, and New Hampshire are significantly cooler overall, with Vermont exhibiting the coolest nighttime temperatures in New England. Notably all three northern New England states showed significant outliers for high nighttime temperatures, particularly in `r nh_LSTNightoutliers`, New Hampshire, `r me_LSTNightoutliers`, Maine, and `r vt_LSTNightoutliers`, Vermont. 

```{r boxplotLSTnightNE, fig.align = "center", fig.cap="Boxplot of July-Aug 2019 nighttime Land Surface Temperatures (LST) across New England at Census Block Group level."}
# Test for statistically significant differences
# First, set up a one-way ANOVA
# amod <- ne_blkgrpLST_sf %>% 
#   mutate(STATE = as.factor(STATE)) %>% 
#   aov(meanNightLST ~ STATE, data = .)
# library(multcomp)
# # Next, set up all-pair comparisons for factor `STATE`. See https://stats.stackexchange.com/questions/341463/test-for-significant-differences-between-groups-in-r for example
# g <- glht(amod, linfct = mcp(STATE = "Tukey"))
# summary(g)

# boxplot of PM2.5 by state for 2016
ne_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  drop_na(meanNightLST) %>% 
  ggplot(aes(x = reorder(STATE, meanNightLST, FUN = "median", .desc = TRUE), 
             y = meanNightLST, fill = STATE)) + 
  geom_boxplot(notch = TRUE) + 
  ggtitle("Average Nighttime Land Surface Temperatures July-Aug 2019") +
  theme_minimal() +
  scale_y_continuous(labels = function(x) paste0(x, "°")) +
  theme(legend.position = "none", axis.text=element_text(size=8)) + 
  xlab(NULL) + 
  ylab("Temperature (°F)") +
  coord_flip()
```


### Nighttime Urban Heat Island Effect

The UHI effect for nighttime average temperatures in New England are most pronounced in the more densely urbanized areas (see Figure \@ref(fig:mapUHInightNE)), particularly around Bridgeport and Hartford, Connecticut, around Boston, Massachusetts. Small pockets of high UHI are also evident around outside of the region's major cities. Temperatures in these areas are up to `r round(max(ne_blkgrpLST_sf$UHINight,na.rm = TRUE),1)`° warmer than the rural background average. By contrast, most of Vemont (outside of the Lake Champlain area), northern New Hampshire, and most of Maine are cooler than rural background average. 

```{r mapUHInightNE, fig.align = "center", fig.cap="Map of July-Aug 2019 nighttime Urban Heat Island (UHI) effect across New England at Census Block Group level. UHI is the difference between Block Group average temperature and average temperature of rural areas.", strip.white=TRUE}
# Map of average LST UHI across New England
m <- tm_shape(ne_blkgrpLST_sf, unit = "mi") + 
  tm_fill("UHINight", palette = "-RdYlGn", midpoint = 0,
          legend.reverse = TRUE,
          title = expression("Temperature " ( degree*F)),
          legend.hist = FALSE,
          colorNA = NULL,
          textNA = NULL,
          legend.format=list(list(digits=1)),
          legend.is.portrait = TRUE) +
  tm_shape(ne_states_sf_cb) + tm_borders(alpha = 0.4, lwd = 0.5) + 
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  # tm_shape(ne_urban_sf) + tm_fill(alpha = 0.5) +
  tm_shape(ne_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", xmod = 0.7, ymod = 0.2) +
  tm_scale_bar(breaks = c(0, 50, 100), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "Nighttime\nUrban Heat Island\nJuly-Aug 2019", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"),
            legend.hist.width = 0.9)

tmap_save(m, "DATA/LST/ne_NightUHI.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/ne_NightUHI.png")
```

Rhode Island, Connecticut, and Massachusetts exhibit the highest UHI effects (see Figure \@ref(fig:boxplotUHINightNE)). In these states, median UHI effect is range from approximately `r round(ne_blkgrpLST_sf %>% as.data.frame() %>% filter(STATE == "Massachusetts") %>% summarize(median(UHINight,na.rm = TRUE)) %>% pull(),1)`° in Massachusetts to `r round(ne_blkgrpLST_sf %>% as.data.frame() %>% filter(STATE == "Rhode Island") %>% summarize(median(UHINight,na.rm = TRUE)) %>% pull(),1)`° above the rural background. Northern New England states show statistically significantly lower UHI effect values. Vermont's median value is below 0°, indicating cooler nighttime temperatures than the rural background temperature.

```{r boxplotUHINightNE, fig.align = "center", fig.cap="Boxplot of July-Aug 2019 nighttime Urban Heat Island (UHI) across New England at Census Block Group level. UHI is the difference between Block Group average temperature and average temperature of rural areas."}
# Test for statistically significant differences
# First, set up a one-way ANOVA
# amod <- ne_blkgrpLST_sf %>% 
#   mutate(STATE = as.factor(STATE)) %>% 
#   aov(UHINight ~ STATE, data = .)
# library(multcomp)
# # Next, set up all-pair comparisons for factor `STATE`. See https://stats.stackexchange.com/questions/341463/test-for-significant-differences-between-groups-in-r for example
# g <- glht(amod, linfct = mcp(STATE = "Tukey"))
# summary(g)

# boxplot of PM2.5 by state for 2016
ne_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  drop_na(UHINight) %>% 
  ggplot(aes(x = reorder(STATE, UHINight, FUN = "median", .desc = TRUE), 
             y = UHINight, fill = STATE)) + 
  geom_boxplot(notch = TRUE) + 
  ggtitle("Nighttime UHI July-Aug 2019") +
  theme_minimal() +
  scale_y_continuous(labels = function(x) paste0(x, "°")) +
  theme(legend.position = "none", axis.text=element_text(size=8)) + 
  xlab(NULL) + 
  ylab("Temperature (°F)") +
  coord_flip()
```


### Nighttime Average Land Surface Temperatures and Populations of Concern

```{r LSTNightPop, include=FALSE}
# Pop Weighted avg of 24hr LST exposure for all Groups in New England relative to NE average
LSTNightPop <- ne_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  # filter(STATE == "Massachusetts") %>%
  dplyr::select(totalpopE,
                nhwhitepopE,
                minorityE,
                nhblackpopE,
                nhamerindpopE,
                nhasianpopE,
                nhnativhpopE,
                nhotherpopE,
                nh2morepopE,
                hisppopE,
                povknownE,
                num2povE, 
                eng_hhE,
                eng_limitE,
                age25upE,
                lthsE, 
                allAgesE, 
                under5E, 
                over64E, 
                meanNightLST) %>% 
  gather(key = Group, value = Pop, totalpopE:over64E) %>% 
  group_by(Group) %>% 
  summarize(LSTNightwMean = weighted.mean(x = meanNightLST, w = Pop, na.rm = TRUE),
            LSTNightMean = mean(meanNightLST, na.rm = TRUE)) %>% 
  filter(Group %in% c("minorityE","eng_limitE","num2povE","lthsE","under5E","over64E")) %>% 
  mutate(Group = case_when(
    Group == "minorityE" ~ "Minority",
    Group == "eng_limitE" ~ "Limited English HH",
    Group == "num2povE" ~ "Low Income",
    Group == "lthsE" ~ "No HS Dip",
    Group == "under5E" ~ "Under 5",
    Group == "over64E" ~ "Over 64"))

# Pop Weighted avg of PM2.5 for all Groups in New England relative to NE average
LSTNightPopTab <- ne_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  # filter(STATE == "Massachusetts") %>%
  dplyr::select(totalpopE,
                nhwhitepopE,
                minorityE,
                nhblackpopE,
                nhamerindpopE,
                nhasianpopE,
                nhnativhpopE,
                nhotherpopE,
                nh2morepopE,
                hisppopE,
                povknownE,
                num2povE, 
                eng_hhE,
                eng_limitE,
                age25upE,
                lthsE, 
                allAgesE, 
                under5E, 
                over64E, 
                meanNightLST) %>% 
  gather(key = Group, value = Pop, totalpopE:over64E) %>% 
  group_by(Group) %>% 
  summarize(LSTNightwMean = weighted.mean(x = meanNightLST, w = Pop, na.rm = TRUE),
            LSTNightMean = mean(meanNightLST, na.rm = TRUE)) %>% 
  spread(key = Group, value = LSTNightwMean) %>% 
  transmute(Minority = (minorityE/LSTNightMean - 1)*100,
         #Minority_NHW = (minorityE/nhwhitepopE - 1)*100,
         `Limited English HH` = (eng_limitE/LSTNightMean - 1)*100,
         `Low Income` = (num2povE/LSTNightMean - 1)*100,
         `No HS Dip` = (lthsE/LSTNightMean - 1)*100,
         `Under 5` = (under5E/LSTNightMean - 1)*100,
         `Over 64` = (over64E/LSTNightMean - 1)*100) %>%
  gather(key = Group, value = Pct) %>% 
  left_join(x = LSTNightPop, y = ., by = "Group") %>% 
  mutate(Difference = LSTNightwMean - LSTNightMean) %>% 
  transmute(Group = Group,
            LSTNightwMean = paste0(round(LSTNightwMean,2),"°"),
            Difference = if_else(Difference > 0,
                                 paste0("+",round(Difference,2),"°"),
                                 paste0(round(Difference,2),"°")),
            Pct = if_else(Pct > 0, paste0("+",round(Pct,2),"%"),
                       paste0(round(Pct,2),"%")))
```

In addition to variations in the general geography of Land Surface Temperatures (LST), temperature exposure also varies demographically. Figure \@ref(fig:plotLSTNightTempPopAvgNE) and Table \@ref(tab:tabLSTNightPopAvgNE) show population-weighted exposures for populations of concern relative to the average nighttime average LST for the region of `r round(mean(ne_blkgrpLST_sf$meanNightLST,na.rm=TRUE),1)`°F. For example, limited English speaking households in New England were exposed to a population-weighted nighttime average LST of `r round(max(LSTNightPop$LSTNightwMean),1)`°, approximately `r round(max(LSTNightPop$LSTNightwMean) - max(LSTNightPop$LSTNightMean),1)`° (`r LSTNightPopTab[1,4]`) above the regional nighttime average LST. Similarly, minorities were exposed to a population-weighted nighttime average LST of `r LSTNightPop %>% filter(Group == "Minority") %>% summarize(round(max(LSTNightwMean),1)) %>% pull()`°, over `r round(LSTNightPop %>% filter(Group == "Minority") %>% summarize(max(LSTNightwMean)) %>% pull() - max(LSTNightPop$LSTNightMean),1)`° (`r LSTNightPopTab[3,4]`) above the regional nighttime average LST. By contrast, persons over age 64 were, on average, exposed to population-weighted nighttime average LSTs `r LSTNightPopTab[5,3]` below the regional average.

```{r plotLSTNightTempPopAvgNE, fig.align = "center", fig.cap="Population-weighted average exposures to average nighttime Land Surfact Temperatures for populations of concern in New England relative to the regional average."}
LSTNightPop %>% 
  ggplot(aes(x = reorder(Group,LSTNightwMean), 
             y = LSTNightwMean)) +
  geom_segment(aes(x = reorder(Group,LSTNightwMean), 
                   xend = reorder(Group,LSTNightwMean),
                   y = LSTNightMean, yend = LSTNightwMean), 
               color = "firebrick1") +
  geom_point(color = "firebrick4", size = 4, alpha = 0.8) +
  coord_flip() + 
  labs(x = "", y = "", title = "Population-Weighted Temperature Exposure to\nNighttime Average LST (°F)") + 
  theme_light() +
  theme(panel.grid.major.y = element_blank(),
        panel.border = element_blank(),
        axis.ticks.y = element_blank()) +
  geom_text(aes(x = Group, y = LSTNightwMean + 0.2 * sign(LSTNightwMean), 
                label = paste0(round(LSTNightwMean,1),"°")), 
            hjust = 1.6, vjust = -.8, size = 3,
            color=rgb(100,100,100, maxColorValue=255)) +
  scale_y_continuous(labels = function(x) paste0(x, "°")) +
  geom_hline(yintercept = LSTNightPop$LSTNightMean, linetype = "dashed") +
  geom_text(aes(x = "Low Income", y = 66.8, label = "Above regional\naverage"),
            color = "gray48") +
  geom_text(aes(x = "Low Income", y = 65.3, label = "Below regional\naverage"),
            color = "gray48") +
  # geom_segment(aes(x = "Under 5", xend = "Under 5", y = 66.2, yend = 67), 
  #              arrow = arrow(length = unit(0.3,"cm"))) +
  expand_limits(y = c(65,68))
```


```{r tabLSTNightPopAvgNE, fig.align = "center", fig.cap="Population-weighted average exposures to nighttime LST for populations of concern in New England relative to the regional average."}
LSTNightPopTab %>% 
  arrange(Group) %>% 
  kableExtra::kable(longtable = T, booktabs = T, digits = 2, 
                    col.names = c("Group","Avg Night LST (°F)",
                                  "Difference from Regional Avg* (°F)",
                                  "Pct Above/Below Regional Avg*"), 
                    caption = "Population-Weighted Temperature Exposure - Nighttime Avg LST", align = "r") %>% 
  kableExtra::column_spec(2:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header")) %>% 
  kableExtra::footnote(., symbol = paste0("Regional average nighttime LST is ", round(mean(ne_blkgrpLST_sf$meanNightLST,na.rm=TRUE),2), "°"))
```

```{r UHINightPop, include=FALSE}
# Pop Weighted avg of 24hr LST exposure for all Groups in New England relative to NE average
UHINightPop <- ne_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  # filter(STATE == "Massachusetts") %>%
  dplyr::select(totalpopE,
                nhwhitepopE,
                minorityE,
                nhblackpopE,
                nhamerindpopE,
                nhasianpopE,
                nhnativhpopE,
                nhotherpopE,
                nh2morepopE,
                hisppopE,
                povknownE,
                num2povE, 
                eng_hhE,
                eng_limitE,
                age25upE,
                lthsE, 
                allAgesE, 
                under5E, 
                over64E, 
                UHINight) %>% 
  gather(key = Group, value = Pop, totalpopE:over64E) %>% 
  group_by(Group) %>% 
  summarize(UHINightwMean = weighted.mean(x = UHINight, w = Pop, na.rm = TRUE),
            UHINightMean = mean(UHINight, na.rm = TRUE)) %>% 
  filter(Group %in% c("minorityE","eng_limitE","num2povE","lthsE","under5E","over64E")) %>% 
  mutate(Group = case_when(
    Group == "minorityE" ~ "Minority",
    Group == "eng_limitE" ~ "Limited English HH",
    Group == "num2povE" ~ "Low Income",
    Group == "lthsE" ~ "No HS Dip",
    Group == "under5E" ~ "Under 5",
    Group == "over64E" ~ "Over 64"))


# Pop Weighted avg of PM2.5 for all Groups in New England relative to NE average
UHINightPopTab <- ne_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  # filter(STATE == "Massachusetts") %>%
  dplyr::select(totalpopE,
                nhwhitepopE,
                minorityE,
                nhblackpopE,
                nhamerindpopE,
                nhasianpopE,
                nhnativhpopE,
                nhotherpopE,
                nh2morepopE,
                hisppopE,
                povknownE,
                num2povE, 
                eng_hhE,
                eng_limitE,
                age25upE,
                lthsE, 
                allAgesE, 
                under5E, 
                over64E, 
                UHINight) %>% 
  gather(key = Group, value = Pop, totalpopE:over64E) %>% 
  group_by(Group) %>% 
  summarize(UHINightwMean = weighted.mean(x = UHINight, w = Pop, na.rm = TRUE),
            UHINightMean = mean(UHINight, na.rm = TRUE)) %>% 
  spread(key = Group, value = UHINightwMean) %>% 
  transmute(Minority = (minorityE/UHINightMean - 1)*100,
         #Minority_NHW = (minorityE/nhwhitepopE - 1)*100,
         `Limited English HH` = (eng_limitE/UHINightMean - 1)*100,
         `Low Income` = (num2povE/UHINightMean - 1)*100,
         `No HS Dip` = (lthsE/UHINightMean - 1)*100,
         `Under 5` = (under5E/UHINightMean - 1)*100,
         `Over 64` = (over64E/UHINightMean - 1)*100) %>%
  gather(key = Group, value = Pct) %>% 
  left_join(x = UHINightPop, y = ., by = "Group") %>% 
  mutate(Difference = UHINightwMean - UHINightMean) %>% 
  transmute(Group = Group,
            UHINightwMean = paste0(round(UHINightwMean,2),"°"),
            Difference = if_else(Difference > 0,
                                 paste0("+",round(Difference,2),"°"),
                                 paste0(round(Difference,2),"°")),
            Pct = if_else(Pct > 0, paste0("+",round(Pct,2),"%"),
                       paste0(round(Pct,2),"%")))
```

Figure \@ref(fig:plotUHINightTempPopAvgNE) and Table \@ref(tab:tabUHINightPopAvgNE) show population-weighted exposures for populations of concern relative to the average daytime average Urban Heat Island (UHI) effect for the region. All populations of concern are exposed to the UHI effect on average. Indeed, the average UHI effect across the region is approximately `r round(mean(ne_blkgrpLST_sf$UHINight,na.rm=TRUE),1)`°. However, some populations experience significantly higher UHI effects. For example, limited English speaking households in New England were exposed to a population-weighted daytime average UHI effect of `r round(max(UHINightPop$UHINightwMean),1)`° above the rural average. This value is more than `r round(max(UHINightPop$UHINightwMean) - max(UHINightPop$UHINightMean),1)`° (`r UHINightPopTab[1,4]`) above the regional daytime average UHI. Similarly, minorities were exposed to a population-weighted daytime average UHI effect of `r UHINightPop %>% filter(Group == "Minority") %>% summarize(round(max(UHINightwMean),1)) %>% pull()`°, approximately `r round(UHINightPop %>% filter(Group == "Minority") %>% summarize(max(UHINightwMean)) %>% pull() - max(UHINightPop$UHINightMean),1)`° (`r UHINightPopTab[3,4]`) above the regional daytime average UHI. By contrast, persons over age 64 were, on average, exposed to population-weighted daytime average UHI of `r UHINightPop %>% filter(Group == "Over 64") %>% summarize(round(max(UHINightwMean),1)) %>% pull()`° (`r UHINightPopTab[5,4]`) below the regional average.

```{r plotUHINightTempPopAvgNE, echo=FALSE, message=FALSE, warning=FALSE, fig.align = "center", fig.cap="Population-weighted average exposures to nighttime Urban Heat Island effect for populations of concern in New England."}
UHINightPop %>% 
  ggplot(aes(x = reorder(Group,UHINightwMean), 
             y = UHINightwMean)) +
  geom_segment(aes(x = reorder(Group,UHINightwMean), 
                   xend = reorder(Group,UHINightwMean),
                   y = 0, yend = UHINightwMean), 
               color = "firebrick1") +
  geom_point(color = "firebrick4", size = 4, alpha = 0.8) +
  coord_flip() + 
  labs(x = "", y = "", title = "Population-Weighted Temperature Exposure to\nNighttime Average UHI (°F)") + 
  theme_light() +
  theme(panel.grid.major.y = element_blank(),
        panel.border = element_blank(),
        axis.ticks.y = element_blank()) +
  geom_text(aes(x = Group, y = UHINightwMean + 0.2 * sign(UHINightwMean), 
                label = paste0(round(UHINightwMean,1),"°")), 
            hjust = 1.6, vjust = -.8, size = 3,
            color=rgb(100,100,100, maxColorValue=255)) +
  scale_y_continuous(labels = function(x) paste0(x, "°")) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  # geom_text(aes(x = "Under 5", y = 8.5, label = "Above\nrural average"),
  #           color = "gray48") +
  # geom_segment(aes(x = "Over 64", xend = "Over 64", y = 6.5, yend = 8.5), 
  #              arrow = arrow(length = unit(0.3,"cm"))) +
  expand_limits(y = c(-1,10))
```


```{r tabUHINightPopAvgNE, fig.align = "center", fig.cap="Population-weighted average exposures to nighttime Urban Heat Island effect for populations of concern in New England relative to the rural average."}
UHINightPopTab %>% 
  arrange(Group) %>% 
  kableExtra::kable(longtable = T, booktabs = T, digits = 2, 
                    col.names = c("Group","Avg Nighttime UHI (°F)",
                                  "Difference from Regional UHI Avg* (°F)",
                                  "Pct Above/Below Regional UHI Avg*"), 
                    caption = "Population-Weighted Temperature Exposure - Nighttime Avg UHI", align = "r") %>% 
  kableExtra::column_spec(2:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header")) %>% 
  kableExtra::footnote(., symbol = paste0("Regional average nighttime UHI effect is ", round(mean(ne_blkgrpLST_sf$UHINight,na.rm=TRUE),2), "°"))
```

Because the risk of heat-related illness is a combination of both exposure and vulnerability, it is important to consider where these factors come together. Figure \@ref(fig:mapLSTNightHighRiskLHE) highlights Census Block Groups across New England which exhibited nighttime average LSTs in the top quintile (i.e., 80th percentile) and also contained high percentages of limited English speaking households in the top quintile for the region.[^Nightquintile] These high risk locations were present across three New England states in the late summer of 2019. Table \@ref(tab:tabLSTNightHighRiskLHE) shows a breakdown of how many households these represent for each state, and the average nighttime LST values for those same Census Block Groups. Note that Massachusetts contains the largest number of these at-risk households, over 62,000, while Rhode Island shows the highest average nighttime LST for at-risk housholds, at 98.4°. Maps and tables for other populations of concern can be found in Figures \@ref(fig:mapLSTNightHighRiskLI) to \@ref(fig:mapLSTNightHighRiskU5) and Tables \@ref(tab:tabLSTNightHighRiskLI) to \@ref(tab:tabLSTNightHighRiskU5) in Appendix B.

```{r mapLSTNightHighRiskLHE, fig.align = "center", fig.cap="Map of Block Groups in the highest quintile (80th percentile) of nighttime average LST and 80th percentile of percentage of Limited English Speaking Households by Census Block Group, July-Aug 2019."}
LSTNightEngLang5 <- ne_blkgrpLST_sf %>%
  filter(LSTNightEngLangScore == 5)

bbox_new <- st_bbox(LSTNightEngLang5) # current bounding box

xrange <- bbox_new$xmax - bbox_new$xmin # range of x values
yrange <- bbox_new$ymax - bbox_new$ymin # range of y values

# bbox_new[1] <- bbox_new[1] - (0.5 * xrange) # xmin - left
bbox_new[3] <- bbox_new[3] + (0.25 * xrange) # xmax - right
# bbox_new[2] <- bbox_new[2] - (0.5 * yrange) # ymin - bottom
bbox_new[4] <- bbox_new[4] + (0.25 * yrange) # ymax - top

bbox_new <- bbox_new %>%  # take the bounding box ...
  st_as_sfc() # ... and make it a sf polygon

m <- LSTNightEngLang5 %>%
  tm_shape(., unit = "mi", bbox = bbox_new) + tm_fill(col = "red", 
                        title = "80th Percentile Temp and Pop",
                        legend.is.portrait = TRUE) +
  tm_shape(ne_states_sf_cb) + tm_polygons(alpha = 0.4, lwd = 0.2) + 
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ne_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", xmod = 0.7, ymod = 0.2) +
  tm_scale_bar(breaks = c(0, 50, 100), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_add_legend(type = "fill", 
                labels = "80th Percentile Block Groups",
                col = "red",
                border.col = "white", border.alpha = 0) + 
  tm_layout(title = "Highest Exposure\nand Vulnerability\nto Nighttime Avg\nTemps for Limited\nEnglish Speaking\nHouseholds\nJuly-Aug 2019", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"),
            legend.hist.width = 0.9)

tmap_save(m, "DATA/LST/ne_LSTNightEngLang5.png",
          height = 5, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/ne_LSTNightEngLang5.png")
```


```{r tabLSTNightHighRiskLHE, fig.align = "center", fig.cap="Limited English Speaking Housheolds in Census Block Groups at or above 80th percentile of nighttime average Land Surface Temperature."}
ne_blkgrpLST_sf %>%
  as.data.frame() %>%
  left_join(., statepopsdf,by="STATE") %>%
  filter(LSTNightEngLangScore == 5) %>%
  group_by(STATE) %>%
  summarize(`Limited English Households` = sum(eng_limitE),
            `Pct of Limited English Households` = round(sum(eng_limitE)/max(TotalLEH)*100,1),
            Min = min(meanNightLST,na.rm=TRUE),
            Mean = mean(meanNightLST,na.rm=TRUE),
            Max = max(meanNightLST,na.rm=TRUE)) %>% 
  mutate(`Pct of Limited English Households` = paste0(`Pct of Limited English Households`,"%"),
         Min = paste0(round(Min,2),"°"),
         Mean = paste0(round(Mean,2),"°"),
         Max = paste0(round(Max,2),"°")) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Limited English Households in Highest Quintile of Nighttime Avg LST", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header")) %>% 
  kableExtra::column_spec(2:3, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 3, "Nighttime LST (°F)" = 3))
```

```{r mapLSTNightHighRiskMIN, fig.align = "center", fig.cap="Map of Block Groups in the highest quintile (80th percentile) of nighttime average LST and 80th percentile of percentage of Minorities by Census Block Group, July-Aug 2019."}
LSTNightMinority5 <- ne_blkgrpLST_sf %>%
  filter(LSTNightMinorScore == 5)

bbox_new <- st_bbox(LSTNightMinority5) # current bounding box

xrange <- bbox_new$xmax - bbox_new$xmin # range of x values
yrange <- bbox_new$ymax - bbox_new$ymin # range of y values

# bbox_new[1] <- bbox_new[1] - (0.5 * xrange) # xmin - left
bbox_new[3] <- bbox_new[3] + (0.2 * xrange) # xmax - right
# bbox_new[2] <- bbox_new[2] - (0.5 * yrange) # ymin - bottom
# bbox_new[4] <- bbox_new[4] + (0.25 * yrange) # ymax - top

bbox_new <- bbox_new %>%  # take the bounding box ...
  st_as_sfc() # ... and make it a sf polygon

m <- LSTNightMinority5 %>%
  tm_shape(., unit = "mi", bbox = bbox_new) + tm_fill(col = "red", 
                        title = "80th Percentile Temp and Pop",
                        legend.is.portrait = TRUE) +
  tm_shape(ne_states_sf_cb) + tm_polygons(alpha = 0.4, lwd = 0.2) + 
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ne_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", xmod = 0.7, ymod = 0.2) +
  tm_scale_bar(breaks = c(0, 50, 100), text.size = 0.5, 
               position = c(0.4,0.005)) +
  tm_add_legend(type = "fill", 
                labels = "80th Percentile Block Groups",
                col = "red",
                border.col = "white", border.alpha = 0) + 
  tm_layout(title = "Highest Exposure\nand Vulnerability\nto Nighttime Avg\nTemps for\nMinorities\nJuly-Aug 2019", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"),
            legend.hist.width = 0.9)

tmap_save(m, "DATA/LST/ne_LSTNightMinority5.png",
          height = 5, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/ne_LSTNightMinority5.png")
```


```{r tabLSTNightHighRiskMIN, fig.align = "center", fig.cap="Minorities in Census Block Groups at or above 80th percentile of Nighttime average Land Surface Temperature."}
ne_blkgrpLST_sf %>%
  as.data.frame() %>%
  left_join(., statepopsdf,by="STATE") %>%
  filter(LSTNightMinorScore == 5) %>%
  group_by(STATE) %>%
  summarize(Minorities = sum(minorityE),
            `Pct of Minorities` = round(sum(minorityE)/max(TotalMin)*100,1),
            Min = min(meanNightLST,na.rm=TRUE),
            Mean = mean(meanNightLST,na.rm=TRUE),
            Max = max(meanNightLST,na.rm=TRUE)) %>% 
  mutate(`Pct of Minorities` = paste0(`Pct of Minorities`,"%"),
         Min = paste0(round(Min,2),"°"),
         Mean = paste0(round(Mean,2),"°"),
         Max = paste0(round(Max,2),"°")) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Minorities in Highest Quintile of Nighttime Avg LST", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header")) %>% 
  kableExtra::column_spec(2:3, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 3, "Nighttime LST (°F)" = 3))
```

[^NWShazards]: See National Weather Service. Weather Related Fatality and Injury Statistics. https://www.weather.gov/hazstat/

[^CDCpicture]: See CDC. Picture of America: Heat-Related Illness Fact Sheet. https://www.cdc.gov/pictureofamerica/pdfs/picture_of_america_heat-related_illness.pdf

[^Samuelson2020]: See Holly Samuelson, Amir Baniassadi, Anne Lin, Pablo Izaga González, Thomas Brawley, and Tushar Narula. 2020. "Housing as a critical determinant of heat vulnerability and health." *Science of The Total Environment* 720. https://doi.org/10.1016/j.scitotenv.2020.137296.

[^Williams2020]: See Augusta A. Williams, Joseph G. Allen, Paul J. Catalano, and John D. Spengler. 2020. "The Role of Individual and Small-Area Social and Environmental Factors on Heat Vulnerability to Mortality Within and Outside of the Home in Boston, MA." *Climate* 8(2): 29. https://doi.org/10.3390/cli8020029

[^Hondula2015]: See D. M. Hondula, R. E. Davis, M. V. Saha, C. R. Wegner, and L. M. Veazey. 2015. "Geographic dimensions of heat-related mortality in seven U.S. cities." *Environmental research* 138: 439–452. https://doi.org/10.1016/j.envres.2015.02.033; Francesco Nordio, Antonella Zanobetti, Elena Colicino, Itai Kloog, and Joel Schwartz. 2015.
"Changing patterns of the temperature–mortality association by time and location in the US, and implications for climate change." *Environment International* 81: 80-86. https://doi.org/10.1016/j.envint.2015.04.009; G. A. Wellenius, M. N. Eliot, K. F. Bush, D. Holt, , R. A. Lincoln,  A. E. Smith, and J. Gold. 2017. "Heat-related morbidity and mortality in New England: Evidence for local policy." *Environmental research* 156: 845–853. https://doi.org/10.1016/j.envres.2017.02.005.

[^IPCCsummary]: See IPCC. 2018. "Summary for Policymakers," in *Global Warming of 1.5°C. An IPCC Special Report on the impacts of global warming of 1.5°C above pre-industrial levels and related global greenhouse gas emission pathways, in the context of strengthening the global response to the threat of climate change, sustainable development, and efforts to eradicate poverty.* V. Masson-Delmotte, P. Zhai, H.-O. Pörtner, D. Roberts, J. Skea, P.R. Shukla, A. Pirani, W. Moufouma-Okia, C. Péan, R. Pidcock, S. Connors, J.B.R. Matthews, Y. Chen, X. Zhou, M.I. Gomis, E. Lonnoy, T. Maycock, M. Tignor, and T. Waterfield (eds.). World Meteorological Organization, Geneva, Switzerland, 32 pp. https://www.ipcc.ch/sr15/chapter/spm/

[^4thNatlAssessment]: See USGCRP. 2018. Impacts, Risks, and Adaptation in the United States: Fourth National Climate Assessment, Volume II. Reidmiller, D.R., C.W. Avery, D.R. Easterling, K.E. Kunkel, K.L.M. Lewis, T.K. Maycock, and B.C. Stewart (eds.). U.S. Global Change Research Program, Washington, DC, USA, 1515 pp. doi: 10.7930/NCA4.2018. https://nca2018.globalchange.gov/chapter/18/

[^MODIS11]: NASA MODIS Land Surface Temperature and Emissivity (MOD11). https://modis.gsfc.nasa.gov/data/dataprod/mod11.php

[^Kloog2014]: See Itai Kloog, Francesco Nordio, Brent A. Coull, and Joel Schwartz. 2014. "Predicting spatiotemporal mean air temperature using MODIS satellite surface temperature measurements across the Northeastern USA." *Remote Sensing of Environment* 150:132-139. https://doi.org/10.1016/j.rse.2014.04.024; Thanh Noi Phan and Martin Kappas. 2018. "Application of MODIS land surface temperature data: a systematic literature review and analysis." *Journal of Applied Remote Sensing* 12(4): 041501. https://doi.org/10.1117/1.JRS.12.041501.

[^Zhou2019]: See  D. Zhou, J. Xiao,  S. Bonafoni, C. Berger,  K. Deilami, Y. Zhou, S. Frolking, R. Yao, Z. Qiao, and J.A. Sobrino.  2019. "Satellite Remote Sensing of Surface Urban Heat Islands: Progress, Challenges, and Perspectives." *Remote Sens.* 11:48. 

[^urban]: Urban areas are defined by the U.S. Census as Census Blocks meeting specific population density thresholds, as well as land use. See U.S. Census Bureau. Urban and Rural. https://www.census.gov/programs-surveys/geography/guidance/geo-areas/urban-rural.html

[^LST24quintile]: The top quintile, or 80th percentile, represents the value that is greater than 80% of all values in the data set. In the case of day-night average Land Surface Temperatures (LSTs) for New England, a temperature of `r round(quantile(ne_blkgrpLST_sf$meanAvgLST,0.8,na.rm=TRUE),2)`° is greater than 80% of all LST values in the region. Similarly for the percentage of limited English speaking households, a Census Block Group with `r round(quantile(ne_blkgrpLST_sf$eng_limit_pctE,0.8,na.rm=TRUE),2)`% of households classified as limited English speaking respresents a percentage that is higher than 80% of all other Block Groups in the region. 

[^Dayquintile]: The top quintile, or 80th percentile, represents the value that is greater than 80% of all values in the data set. In the case of daytime average Land Surface Temperatures (LSTs) for New England, a temperature of `r round(quantile(ne_blkgrpLST_sf$meanDayLST,0.8,na.rm=TRUE),2)`° is greater than 80% of all LST values in the region. Similarly for the percentage of limited English speaking households, a Census Block Group with `r round(quantile(ne_blkgrpLST_sf$eng_limit_pctE,0.8,na.rm=TRUE),2)`% of households classified as limited English speaking respresents a percentage that is higher than 80% of all other Block Groups in the region. 

[^Gabriel2011]: See Katharina M.A. Gabriel and Wilfried R. Endlicher. 2011. "Urban and rural mortality rates during heat waves in Berlin and Brandenburg, Germany." *Environmental Pollution* 159(8–9): 2044-2050. https://doi.org/10.1016/j.envpol.2011.01.016; Karine Laaidi, Abdelkrim Zeghnoun, Bénédicte Dousset, Philippe Bretin, Stéphanie Vandentorren, Emmanuel Giraudet, and Pascal Beaudeau. 2012. "The Impact of Heat Islands on Mortality in Paris during the August 2003 Heat Wave." *Environmental Health Perspectives* 120:2. https://doi.org/10.1289/ehp.1103532; Peninah Murage, Shakoor Hajat, and R. Sari Kovats. 2017. "Effect of night-time temperatures on cause and age-specific mortality in London." *Environmental Epidemiology* 1(2):e005. 

[^Nightquintile]: The top quintile, or 80th percentile, represents the value that is greater than 80% of all values in the data set. In the case of daytime average Land Surface Temperatures (LSTs) for New England, a temperature of `r round(quantile(ne_blkgrpLST_sf$meanNightLST,0.8,na.rm=TRUE),2)`° is greater than 80% of all LST values in the region. Similarly for the percentage of limited English speaking households, a Census Block Group with `r round(quantile(ne_blkgrpLST_sf$eng_limit_pctE,0.8,na.rm=TRUE),2)`% of households classified as limited English speaking respresents a percentage that is higher than 80% of all other Block Groups in the region.


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

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

* MODIS
  + Daytime Land Surface Temperature
  + Nighttime Land Surface Temperature
  + Day-Night Average Land Surface Temperature
  + Validation of Land Surface Temperature (LST) values
  + Urban Heat Island effect
* American Community Survey 5-year Estimates
  + Population demographics

## MODIS {-}

The heat risk analysis presented here is based on Land Surface Temperatures (LST) derived from NASA's Moderate-resolution Imaging Spectroradiometer (MODIS) satellite sensor.[^MODIS11] The MODIS sensor is carried aboard the Terra and Aqua satellites operataed by the U.S. National Aeronautics and Space Administration (NASA). These satellites view the entire earth every 1 to 2 days as they circuit the earth in polar orbits (i.e., orbiting north-south over the poles). They are synchronized so that Terra crosses the earth's equator (from north to south) in the morning, while Aqua crosses it in the afternoon (from south to north). Any given area of land is therefore imaged four times each day. MODIS captures images of the earth's surface in 36 spectral bands, or wavelengths of light, ranging from 0.4$\mu$m (ultraviolet) to 14.4$\mu$m (infrared). Land Surface Temperature (LST) values are derived from MODIS bands 31 and 32 (10.78$\mu$m to 12.27$\mu$m), which capture thermal infrared energy (i.e., heat energy) emitted by the earth's surface. These images are corrected for atmospheric effects, such as water vapor and atmospheric density. 

MODIS/Aqua Land Surface Temperature/Emissivity 8-Day L3 Global 1 km SIN Grid data (MOD11A2 v006) were downloaded through the USGS Earth Explorer interface (https://earthexplorer.usgs.gov/). Specifically, two scenes - h12v4 and h13v4 - covering northeastern North America for the period July 28 to August 4, 2019 were acquired in HDF format. Day and night LST scientific data sets (SDS) were extracted from the HDF files, converted to raster format, and reprojected to WGS NAD83. LST raster values were converted from Kelvin to Fahrenheit using the formula: °F = °K*9/5 - 459.67

The LST rasters were then mosaiced to create a single scene for each time period and clipped to New England. 

### Daytime Land Surface Temperature {-}

Daytime Land Surface Temperatures represent LST values collected during the day (11:48am - 2pm) for the period July 28 to August 4, 2019. Values are presented in degrees of Fahrenheit (°F). 

### Nighttime Land Surface Temperature {-}

Nighttime Land Surface Temperatures represent LST values collected during the night (12am - 3:06am) for the period July 28 to August 4, 2019. Values are presented in degrees of Fahrenheit (°F). 

### Day-Night Average Land Surface Temperature {-}

Day-Night Average Land Surface Temperatures (LST) represent a simple average of LST values collected during the day (11:48am - 2pm) and during the night (12am - 3:06am) over eight days, from July 28 to August 4, 2019. Values are presented in degrees of Fahrenheit (°F). 

### Validation of Land Surface Temperature (LST) values {-}
```{r loadISD, include=FALSE}
# load the previously processed data
isd_ne_DayNight_sf <- readRDS("DATA/LST/isd_ne_DayNight_sf.Rds")
```

In order to validate land surface temperature (LST) values, MODIS LST data was compared to ground-based air temperature measurements for the same time periods across New England. Weather data was acquired from NOAA's Integrated Surface Data ('ISD') files from the NOAA National Centers for Environmental Information (https://www.ncdc.noaa.gov/isd/data-access) using the Commmon Access tool (https://www.ncei.noaa.gov/access/search/data-search/global-hourly). This yielded `r nrow(isd_ne_DayNight_sf)` weather stations across New England with air temperature data available for the same period as the MODIS LST data (i.e., July 28 - August 4, 2019; see Figure \@ref(fig:ISDstations)). 

```{r ISDstations, fig.align = "center", fig.cap="NOAA Integrated Surface Data Weather Stations across New England."}
tm_shape(ne_states_sf_cb, unit = "mi") + tm_borders(lwd = 0.2) +
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(isd_ne_DayNight_sf) + tm_dots(shape = 4, size = 0.05) +
   tm_scale_bar(breaks = c(0, 50, 100), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "NOAA ISD Weather Stations",
            frame = FALSE, main.title.size = 0.8,
            title.position = c("right", "top"),
            inner.margins = c(0.02,0.02,0.1,0.02))
```

Air temperature data for the specific hours covered by the MODIS data were extracted. Weather station air temperature values were then compared to coincident LST pixel values for the same time periods.   

The validation analysis shows that LST values are positively correlated with weather station air temperatures (see Table \@ref(tab:tabCorrelations)). 

```{r tabCorrelations}
# Compare day average LST to weather station air temp data for same time period
cor_day <- cor.test(x = isd_ne_DayNight_sf$gridDay, y = isd_ne_DayNight_sf$tempFday, use = "pairwise.complete.obs", method = "spearman")
cor_day_rho <- cor_day$estimate
cor_day_p <- cor_day$p.value

# Compare night average LST to weather station air temp data for same time period
cor_night <- cor.test(x = isd_ne_DayNight_sf$gridNight, y = isd_ne_DayNight_sf$tempFnight, use = "pairwise.complete.obs", method = "spearman")
cor_night_rho <- cor_night$estimate
cor_night_p <- cor_night$p.value

# Compare day-night average LST to weather station air temp data for same time period
cor_daynight <- cor.test(x = isd_ne_DayNight_sf$gridDayNight, y = isd_ne_DayNight_sf$tempFdaynight, use = "pairwise.complete.obs", method = "spearman")
cor_daynight_rho <- cor_daynight$estimate
cor_daynight_p <- cor_daynight$p.value

# Assemble into table
Time <- c("Day","Night","Day-Night")
Rho <- c(cor_day_rho,cor_night_rho,cor_daynight_rho)
pvalue <- c(cor_day_p,cor_night_p,cor_daynight_p)
data.frame(Time,Rho,pvalue) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 3,
                    col.names = c("Time of Day","Spearman's Rho","p-value"),
                    caption = "Spearman's Rho Correlations between Weather Stations and MODIS LST", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header"))
```

```{r lmModels, include=FALSE}
# Compare day average LST to weather station air temp data for same time period
lm_day <- lm(tempFday ~ gridDay, data = isd_ne_DayNight_sf)
RMSE_day <- sqrt(mean(lm_day$residuals^2))

# Compare night average LST to weather station air temp data for same time period
lm_night <- lm(tempFnight ~ gridNight, data = isd_ne_DayNight_sf)
RMSE_night <- sqrt(mean(lm_night$residuals^2))

# Compare day-night average LST to weather station air temp data for same time period
lm_daynight <- lm(tempFdaynight ~ gridDayNight, data = isd_ne_DayNight_sf)
RMSE_daynight <- sqrt(mean(lm_daynight$residuals^2))
```

Linear models show LST values are significantly predictive of air temperatures (see Figure \@ref(fig:LSTvalidation)). Night and Day-Night Average LST values are better predictors of air temperature than Day LST. The latter has a Root Mean Square Error (RMSE) of `r round(RMSE_day,2)`°, while Day-Night Avg LST has a RMSE of `r round(RMSE_daynight,2)`°. Most of the prediction error is due to one especially low temperature from a single weather station. These prediction error levels are consistent with published research on the use of MODIS LST as a proxy for air temperature.

```{r LSTvalidation, warning=FALSE, message=FALSE, fig.align = "center", fig.cap="Comparison of MODIS LST values and air temperatures for same time periods from NOAA Integrated Surface Data Weather Stations across New England."}
day <- isd_ne_DayNight_sf %>% 
  ggplot(., aes(x = gridDay, y = tempFday)) + 
  geom_point(alpha = 0.5, size = 3, col = "firebrick4") + 
  theme_light() +
  scale_y_continuous(labels = function(x) paste0(x, "°")) +
  scale_x_continuous(labels = function(x) paste0(x, "°")) +
  # ggtitle("Day LST vs Weather Station Air Temp") +
  xlab("MODIS LST Day") + ylab("Weather Station Day Air Temp") +
  geom_smooth(method = "lm") +
  annotate("text",label=paste0("R-squared = ", round(summary(lm_day)$r.squared,3),"\nRMSE = ",round(RMSE_day,2)),x=85,y=61)

night <- isd_ne_DayNight_sf %>% 
  ggplot(., aes(x = gridNight, y = tempFnight)) + 
  geom_point(alpha = 0.5, size = 3, col = "darkorchid4") + 
  theme_light() +
  scale_y_continuous(labels = function(x) paste0(x, "°")) +
  scale_x_continuous(labels = function(x) paste0(x, "°")) +
  # ggtitle("Night LST vs Weather Station Air Temp") +
  xlab("MODIS LST Night") + ylab("Weather Station Night Air Temp") +
  geom_smooth(method = "lm") +
  annotate("text",label=paste0("R-squared = ", round(summary(lm_night)$r.squared,3),"\nRMSE = ",round(RMSE_night,2)),x=63,y=55)

daynight <- isd_ne_DayNight_sf %>% 
  ggplot(., aes(x = gridDayNight, y = tempFdaynight)) + 
  geom_point(alpha = 0.5, size = 3, col = "darkorange4") + 
  theme_light() +
  scale_y_continuous(labels = function(x) paste0(x, "°")) +
  scale_x_continuous(labels = function(x) paste0(x, "°")) +
  # ggtitle("Day-Night Avg LST vs Weather Station Air Temp") +
  xlab("MODIS LST Day-Night") + ylab("Weather Station Day-Night Air Temp") +
  geom_smooth(method = "lm") +
  annotate("text",label=paste0("R-squared = ", round(summary(lm_daynight)$r.squared,3),"\nRMSE = ",round(RMSE_daynight,2)),x=75,y=58)

library(gridExtra)
# grid.arrange(day,night,daynight,nrow=1,top="LST vs Weather Station Air Temp")

g <- arrangeGrob(day,night,daynight,ncol=1,
                 top="LST vs Weather Station Air Temp",
                 heights = unit(rep(2.5,3),rep("inches",3)),
                 widths = unit(5,"inches"))
ggsave(g, filename = "DATA/LST/ne_validate.png", 
       width = 5, height = 8, units = "in")
knitr::include_graphics("DATA/LST/ne_validate.png")
```

 

### Urban Heat Island effects {-}

Urban Heat Island (UHI) effects are measured as the difference between LST raster values across New England for each time period and the average LST for rural regions for the same time periods. It is represented by the following formula:

*UHI effect = LST - Mean Rural LST*

For example, to calculate the UHI effects for daytime temperatures, daytime LST pixels falling within the rural regions of New England were all averaged to provide a mean rural LST benchmark value. This mean rural LST value was then subtracted from the daytime LST values of all pixels across the region. The resulting value for each pixel indicates the magnitude or size of the difference, and the sign (+/-) indicates whether it is warmer or cooler than the mean rural LST.

Rural regions of New England were identified according to the U.S. Census Bureau's urban-rural classification. The Census Bureau identifies two types of urban areas:

* Urbanized Areas (UAs) of 50,000 or more people;
* Urban Clusters (UCs) of at least 2,500 and less than 50,000 people.

“Rural” encompasses all population, housing, and territory not included within an urban area.[^urban] 

Urban areas TIGER/Line spatial files were downloaded from the Census Bureau via API using the `tigris` package in R. Areas of New England not within the urban areas polygons were then extracted to represent rural regions. LST raster pixels falling within these non-urban or rual polygons were extracted and a mean LST value calculated for each time period (i.e., day, night, and day-night average) to represent the mean rural LST benchmark for calculating the UHI effect.


## 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",
                                  "Minority",
                                  "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 Minority",
                                  "RI Income",
                                  "RI Minority"),
                     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 Minority threshold: 25% or more are persons of Hispanic or Latino origin or persons who are not White",
                                     "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 Minority threshold: highest 15% of block groups with persons of Hispanic or Latino origin or persons who are not White"),
                     'ACS Table ID' = c("B03002: Hispanic or Latino Origin by Race",
                                "B03002: Hispanic or Latino Origin by Race",
                                "C17002: Ratio of Income to Poverty Level",
                                "C16002: Household Language by Household Limited English Speaking Status",
                                "B15002: Sex by Educational Attainment",
                                "B01001: Sex by Age",
                                "B01001: Sex by Age",
                                "B01001: Sex by Age",
                                "B18101: Sex by Age by Disability Status",
                                "B08201: Household Size by Vehicles Available",
                                "C17002: Ratio of Income to Poverty Level",
                                "B19001: Household Income",
                                "C16002: Household Language by Household Limited English Speaking Status",
                                "B03002: Hispanic or Latino Origin by Race",
                                "C17002: Ratio of Income to Poverty Level",
                                "B03002: Hispanic or Latino Origin by Race"),
                     Geography = c("Block Group",
                                   "Block Group",
                                   "Block Group",
                                   "Block Group",
                                   "Block Group",
                                   "Block Group",
                                   "Block Group",
                                   "Block Group",
                                   "Tract",
                                   "Tract",
                                   "Block Group",
                                   "Block Group",
                                   "Block Group",
                                   "Block Group",
                                   "Block Group",
                                   "Block Group"))

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

## Population-weighted averages {-}

Wherever feasible, population exposure to temperature is reported as a population-weighted average. A population weighted-average is equivalent to a weighted mean in which the raw values for which a mean (or average) is calculated are multiplied by a weight factor. For example, we are interested in knowing whether Minorities are exposed to higher average temperatures than the general or Total Population. Consider the table below.

```{r wmean1}
data.frame(BG = c("BG1","BG2","BG3","BG4"),
           Temp = c(65,75,85,70),
           Minority = c(5,10,20,10),
           NonMinority = c(20,12,5,10)) %>% 
  mutate(TotalPop = Minority + NonMinority) %>% 
  kableExtra::kable(longtable = T, booktabs = T, 
                    caption = "Block Group Populations and Temperature", align = "l") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header"))
```

Population numbers of minorities, as well as the total population, and temperatures, are each reported by Block Group. Since each Block Group is associated with one temperature value, we might assume (incorrectly) that everyone is equally exposed to the average temperatures of all Block Groups - `r round(mean(c(65,75,85,70)),2)`. However, not all Block Groups have the same number or categories of people, which means that each temperature value is associated with a different number of people. Do Minorities occupy Block Groups with higher temperatures than the simple average would indicate?

To calculate the population-weighted average temperature exposure for Minorities, the number of Minorities in each Block Group is used as a 'weight'. The temperature of each Block Group is multiplied by its respective number of Minorities. See the table below. 

```{r wmean2}
data.frame(BG = c("BG1","BG2","BG3","BG4"),
           Temp = c(65,75,85,70),
           Minority = c(5,10,20,10),
           NonMinority = c(20,12,5,10)) %>% 
  mutate(TotalPop = Minority + NonMinority,
         TempxMinority = Temp*Minority) %>% 
  bind_rows(summarize_all(., funs(if(is.numeric(.)) sum(.) else "Total"))) %>% 
  kableExtra::kable(longtable = T, booktabs = T, 
                    caption = "Block Group Populations and Temperature", 
                    align = "l") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header")) %>% 
  kableExtra::row_spec(5, bold = TRUE) %>% 
  kableExtra::column_spec(6, bold = TRUE, background = "#d3d3d3")
```

The total or sum of the products (i.e., Minority pop x Temperature) is then divided by the sum of the weights (i.e., total Minority population), so that 3475/45 = `r round(3475/45,2)`. The result is a weighted average temperature that is influenced by the number of minorities. 

This process is repeated for the Total Population so that the two population-weighted average temperatures can be compared. Below is the calculation of population-weighted calculation for the total population. 

```{r wmean3}
data.frame(BG = c("BG1","BG2","BG3","BG4"),
           Temp = c(65,75,85,70),
           Minority = c(5,10,20,10),
           NonMinority = c(20,12,5,10)) %>% 
  mutate(TotalPop = Minority + NonMinority,
         TempxTotalPop = Temp*TotalPop) %>% 
  bind_rows(summarize_all(., funs(if(is.numeric(.)) sum(.) else "Total"))) %>% 
  kableExtra::kable(longtable = T, booktabs = T, 
                    caption = "Block Group Populations and Temperature", 
                    align = "l") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header")) %>% 
  kableExtra::row_spec(5, bold = TRUE) %>% 
  kableExtra::column_spec(6, bold = TRUE, background = "#d3d3d3")
```

For the TotalPop, the population-weighted average temperature is 6800/92 = `r round(6800/92,2)`. Thus we can see that Minorities experience a higher population-weighted average temperature than the general or total population. 

\pagebreak
# Appendix B: Supplementary Figures {-}

## Day-Night Highest Exposure Populations of Concern {-}
```{r mapLST24hrHighRiskLI, fig.align = "center", fig.cap="Map of Block Groups in the highest quintile (80th percentile) of day-night average LST and 80th percentile of percentage of Low Income persons by Census Block Group, July-Aug 2019."}
LST24LowInc5 <- ne_blkgrpLST_sf %>%
  filter(LST24LowIncScore == 5)

bbox_new <- st_bbox(LST24LowInc5) # current bounding box

xrange <- bbox_new$xmax - bbox_new$xmin # range of x values
yrange <- bbox_new$ymax - bbox_new$ymin # range of y values

# bbox_new[1] <- bbox_new[1] - (0.5 * xrange) # xmin - left
bbox_new[3] <- bbox_new[3] + (0.2 * xrange) # xmax - right
# bbox_new[2] <- bbox_new[2] - (0.5 * yrange) # ymin - bottom
bbox_new[4] <- bbox_new[4] + (0.2 * yrange) # ymax - top

bbox_new <- bbox_new %>%  # take the bounding box ...
  st_as_sfc() # ... and make it a sf polygon

m <- tm_shape(LST24LowInc5, unit = "mi" , bbox = bbox_new) + tm_fill(col = "red", title = "80th Percentile", legend.is.portrait = TRUE) +
  tm_shape(ne_states_sf_cb) + tm_polygons(alpha = 0.4, lwd = 0.2) +
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ne_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", xmod = 0.7, ymod = 0.2) +
  tm_scale_bar(breaks = c(0, 50, 100), text.size = 0.5, 
               position = c(0.4,0.005)) +
  tm_add_legend(type = "fill", 
                labels = "80th Percentile Block Groups",
                col = "red",
                border.col = "white", border.alpha = 0) +
  tm_layout(title = "Highest Exposure\nand Vulnerability\nto Day-Night Avg\nTemps for Low\nIncome Persons\nJuly-Aug 2019", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"),
            legend.hist.width = 0.9)

tmap_save(m, "DATA/LST/ne_LST24LowInc5.png",
          height = 5, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/ne_LST24LowInc5.png")
```


```{r tabLST24hrHighRiskLI, fig.align = "center", fig.cap="Low Income Persons in Census Block Groups at or above 80th percentile of day-night average Land Surface Temperature."}
ne_blkgrpLST_sf %>%
  as.data.frame() %>%
  left_join(., statepopsdf,by="STATE") %>% 
  filter(LST24LowIncScore == 5) %>%
  group_by(STATE) %>%
  summarize(`Low Income Persons` = sum(eng_limitE),
            `Pct of Low Income Persons` = round(sum(num2povE)/max(TotalLowInc)*100,1),
            Min = min(meanAvgLST,na.rm=TRUE),
            Mean = mean(meanAvgLST,na.rm=TRUE),
            Max = max(meanAvgLST,na.rm=TRUE)) %>% 
  mutate(`Pct of Low Income Persons` = paste0(`Pct of Low Income Persons`,"%"),
         Min = paste0(round(Min,2),"°"),
         Mean = paste0(round(Mean,2),"°"),
         Max = paste0(round(Max,2),"°")) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Low Income Persons in Highest Quintile of Day-Night Avg LST", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header")) %>% 
  kableExtra::column_spec(2:3, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 3, "Day-Night LST (°F)" = 3))
```


```{r mapLST24hrHighRiskNOHS, fig.align = "center", fig.cap="Map of Block Groups in the highest quintile (80th percentile) of day-night average LST and 80th percentile of percentage of adults with less than a high school diploma by Census Block Group, July-Aug 2019."}
LST24NoHSDip5 <- ne_blkgrpLST_sf %>%
  filter(LST24NoHSDipScore == 5)

bbox_new <- st_bbox(LST24NoHSDip5) # current bounding box

xrange <- bbox_new$xmax - bbox_new$xmin # range of x values
yrange <- bbox_new$ymax - bbox_new$ymin # range of y values

# bbox_new[1] <- bbox_new[1] - (0.5 * xrange) # xmin - left
bbox_new[3] <- bbox_new[3] + (0.2 * xrange) # xmax - right
# bbox_new[2] <- bbox_new[2] - (0.5 * yrange) # ymin - bottom
bbox_new[4] <- bbox_new[4] + (0.2 * yrange) # ymax - top

bbox_new <- bbox_new %>%  # take the bounding box ...
  st_as_sfc() # ... and make it a sf polygon

m <- tm_shape(LST24NoHSDip5, unit = "mi" , bbox = bbox_new) + tm_fill(col = "red", title = "80th Percentile", legend.is.portrait = TRUE) +
  tm_shape(ne_states_sf_cb) + tm_polygons(alpha = 0.4, lwd = 0.2) +
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ne_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", xmod = 0.7, ymod = 0.2) +
  tm_scale_bar(breaks = c(0, 50, 100), text.size = 0.5, 
               position = c(0.4,0.005)) +
  tm_add_legend(type = "fill", 
                labels = "80th Percentile Block Groups",
                col = "red",
                border.col = "white", border.alpha = 0) +
  tm_layout(title = "Highest Exposure\nand Vulnerability\nto Day-Night Avg\nTemps for Adults\nwithout a High\nSchool Diploma\nJuly-Aug 2019", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"),
            legend.hist.width = 0.9)

tmap_save(m, "DATA/LST/ne_LST24NoHSDip5.png",
          height = 5, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/ne_LST24NoHSDip5.png")
```


```{r tabLST24hrHighRiskNOHS, fig.align = "center", fig.cap="Adults without a High School Diploma in Census Block Groups at or above 80th percentile of day-night average Land Surface Temperature."}
ne_blkgrpLST_sf %>%
  as.data.frame() %>%
  left_join(., statepopsdf,by="STATE") %>% 
  filter(LST24NoHSDipScore == 5) %>%
  group_by(STATE) %>%
  summarize(`No HS Dip` = sum(lthsE),
            `Pct of No HS Dip` = round(sum(lthsE)/max(TotalNoHS)*100,1),
            Min = min(meanAvgLST,na.rm=TRUE),
            Mean = mean(meanAvgLST,na.rm=TRUE),
            Max = max(meanAvgLST,na.rm=TRUE)) %>% 
  mutate(`Pct of No HS Dip` = paste0(`Pct of No HS Dip`,"%"),
         Min = paste0(round(Min,2),"°"),
         Mean = paste0(round(Mean,2),"°"),
         Max = paste0(round(Max,2),"°")) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Adults without a High School Diploma in Highest Quintile of Day-Night Avg LST", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header")) %>% 
  kableExtra::column_spec(2:3, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 3, "Day-Night LST (°F)" = 3))
```


```{r mapLST24hrHighRiskO64, fig.align = "center", fig.cap="Map of Block Groups in the highest quintile (80th percentile) of day-night average LST and 80th percentile of percentage of persons over 64 by Census Block Group, July-Aug 2019."}
LST24Over645 <- ne_blkgrpLST_sf %>%
  filter(LST24Over64Score == 5)

bbox_new <- st_bbox(LST24Over645) # current bounding box

xrange <- bbox_new$xmax - bbox_new$xmin # range of x values
yrange <- bbox_new$ymax - bbox_new$ymin # range of y values

# bbox_new[1] <- bbox_new[1] - (0.5 * xrange) # xmin - left
bbox_new[3] <- bbox_new[3] + (0.2 * xrange) # xmax - right
# bbox_new[2] <- bbox_new[2] - (0.5 * yrange) # ymin - bottom
bbox_new[4] <- bbox_new[4] + (0.2 * yrange) # ymax - top

bbox_new <- bbox_new %>%  # take the bounding box ...
  st_as_sfc() # ... and make it a sf polygon

m <- tm_shape(LST24Over645, unit = "mi" , bbox = bbox_new) + tm_fill(col = "red", title = "80th Percentile", legend.is.portrait = TRUE) +
  tm_shape(ne_states_sf_cb) + tm_polygons(alpha = 0.4, lwd = 0.2) +
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ne_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", xmod = 0.7, ymod = 0.2) +
  tm_scale_bar(breaks = c(0, 50, 100), text.size = 0.5, 
               position = c(0.4,0.005)) +
  tm_add_legend(type = "fill", 
                labels = "80th Percentile Block Groups",
                col = "red",
                border.col = "white", border.alpha = 0) +
  tm_layout(title = "Highest Exposure\nand Vulnerability\nto Day-Night Avg\nTemps for Persons\nOver Age 64\nJuly-Aug 2019", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"),
            legend.hist.width = 0.9)

tmap_save(m, "DATA/LST/ne_LST24Over645.png",
          height = 5, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/ne_LST24Over645.png")
```


```{r tabLST24hrHighRiskO64, fig.align = "center", fig.cap="Persons over age 64 in Census Block Groups at or above 80th percentile of day-night average Land Surface Temperature."}
ne_blkgrpLST_sf %>%
  as.data.frame() %>%
  left_join(., statepopsdf,by="STATE") %>% 
  filter(LST24Over64Score == 5) %>%
  group_by(STATE) %>%
  summarize(`Over 64` = sum(over64E),
            `Pct of Over 64` = round(sum(over64E)/max(TotalOver64)*100,1),
            Min = min(meanAvgLST,na.rm=TRUE),
            Mean = mean(meanAvgLST,na.rm=TRUE),
            Max = max(meanAvgLST,na.rm=TRUE)) %>% 
  mutate(`Pct of Over 64` = paste0(`Pct of Over 64`,"%"),
         Min = paste0(round(Min,2),"°"),
         Mean = paste0(round(Mean,2),"°"),
         Max = paste0(round(Max,2),"°")) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Persons Over Age 64 in Highest Quintile of Day-Night Avg LST", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header")) %>% 
  kableExtra::column_spec(2:3, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 3, "Day-Night LST (°F)" = 3))
```


```{r mapLST24hrHighRiskU5, fig.align = "center", fig.cap="Map of Block Groups in the highest quintile (80th percentile) of day-night average LST and 80th percentile of percentage of children under 5 by Census Block Group, July-Aug 2019."}
LST24Under55 <- ne_blkgrpLST_sf %>%
  filter(LST24Under5Score == 5)

bbox_new <- st_bbox(LST24Under55) # current bounding box

xrange <- bbox_new$xmax - bbox_new$xmin # range of x values
yrange <- bbox_new$ymax - bbox_new$ymin # range of y values

# bbox_new[1] <- bbox_new[1] - (0.5 * xrange) # xmin - left
bbox_new[3] <- bbox_new[3] + (0.2 * xrange) # xmax - right
# bbox_new[2] <- bbox_new[2] - (0.5 * yrange) # ymin - bottom
bbox_new[4] <- bbox_new[4] + (0.2 * yrange) # ymax - top

bbox_new <- bbox_new %>%  # take the bounding box ...
  st_as_sfc() # ... and make it a sf polygon

m <- tm_shape(LST24Under55, unit = "mi" , bbox = bbox_new) + tm_fill(col = "red", title = "80th Percentile", legend.is.portrait = TRUE) +
  tm_shape(ne_states_sf_cb) + tm_polygons(alpha = 0.4, lwd = 0.2) +
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ne_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", xmod = 0.7, ymod = 0.2) +
  tm_scale_bar(breaks = c(0, 50, 100), text.size = 0.5, 
               position = c(0.4,0.005)) +
  tm_add_legend(type = "fill", 
                labels = "80th Percentile Block Groups",
                col = "red",
                border.col = "white", border.alpha = 0) +
  tm_layout(title = "Highest Exposure\nand Vulnerability\nto Day-Night Avg\nTemps for Children\nUnder Age 5\nJuly-Aug 2019", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"),
            legend.hist.width = 0.9)

tmap_save(m, "DATA/LST/ne_LST24Under55.png",
          height = 5, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/ne_LST24Under55.png")
```


```{r tabLST24hrHighRiskU5, fig.align = "center", fig.cap="Children Under 5 in Census Block Groups at or above 80th percentile of day-night average Land Surface Temperature."}
ne_blkgrpLST_sf %>%
  as.data.frame() %>%
  left_join(., statepopsdf,by="STATE") %>% 
  filter(LST24Under5Score == 5) %>%
  group_by(STATE) %>%
  summarize(`Under 5` = sum(over64E),
            `Pct of Under 5` = round(sum(under5E)/max(TotalUnder5)*100,1),
            Min = min(meanAvgLST,na.rm=TRUE),
            Mean = mean(meanAvgLST,na.rm=TRUE),
            Max = max(meanAvgLST,na.rm=TRUE)) %>% 
  mutate(`Pct of Under 5` = paste0(`Pct of Under 5`,"%"),
         Min = paste0(round(Min,2),"°"),
         Mean = paste0(round(Mean,2),"°"),
         Max = paste0(round(Max,2),"°")) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Children Under 5 in Highest Quintile of Day-Night Avg LST", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header")) %>% 
  kableExtra::column_spec(2:3, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 3, "Day-Night LST (°F)" = 3))
```


## Daytime Highest Exposure Populations of Concern {-}

```{r mapLSTDayHighRiskLI, fig.align = "center", fig.cap="Map of Block Groups in the highest quintile (80th percentile) of daytime average LST and 80th percentile of percentage of Low Income persons by Census Block Group, July-Aug 2019."}
LSTDayLowInc5 <- ne_blkgrpLST_sf %>%
  filter(LSTDayLowIncScore == 5)

bbox_new <- st_bbox(LSTDayLowInc5) # current bounding box

xrange <- bbox_new$xmax - bbox_new$xmin # range of x values
yrange <- bbox_new$ymax - bbox_new$ymin # range of y values

# bbox_new[1] <- bbox_new[1] - (0.5 * xrange) # xmin - left
bbox_new[3] <- bbox_new[3] + (0.2 * xrange) # xmax - right
# bbox_new[2] <- bbox_new[2] - (0.5 * yrange) # ymin - bottom
bbox_new[4] <- bbox_new[4] + (0.2 * yrange) # ymax - top

bbox_new <- bbox_new %>%  # take the bounding box ...
  st_as_sfc() # ... and make it a sf polygon

m <- tm_shape(LSTDayLowInc5, unit = "mi" , bbox = bbox_new) + tm_fill(col = "red", title = "80th Percentile", legend.is.portrait = TRUE) +
  tm_shape(ne_states_sf_cb) + tm_polygons(alpha = 0.4, lwd = 0.2) +
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ne_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", xmod = 0.7, ymod = 0.2) +
  tm_scale_bar(breaks = c(0, 50, 100), text.size = 0.5, 
               position = c(0.4,0.005)) +
  tm_add_legend(type = "fill", 
                labels = "80th Percentile Block Groups",
                col = "red",
                border.col = "white", border.alpha = 0) +
  tm_layout(title = "Highest Exposure\nand Vulnerability\nto Daytime Avg\nTemps for Low\nIncome Persons\nJuly-Aug 2019", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"),
            legend.hist.width = 0.9)

tmap_save(m, "DATA/LST/ne_LSTDayLowInc5.png",
          height = 5, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/ne_LSTDayLowInc5.png")
```


```{r tabLSTDayHighRiskLI, fig.align = "center", fig.cap="Low Income Persons in Census Block Groups at or above 80th percentile of daytime average Land Surface Temperature."}
ne_blkgrpLST_sf %>%
  as.data.frame() %>%
  left_join(., statepopsdf,by="STATE") %>% 
  filter(LSTDayLowIncScore == 5) %>%
  group_by(STATE) %>%
  summarize(`Low Income Persons` = sum(eng_limitE),
            `Pct of Low Income Persons` = round(sum(num2povE)/max(TotalLowInc)*100,1),
            Min = min(meanAvgLST,na.rm=TRUE),
            Mean = mean(meanAvgLST,na.rm=TRUE),
            Max = max(meanAvgLST,na.rm=TRUE)) %>% 
  mutate(`Pct of Low Income Persons` = paste0(`Pct of Low Income Persons`,"%"),
         Min = paste0(round(Min,2),"°"),
         Mean = paste0(round(Mean,2),"°"),
         Max = paste0(round(Max,2),"°")) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Low Income Persons in Highest Quintile of Daytime Avg LST", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header")) %>% 
  kableExtra::column_spec(2:3, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 3, "Daytime LST (°F)" = 3))
```


```{r mapLSTDayHighRiskNOHS, fig.align = "center", fig.cap="Map of Block Groups in the highest quintile (80th percentile) of daytime average LST and 80th percentile of percentage of adults with less than a high school diploma by Census Block Group, July-Aug 2019."}
LSTDayNoHSDip5 <- ne_blkgrpLST_sf %>%
  filter(LSTDayNoHSDipScore == 5)

bbox_new <- st_bbox(LSTDayNoHSDip5) # current bounding box

xrange <- bbox_new$xmax - bbox_new$xmin # range of x values
yrange <- bbox_new$ymax - bbox_new$ymin # range of y values

# bbox_new[1] <- bbox_new[1] - (0.5 * xrange) # xmin - left
bbox_new[3] <- bbox_new[3] + (0.2 * xrange) # xmax - right
# bbox_new[2] <- bbox_new[2] - (0.5 * yrange) # ymin - bottom
bbox_new[4] <- bbox_new[4] + (0.2 * yrange) # ymax - top

bbox_new <- bbox_new %>%  # take the bounding box ...
  st_as_sfc() # ... and make it a sf polygon

m <- tm_shape(LSTDayNoHSDip5, unit = "mi" , bbox = bbox_new) + tm_fill(col = "red", title = "80th Percentile", legend.is.portrait = TRUE) +
  tm_shape(ne_states_sf_cb) + tm_polygons(alpha = 0.4, lwd = 0.2) +
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ne_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", xmod = 0.7, ymod = 0.2) +
  tm_scale_bar(breaks = c(0, 50, 100), text.size = 0.5, 
               position = c(0.4,0.005)) +
  tm_add_legend(type = "fill", 
                labels = "80th Percentile Block Groups",
                col = "red",
                border.col = "white", border.alpha = 0) +
  tm_layout(title = "Highest Exposure\nand Vulnerability\nto Daytime Avg\nTemps for Adults\nwithout a High\nSchool Diploma\nJuly-Aug 2019", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"),
            legend.hist.width = 0.9)

tmap_save(m, "DATA/LST/ne_LSTDayNoHSDip5.png",
          height = 5, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/ne_LSTDayNoHSDip5.png")
```


```{r tabLSTDayHighRiskNOHS, fig.align = "center", fig.cap="Adults without a High School Diploma in Census Block Groups at or above 80th percentile of daytime average Land Surface Temperature."}
ne_blkgrpLST_sf %>%
  as.data.frame() %>%
  left_join(., statepopsdf,by="STATE") %>% 
  filter(LSTDayNoHSDipScore == 5) %>%
  group_by(STATE) %>%
  summarize(`No HS Dip` = sum(lthsE),
            `Pct of No HS Dip` = round(sum(lthsE)/max(TotalNoHS)*100,1),
            Min = min(meanAvgLST,na.rm=TRUE),
            Mean = mean(meanAvgLST,na.rm=TRUE),
            Max = max(meanAvgLST,na.rm=TRUE)) %>% 
  mutate(`Pct of No HS Dip` = paste0(`Pct of No HS Dip`,"%"),
         Min = paste0(round(Min,2),"°"),
         Mean = paste0(round(Mean,2),"°"),
         Max = paste0(round(Max,2),"°")) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Adults without a High School Diploma in Highest Quintile of Daytime Avg LST", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header")) %>% 
  kableExtra::column_spec(2:3, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 3, "Daytime LST (°F)" = 3))
```


```{r mapLSTDayHighRiskO64, fig.align = "center", fig.cap="Map of Block Groups in the highest quintile (80th percentile) of daytime average LST and 80th percentile of percentage of persons over 64 by Census Block Group, July-Aug 2019."}
LSTDayOver645 <- ne_blkgrpLST_sf %>%
  filter(LSTDayOver64Score == 5)

bbox_new <- st_bbox(LSTDayOver645) # current bounding box

xrange <- bbox_new$xmax - bbox_new$xmin # range of x values
yrange <- bbox_new$ymax - bbox_new$ymin # range of y values

# bbox_new[1] <- bbox_new[1] - (0.5 * xrange) # xmin - left
bbox_new[3] <- bbox_new[3] + (0.2 * xrange) # xmax - right
# bbox_new[2] <- bbox_new[2] - (0.5 * yrange) # ymin - bottom
bbox_new[4] <- bbox_new[4] + (0.2 * yrange) # ymax - top

bbox_new <- bbox_new %>%  # take the bounding box ...
  st_as_sfc() # ... and make it a sf polygon

m <- tm_shape(LSTDayOver645, unit = "mi" , bbox = bbox_new) + tm_fill(col = "red", title = "80th Percentile", legend.is.portrait = TRUE) +
  tm_shape(ne_states_sf_cb) + tm_polygons(alpha = 0.4, lwd = 0.2) +
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ne_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", xmod = 0.7, ymod = 0.2) +
  tm_scale_bar(breaks = c(0, 50, 100), text.size = 0.5, 
               position = c(0.4,0.005)) +
  tm_add_legend(type = "fill", 
                labels = "80th Percentile Block Groups",
                col = "red",
                border.col = "white", border.alpha = 0) +
  tm_layout(title = "Highest Exposure\nand Vulnerability\nto Daytime Avg\nTemps for Persons\nOver Age 64\nJuly-Aug 2019", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"),
            legend.hist.width = 0.9)

tmap_save(m, "DATA/LST/ne_LSTDayOver645.png",
          height = 5, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/ne_LSTDayOver645.png")
```


```{r tabLSTDayHighRiskO64, fig.align = "center", fig.cap="Persons over age 64 in Census Block Groups at or above 80th percentile of daytime average Land Surface Temperature."}
ne_blkgrpLST_sf %>%
  as.data.frame() %>%
  left_join(., statepopsdf,by="STATE") %>% 
  filter(LSTDayOver64Score == 5) %>%
  group_by(STATE) %>%
  summarize(`Over 64` = sum(over64E),
            `Pct of Over 64` = round(sum(over64E)/max(TotalOver64)*100,1),
            Min = min(meanAvgLST,na.rm=TRUE),
            Mean = mean(meanAvgLST,na.rm=TRUE),
            Max = max(meanAvgLST,na.rm=TRUE)) %>% 
  mutate(`Pct of Over 64` = paste0(`Pct of Over 64`,"%"),
         Min = paste0(round(Min,2),"°"),
         Mean = paste0(round(Mean,2),"°"),
         Max = paste0(round(Max,2),"°")) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Persons Over Age 64 in Highest Quintile of Daytime Avg LST", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header")) %>% 
  kableExtra::column_spec(2:3, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 3, "Daytime LST (°F)" = 3))
```


```{r mapLSTDayHighRiskU5, fig.align = "center", fig.cap="Map of Block Groups in the highest quintile (80th percentile) of daytime average LST and 80th percentile of percentage of children under 5 by Census Block Group, July-Aug 2019."}
LSTDayUnder55 <- ne_blkgrpLST_sf %>%
  filter(LSTDayUnder5Score == 5)

bbox_new <- st_bbox(LSTDayUnder55) # current bounding box

xrange <- bbox_new$xmax - bbox_new$xmin # range of x values
yrange <- bbox_new$ymax - bbox_new$ymin # range of y values

# bbox_new[1] <- bbox_new[1] - (0.5 * xrange) # xmin - left
bbox_new[3] <- bbox_new[3] + (0.2 * xrange) # xmax - right
# bbox_new[2] <- bbox_new[2] - (0.5 * yrange) # ymin - bottom
bbox_new[4] <- bbox_new[4] + (0.2 * yrange) # ymax - top

bbox_new <- bbox_new %>%  # take the bounding box ...
  st_as_sfc() # ... and make it a sf polygon

m <- tm_shape(LSTDayUnder55, unit = "mi" , bbox = bbox_new) + tm_fill(col = "red", title = "80th Percentile", legend.is.portrait = TRUE) +
  tm_shape(ne_states_sf_cb) + tm_polygons(alpha = 0.4, lwd = 0.2) +
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ne_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", xmod = 0.7, ymod = 0.2) +
  tm_scale_bar(breaks = c(0, 50, 100), text.size = 0.5, 
               position = c(0.4,0.005)) +
  tm_add_legend(type = "fill", 
                labels = "80th Percentile Block Groups",
                col = "red",
                border.col = "white", border.alpha = 0) +
  tm_layout(title = "Highest Exposure\nand Vulnerability\nto Daytime Avg\nTemps for Children\nUnder Age 5\nJuly-Aug 2019", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"),
            legend.hist.width = 0.9)

tmap_save(m, "DATA/LST/ne_LSTDayUnder55.png",
          height = 5, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/ne_LSTDayUnder55.png")
```


```{r tabLSTDayHighRiskU5, fig.align = "center", fig.cap="Children Under 5 in Census Block Groups at or above 80th percentile of daytime average Land Surface Temperature."}
ne_blkgrpLST_sf %>%
  as.data.frame() %>%
  left_join(., statepopsdf,by="STATE") %>% 
  filter(LSTDayUnder5Score == 5) %>%
  group_by(STATE) %>%
  summarize(`Under 5` = sum(under5E),
            `Pct of Under 5` = round(sum(under5E)/max(TotalUnder5)*100,1),
            Min = min(meanAvgLST,na.rm=TRUE),
            Mean = mean(meanAvgLST,na.rm=TRUE),
            Max = max(meanAvgLST,na.rm=TRUE)) %>% 
  mutate(`Pct of Under 5` = paste0(`Pct of Under 5`,"%"),
         Min = paste0(round(Min,2),"°"),
         Mean = paste0(round(Mean,2),"°"),
         Max = paste0(round(Max,2),"°")) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Children Under 5 in Highest Quintile of Daytime Avg LST", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header")) %>% 
  kableExtra::column_spec(2:3, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 3, "Daytime LST (°F)" = 3))
```


## Nighttime Highest Exposure Populations of Concern {-}

```{r mapLSTNightHighRiskLI, fig.align = "center", fig.cap="Map of Block Groups in the highest quintile (80th percentile) of nighttime average LST and 80th percentile of percentage of Low Income persons by Census Block Group, July-Aug 2019."}
LSTNightLowInc5 <- ne_blkgrpLST_sf %>%
  filter(LSTNightLowIncScore == 5)

bbox_new <- st_bbox(LSTNightLowInc5) # current bounding box

xrange <- bbox_new$xmax - bbox_new$xmin # range of x values
yrange <- bbox_new$ymax - bbox_new$ymin # range of y values

# bbox_new[1] <- bbox_new[1] - (0.5 * xrange) # xmin - left
bbox_new[3] <- bbox_new[3] + (0.2 * xrange) # xmax - right
# bbox_new[2] <- bbox_new[2] - (0.5 * yrange) # ymin - bottom
bbox_new[4] <- bbox_new[4] + (0.2 * yrange) # ymax - top

bbox_new <- bbox_new %>%  # take the bounding box ...
  st_as_sfc() # ... and make it a sf polygon

m <- tm_shape(LSTNightLowInc5, unit = "mi" , bbox = bbox_new) + tm_fill(col = "red", title = "80th Percentile", legend.is.portrait = TRUE) +
  tm_shape(ne_states_sf_cb) + tm_polygons(alpha = 0.4, lwd = 0.2) +
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ne_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", xmod = 0.7, ymod = 0.2) +
  tm_scale_bar(breaks = c(0, 50, 100), text.size = 0.5, 
               position = c(0.4,0.005)) +
  tm_add_legend(type = "fill", 
                labels = "80th Percentile Block Groups",
                col = "red",
                border.col = "white", border.alpha = 0) +
  tm_layout(title = "Highest Exposure\nand Vulnerability\nto Nighttime Avg\nTemps for Low\nIncome Persons\nJuly-Aug 2019", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"),
            legend.hist.width = 0.9)

tmap_save(m, "DATA/LST/ne_LSTNightLowInc5.png",
          height = 5, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/ne_LSTNightLowInc5.png")
```


```{r tabLSTNightHighRiskLI, fig.align = "center", fig.cap="Low Income Persons in Census Block Groups at or above 80th percentile of nighttime average Land Surface Temperature."}
ne_blkgrpLST_sf %>%
  as.data.frame() %>%
  left_join(., statepopsdf,by="STATE") %>% 
  filter(LSTNightLowIncScore == 5) %>%
  group_by(STATE) %>%
  summarize(`Low Income Persons` = sum(eng_limitE),
            `Pct of Low Income Persons` = round(sum(num2povE)/max(TotalLowInc)*100,1),
            Min = min(meanAvgLST,na.rm=TRUE),
            Mean = mean(meanAvgLST,na.rm=TRUE),
            Max = max(meanAvgLST,na.rm=TRUE)) %>% 
  mutate(`Pct of Low Income Persons` = paste0(`Pct of Low Income Persons`,"%"),
         Min = paste0(round(Min,2),"°"),
         Mean = paste0(round(Mean,2),"°"),
         Max = paste0(round(Max,2),"°")) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Low Income Persons in Highest Quintile of Nighttime Avg LST", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header")) %>% 
  kableExtra::column_spec(2:3, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 3, "Nighttime LST (°F)" = 3))
```


```{r mapLSTNightHighRiskNOHS, fig.align = "center", fig.cap="Map of Block Groups in the highest quintile (80th percentile) of nighttime average LST and 80th percentile of percentage of adults with less than a high school diploma by Census Block Group, July-Aug 2019."}
LSTNightNoHSDip5 <- ne_blkgrpLST_sf %>%
  filter(LSTNightNoHSDipScore == 5)

bbox_new <- st_bbox(LSTNightNoHSDip5) # current bounding box

xrange <- bbox_new$xmax - bbox_new$xmin # range of x values
yrange <- bbox_new$ymax - bbox_new$ymin # range of y values

# bbox_new[1] <- bbox_new[1] - (0.5 * xrange) # xmin - left
bbox_new[3] <- bbox_new[3] + (0.2 * xrange) # xmax - right
# bbox_new[2] <- bbox_new[2] - (0.5 * yrange) # ymin - bottom
bbox_new[4] <- bbox_new[4] + (0.2 * yrange) # ymax - top

bbox_new <- bbox_new %>%  # take the bounding box ...
  st_as_sfc() # ... and make it a sf polygon

m <- tm_shape(LSTNightNoHSDip5, unit = "mi" , bbox = bbox_new) + tm_fill(col = "red", title = "80th Percentile", legend.is.portrait = TRUE) +
  tm_shape(ne_states_sf_cb) + tm_polygons(alpha = 0.4, lwd = 0.2) +
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ne_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", xmod = 0.7, ymod = 0.2) +
  tm_scale_bar(breaks = c(0, 50, 100), text.size = 0.5, 
               position = c(0.4,0.005)) +
  tm_add_legend(type = "fill", 
                labels = "80th Percentile Block Groups",
                col = "red",
                border.col = "white", border.alpha = 0) +
  tm_layout(title = "Highest Exposure\nand Vulnerability\nto Nighttime Avg\nTemps for Adults\nwithout a High\nSchool Diploma\nJuly-Aug 2019", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"),
            legend.hist.width = 0.9)

tmap_save(m, "DATA/LST/ne_LSTNightNoHSDip5.png",
          height = 5, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/ne_LSTNightNoHSDip5.png")
```


```{r tabLSTNightHighRiskNOHS, fig.align = "center", fig.cap="Adults without a High School Diploma in Census Block Groups at or above 80th percentile of nighttime average Land Surface Temperature."}
ne_blkgrpLST_sf %>%
  as.data.frame() %>%
  left_join(., statepopsdf,by="STATE") %>% 
  filter(LSTNightNoHSDipScore == 5) %>%
  group_by(STATE) %>%
  summarize(`No HS Dip` = sum(lthsE),
            `Pct of No HS Dip` = round(sum(lthsE)/max(TotalNoHS)*100,1),
            Min = min(meanAvgLST,na.rm=TRUE),
            Mean = mean(meanAvgLST,na.rm=TRUE),
            Max = max(meanAvgLST,na.rm=TRUE)) %>% 
  mutate(`Pct of No HS Dip` = paste0(`Pct of No HS Dip`,"%"),
         Min = paste0(round(Min,2),"°"),
         Mean = paste0(round(Mean,2),"°"),
         Max = paste0(round(Max,2),"°")) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Adults without a High School Diploma in Highest Quintile of Nighttime Avg LST", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header")) %>% 
  kableExtra::column_spec(2:3, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 3, "Nighttime LST (°F)" = 3))
```


```{r mapLSTNightHighRiskO64, fig.align = "center", fig.cap="Map of Block Groups in the highest quintile (80th percentile) of nighttime average LST and 80th percentile of percentage of persons over 64 by Census Block Group, July-Aug 2019."}
LSTNightOver645 <- ne_blkgrpLST_sf %>%
  filter(LSTNightOver64Score == 5)

bbox_new <- st_bbox(LSTNightOver645) # current bounding box

xrange <- bbox_new$xmax - bbox_new$xmin # range of x values
yrange <- bbox_new$ymax - bbox_new$ymin # range of y values

# bbox_new[1] <- bbox_new[1] - (0.5 * xrange) # xmin - left
bbox_new[3] <- bbox_new[3] + (0.2 * xrange) # xmax - right
# bbox_new[2] <- bbox_new[2] - (0.5 * yrange) # ymin - bottom
bbox_new[4] <- bbox_new[4] + (0.2 * yrange) # ymax - top

bbox_new <- bbox_new %>%  # take the bounding box ...
  st_as_sfc() # ... and make it a sf polygon

m <- tm_shape(LSTNightOver645, unit = "mi" , bbox = bbox_new) + tm_fill(col = "red", title = "80th Percentile", legend.is.portrait = TRUE) +
  tm_shape(ne_states_sf_cb) + tm_polygons(alpha = 0.4, lwd = 0.2) +
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ne_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", xmod = 0.7, ymod = 0.2) +
  tm_scale_bar(breaks = c(0, 50, 100), text.size = 0.5, 
               position = c(0.4,0.005)) +
  tm_add_legend(type = "fill", 
                labels = "80th Percentile Block Groups",
                col = "red",
                border.col = "white", border.alpha = 0) +
  tm_layout(title = "Highest Exposure\nand Vulnerability\nto Nighttime Avg\nTemps for Persons\nOver Age 64\nJuly-Aug 2019", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"),
            legend.hist.width = 0.9)

tmap_save(m, "DATA/LST/ne_LSTNightOver645.png",
          height = 5, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/ne_LSTNightOver645.png")
```


```{r tabLSTNightHighRiskO64, fig.align = "center", fig.cap="Persons over age 64 in Census Block Groups at or above 80th percentile of nighttime average Land Surface Temperature."}
ne_blkgrpLST_sf %>%
  as.data.frame() %>%
  left_join(., statepopsdf,by="STATE") %>% 
  filter(LSTNightOver64Score == 5) %>%
  group_by(STATE) %>%
  summarize(`Over 64` = sum(over64E),
            `Pct of Over 64` = round(sum(over64E)/max(TotalOver64)*100,1),
            Min = min(meanAvgLST,na.rm=TRUE),
            Mean = mean(meanAvgLST,na.rm=TRUE),
            Max = max(meanAvgLST,na.rm=TRUE)) %>% 
  mutate(`Pct of Over 64` = paste0(`Pct of Over 64`,"%"),
         Min = paste0(round(Min,2),"°"),
         Mean = paste0(round(Mean,2),"°"),
         Max = paste0(round(Max,2),"°")) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Persons Over Age 64 in Highest Quintile of Nighttime Avg LST", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header")) %>% 
  kableExtra::column_spec(2:3, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 3, "Nighttime LST (°F)" = 3))
```


```{r mapLSTNightHighRiskU5, fig.align = "center", fig.cap="Map of Block Groups in the highest quintile (80th percentile) of nighttime average LST and 80th percentile of percentage of children under 5 by Census Block Group, July-Aug 2019."}
LSTNightUnder55 <- ne_blkgrpLST_sf %>%
  filter(LSTNightUnder5Score == 5)

bbox_new <- st_bbox(LSTNightUnder55) # current bounding box

xrange <- bbox_new$xmax - bbox_new$xmin # range of x values
yrange <- bbox_new$ymax - bbox_new$ymin # range of y values

# bbox_new[1] <- bbox_new[1] - (0.5 * xrange) # xmin - left
bbox_new[3] <- bbox_new[3] + (0.2 * xrange) # xmax - right
# bbox_new[2] <- bbox_new[2] - (0.5 * yrange) # ymin - bottom
bbox_new[4] <- bbox_new[4] + (0.2 * yrange) # ymax - top

bbox_new <- bbox_new %>%  # take the bounding box ...
  st_as_sfc() # ... and make it a sf polygon

m <- tm_shape(LSTNightUnder55, unit = "mi" , bbox = bbox_new) + tm_fill(col = "red", title = "80th Percentile", legend.is.portrait = TRUE) +
  tm_shape(ne_states_sf_cb) + tm_polygons(alpha = 0.4, lwd = 0.2) +
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ne_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", xmod = 0.7, ymod = 0.2) +
  tm_scale_bar(breaks = c(0, 50, 100), text.size = 0.5, 
               position = c(0.4,0.005)) +
  tm_add_legend(type = "fill", 
                labels = "80th Percentile Block Groups",
                col = "red",
                border.col = "white", border.alpha = 0) +
  tm_layout(title = "Highest Exposure\nand Vulnerability\nto Nighttime Avg\nTemps for Children\nUnder Age 5\nJuly-Aug 2019", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"),
            legend.hist.width = 0.9)

tmap_save(m, "DATA/LST/ne_LSTNightUnder55.png",
          height = 5, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/ne_LSTNightUnder55.png")
```


```{r tabLSTNightHighRiskU5, fig.align = "center", fig.cap="Children Under 5 in Census Block Groups at or above 80th percentile of nighttime average Land Surface Temperature."}
ne_blkgrpLST_sf %>%
  as.data.frame() %>%
  left_join(., statepopsdf,by="STATE") %>% 
  filter(LSTNightUnder5Score == 5) %>%
  group_by(STATE) %>%
  summarize(`Under 5` = sum(under5E),
            `Pct of Under 5` = round(sum(under5E)/max(TotalUnder5)*100,1),
            Min = min(meanAvgLST,na.rm=TRUE),
            Mean = mean(meanAvgLST,na.rm=TRUE),
            Max = max(meanAvgLST,na.rm=TRUE)) %>% 
  mutate(`Pct of Under 5` = paste0(`Pct of Under 5`,"%"),
         Min = paste0(round(Min,2),"°"),
         Mean = paste0(round(Mean,2),"°"),
         Max = paste0(round(Max,2),"°")) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Children Under 5 in Highest Quintile of Nighttime Avg LST", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header")) %>% 
  kableExtra::column_spec(2:3, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 3, "Nighttime LST (°F)" = 3))
```



## Correlations {-}

```{r cormatrixLSTpopsNE, echo=FALSE, message=FALSE, warning=FALSE, fig.align = "center", fig.cap="Spearman's correlation matrix of temperature and the proportions of populations of concern by Census Block Group."}
# Create a correlation matrix of PM25 and populations of concern
corrLSTpops <- ne_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  dplyr::transmute(`Day-Night Avg LST` = meanAvgLST,
                   `Day LST` = meanDayLST,
                   `Night LST` = meanNightLST, 
                Minority = minority_pctE, 
                `Low Income` = pct2povE, 
                `Lang Isol` = eng_limit_pctE, 
                `No HS Dip` = pct_lthsE, 
                `Under 5` = pct_under5E, 
                `Over 64` = pct_over64E) %>% 
  drop_na() %>% 
  cor(method = "spearman")
# corrPM25
ggcorrplot(corrLSTpops, hc.order = TRUE, type = "lower", lab = TRUE, 
           title = "Correlation Matrix of LST \nand Populations for New England", 
           legend.title = "Spearman's \nCorrelation\nCoefficient",
           lab_size = 2)
```



```{r cormatrixLSTAvgNEpollutants, echo=FALSE, message=FALSE, warning=FALSE, fig.align = "center", fig.cap="Spearman's correlation matrix of temperature and pollution by Census Block Group."}
# Create a correlation matrix of PM25 and populations of concern
corrLSTemissions <- ne_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  dplyr::transmute(`Day-Night Avg LST` = meanAvgLST,
                   `Day LST` = meanDayLST,
                   `Night LST` = meanNightLST,
                `Diesel PM` = DSLPM_19, 
                Ozone = OZONE_19, 
                PM2.5 = PM25_19,
                `Cancer Risk` = CANCER_19,
                `Resp Risk` = RESP_19) %>% 
  drop_na() %>% 
  cor(method = "spearman")
# corrPM25
ggcorrplot(corrLSTemissions, hc.order = TRUE, type = "lower", lab = TRUE, 
           title = "Correlation Matrix of LST \nand Air Pollution for New England",
           legend.title = "Spearman's \nCorrelation\nCoefficient", 
           lab_size = 2)
```

```{r LSTemissionsGWR}
# library(spgwr)
# 
# # convert to sp
# ne_blkgrpLST_sp <- as_Spatial(ne_blkgrpLST_sf)
# 
# #calculate kernel bandwidth
# GWRbandwidth <- gwr.sel(ne_blkgrpLST_sp$PM25_19 ~ ne_blkgrpLST_sp$meanAvgLST, data=ne_blkgrpLST_sp, adapt=T)


```