LST_VermontPDF.Rmd

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

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo=FALSE, message=FALSE, warning=FALSE)
```
\pagebreak
# Heat risk in Vermont

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 Vermont experiences a generally cool climate, research has shown that people in this state are actually more sensitive to elevated temperatures than those living in warmer states 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 Vermont 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 state each year.[^4thNatlAssessment] The combination of rapidly increasing temperatures and higher sensitivity puts people in Vermont 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 state 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 state. That week constituted the conclusion of a historically warm July for Vermont 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.
# vt_bmap <- ne_states_sf_cb %>% 
#   st_transform(., crs = 4326) %>% 
#   st_bbox(.) %>% 
#   read_osm(., type = "stamen-watercolor")
# 
# # map basemap
# tm_shape(vt_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 Vermont
# # isolate data for Vermont
# vt_blkgrp_sf <- ne_blkgrp_sf %>%
#   filter(STATE == "Vermont") %>%
#   st_transform(., crs = 2852)
# 
# vt_state_sf_cb <- ne_states_sf_cb %>%
#   filter(NAME == "Vermont") %>%
#   st_transform(., crs = 2852)
# 
# # define new projection
# newproj <- st_crs(vt_blkgrp_sf)[[2]]
# 
# LSTavgJuly2019r_ma <- LSTavgJuly2019r %>%
#   projectRaster(., crs = newproj) %>%
#   crop(., vt_state_sf_cb) %>%
#   mask(., vt_state_sf_cb)
# 
# LSTdayJuly2019r_ma <- LSTdayJuly2019r %>%
#   projectRaster(., crs = newproj) %>%
#   crop(., vt_state_sf_cb) %>%
#   mask(., vt_state_sf_cb)
# 
# LSTnightJuly2019r_ma <- LSTnightJuly2019r %>%
#   projectRaster(., crs = newproj) %>%
#   crop(., vt_state_sf_cb) %>%
#   mask(., vt_state_sf_cb)
# 
# LSTDayNighDiffJuly2019r_ma <- LSTDayNighDiffJuly2019r %>%
#   projectRaster(., crs = newproj) %>%
#   crop(., vt_state_sf_cb) %>%
#   mask(., vt_state_sf_cb)
# 
# # download Census-designated urban areas to contrast with rural areas
# us_urban_sf <- urban_areas(cb = TRUE)
# 
# # create urban area polygon for Vermont
# vt_urban_sf <- us_urban_sf %>%
#   st_centroid(.) %>%
#   st_transform(., crs = 2852) %>%
#   st_join(., vt_state_sf_cb, join = st_within, left = FALSE) %>%
#   as.data.frame() %>%
#   dplyr::select(GEOID10) %>%
#   inner_join(us_urban_sf,.,by="GEOID10") %>%
#   st_transform(., crs = 2852) %>%
#   group_by() %>%
#   summarize()
# 
# # create rural area polygon
# vt_rural_sf <- vt_state_sf_cb %>%
#   group_by() %>%
#   summarize() %>%
#   st_difference(., vt_urban_sf) %>%
#   st_make_valid()
# 
# # # select blkgrps that do not intersect wtih urban polygons - rural blkgrps
# # vt_blkgrpLST_rural <- st_disjoint(vt_blkgrp_sf,
# #                                   vt_urban_sf,sparse = FALSE) %>%
# #   apply(., 1, all) %>% # convert dense matrix to logical vector
# #   vt_blkgrpLST_sf[.,] # subset blkgrps where there is NO intersection
# #
# # # select blkgrps that intersect with urban polygons - urban blkgrps
# # vt_blkgrpLST_urban <- st_intersects(vt_blkgrp_sf,
# #                                     vt_urban_sf,sparse = FALSE) %>%
# #   apply(., 1, any) %>% # convert dense matrix to logical vector
# #   vt_blkgrpLST_sf[.,] # subset blkgrps where there is ANY intersection
# #
# # vt_blkgrp_urban_simple <- vt_blkgrpLST_urban %>%
# #   group_by() %>%
# #   summarize()
# #
# # vt_blkgrp_rural_simple <- vt_blkgrpLST_rural %>%
# #   group_by() %>%
# #   summarize()
# 
# # crop LST to urban
# # LSTavgJuly2019r_maUrban <- LSTavgJuly2019r %>%
# #   projectRaster(., crs = newproj) %>%
# #   crop(., vt_urban_sf) %>%
# #   mask(., vt_urban_sf)
# #
# # LSTdayJuly2019r_maUrban <- LSTdayJuly2019r %>%
# #   projectRaster(., crs = newproj) %>%
# #   crop(., vt_urban_sf) %>%
# #   mask(., vt_urban_sf)
# #
# # LSTnightJuly2019r_maUrban <- LSTnightJuly2019r %>%
# #   projectRaster(., crs = newproj) %>%
# #   crop(., vt_urban_sf) %>%
# #   mask(., vt_urban_sf)
# #
# # LSTDayNighDiffJuly2019r_maUrban <- LSTDayNighDiffJuly2019r %>%
# #   projectRaster(., crs = newproj) %>%
# #   crop(., vt_urban_sf) %>%
# #   mask(., vt_urban_sf)
# 
# # crop LST to rural
# LSTavgJuly2019r_maRural <- LSTavgJuly2019r %>%
#   projectRaster(., crs = newproj) %>%
#   crop(., vt_rural_sf) %>%
#   mask(., vt_rural_sf)
# 
# LSTdayJuly2019r_maRural <- LSTdayJuly2019r %>%
#   projectRaster(., crs = newproj) %>%
#   crop(., vt_rural_sf) %>%
#   mask(., vt_rural_sf)
# 
# LSTnightJuly2019r_maRural <- LSTnightJuly2019r %>%
#   projectRaster(., crs = newproj) %>%
#   crop(., vt_rural_sf) %>%
#   mask(., vt_rural_sf)
# 
# LSTDayNighDiffJuly2019r_maRural <- LSTDayNighDiffJuly2019r %>%
#   projectRaster(., crs = newproj) %>%
#   crop(., vt_rural_sf) %>%
#   mask(., vt_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_maRural, stat = "mean")
# # LSTdayJuly2019mean <- cellStats(x = LSTdayJuly2019r_ne, stat = "mean")
# # LSTdayJuly2019Urbanmean <- cellStats(x = LSTdayJuly2019r_neUrban, stat = "mean")
# LSTdayJuly2019Ruralmean <- cellStats(x = LSTdayJuly2019r_maRural, stat = "mean")
# # LSTnightJuly2019mean <- cellStats(x = LSTnightJuly2019r_ne, stat = "mean")
# # LSTnightJuly2019Urbanmean <- cellStats(x = LSTnightJuly2019r_neUrban, stat = "mean")
# LSTnightJuly2019Ruralmean <- cellStats(x = LSTnightJuly2019r_maRural, 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(vt_blkgrp_sf)))
# # identiy empty geometries
# empty_geo <- st_is_empty(vt_blkgrp_sf)
# # filter out empty geometries
# vt_blkgrp_sf <- vt_blkgrp_sf[!empty_geo,]
# # clean up
# rm(empty_geo)
# 
# # Extract mean LST values within each block group to a df
# meanLSTavg_df <- vt_blkgrp_sf %>%
#   dplyr::select(GEOID) %>%
#   as_Spatial() %>%
#   extract(LSTavgJuly2019r_ma, .,
#           fun=mean, sp=TRUE, na.rm=TRUE, small=TRUE) %>%
#   as.data.frame() %>%
#   rename(meanAvgLST = layer)
# 
# meanLSTday_df <- vt_blkgrp_sf %>%
#   dplyr::select(GEOID) %>%
#   as_Spatial() %>%
#   extract(LSTdayJuly2019r_ma, .,
#           fun=mean, sp=TRUE, na.rm=TRUE, small=TRUE) %>%
#   as.data.frame() %>%
#   rename(meanDayLST = layer)
# 
# meanLSTnight_df <- vt_blkgrp_sf %>%
#   dplyr::select(GEOID) %>%
#   as_Spatial() %>%
#   extract(LSTnightJuly2019r_ma, .,
#           fun=mean, sp=TRUE, na.rm=TRUE, small=TRUE) %>%
#   as.data.frame() %>%
#   rename(meanNightLST = layer)
# 
# meanLSTDayNightDiff_df <- vt_blkgrp_sf %>%
#   dplyr::select(GEOID) %>%
#   as_Spatial() %>%
#   extract(LSTDayNighDiffJuly2019r_ma,.,
#           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 <- vt_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
# vt_blkgrpLST_sf <- vt_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(vt_blkgrpLST_sf, vt_urban_sf, vt_rural_sf, file = "DATA/LST/vt_blkgrpLST_sf.rds")
# 
# # clean up rasters
# rm(list = ls(pattern = "July2019r"))

# load previously processed LST data
load("DATA/LST/vt_blkgrpLST_sf.rds")


vt_towns_sf_pts <- county_subdivisions(state = "VT", cb = TRUE) %>% 
  filter(NAME %in% c("Brattleboro",
                     "Manchester",
                     "Rutland",
                     "Middlebury",
                     "Montpelier",
                     "Burlington",
                     "Stowe",
                     "Greensboro",
                     "Addison",
                     "Swanton",
                     "Berlin",
                     "Bennington",
                     "St. Albans"),
         !GEOID %in% c("5002161225","5001161675")) %>% # remove duplicate Rutland and St. Albans points 
  st_transform(., crs = 2852) %>% 
  st_centroid(of_largest_polygon = TRUE)

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

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

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

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

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

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

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

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

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

# Create custom icons of highway shields
I89 <- tmap_icons(file = "https://upload.wikimedia.org/wikipedia/commons/thumb/f/f4/I-89.svg/200px-I-89.svg.png")
I91 <- tmap_icons("https://upload.wikimedia.org/wikipedia/commons/thumb/9/90/I-91.svg/200px-I-91.svg.png")
I93 <- tmap_icons("https://upload.wikimedia.org/wikipedia/commons/thumb/0/0d/I-93.svg/200px-I-93.svg.png")
Hwy2 <- tmap_icons("https://upload.wikimedia.org/wikipedia/commons/thumb/2/27/US_2.svg/200px-US_2.svg.png")
Hwy4 <- tmap_icons("https://upload.wikimedia.org/wikipedia/commons/thumb/7/71/US_4.svg/200px-US_4.svg.png")
Hwy7 <- tmap_icons("https://upload.wikimedia.org/wikipedia/commons/thumb/9/95/US_7.svg/200px-US_7.svg.png")

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

# Calculate percentiles for pops of concern and temperature and create ranks for comparison
vt_blkgrpLST_sf <- vt_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 priority populations for use in computing percentages in later tables
statepopsdf <- vt_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  group_by(STATE) %>% 
  summarize(TotalPop = sum(totalpopE, na.rm = TRUE),
            TotalHH = sum(householdsE, na.rm = TRUE),
            TotalLEH = sum(eng_limitE, na.rm = TRUE),
            TotalMin = sum(minorityE, na.rm = TRUE),
            TotalLowInc = sum(num2povE, na.rm = TRUE),
            TotalNoHS = sum(lthsE, na.rm = TRUE),
            TotalOver64 = sum(over64E, na.rm = TRUE),
            TotalUnder5 = sum(under5E, na.rm = TRUE))

# vt_towns_sf <- ne_towns_sf %>% 
#   dplyr::select(GEOID,NAME) %>% 
#   st_transform(., crs = 2852)

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

townpopsdf <- vt_blkgrpLST_sf %>% 
  dplyr::select(-GEOID) %>% 
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  group_by(NAME) %>% 
  summarize(GEOID = first(GEOID),
            TotalPop = sum(totalpopE, na.rm = TRUE),
            TotalHH = sum(householdsE, na.rm = TRUE),
            TotalLEH = sum(eng_limitE, na.rm = TRUE),
            TotalMin = sum(minorityE, na.rm = TRUE),
            TotalLowInc = sum(num2povE, na.rm = TRUE),
            TotalNoHS = sum(lthsE, na.rm = TRUE),
            TotalOver64 = sum(over64E, na.rm = TRUE),
            TotalUnder5 = sum(under5E, na.rm = TRUE))

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


## 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. These day-night average temperatures varied significantly across Vermont, with highest temperatures apparent along the shore of Lake Champlain and in the southeast along the I-91 corridor (see Figure \@ref(fig:mapLSTavgVT)). 

```{r mapLSTavgVT, fig.align = "center", fig.cap="Map of July-Aug 2019 Day-Night Average Land Surface Temperatures (LST) across Vermont at Census Block Group level.", strip.white=TRUE}
# Map of average LST across Vermont
m <- tm_shape(vt_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(lwd = 0.2, alpha = 0.8) + 
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ny_state_sf) + tm_borders(lwd = 0.2, alpha = 0.8) +
  tm_shape(vt_highways) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(vt_highways2nd) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.2) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.2) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.2) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.2) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), 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/vt_24hrLST.png",
          height = 8, width = 8, units = "in", dpi = 600)

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

```{r LST24outliers, include=FALSE}
# Identify municipalities with high heat outliers in Vermont
vt_LST24outliers <- vt_blkgrpLST_sf %>% 
  arrange(desc(meanAvgLST)) %>% 
  slice(1:10) %>% 
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  group_by(NAME) %>%
  summarize(unique(NAME)) %>% 
  pull(NAME) %>% 
  paste(.,sep = " ", collapse = ", ")
```

Figure \@ref(fig:boxplotLSTavgVT) is a boxplot of average day-night temperatures by Block Group . The box represents temperatures ranging between the 25th and 75th percentiles. The line that divides the box into 2 parts represents the median temperature for all Block Groups, which in this case is `r round(median(vt_blkgrpLST_sf$meanAvgLST, na.rm = TRUE),2)`°. Half of the state's Block Groups are below the median and half are above the median. Each dot represents an individual Block Group. The top ten highest outlier Block Groups (far right) were found in the following municipalities: `r vt_LST24outliers`. 

```{r boxplotLSTavgVT, fig.align = "center", fig.cap="Boxplot of July-Aug 2019 day-night average Land Surface Temperatures (LST) across Vermont at Census Block Group level."}
# boxplot of PM2.5 by state for 2016
vt_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  drop_na(meanAvgLST) %>% 
  ggplot(aes(x = STATE, y = meanAvgLST, fill = STATE)) + 
  geom_boxplot() + 
  geom_jitter(color = "black", size = 0.4, alpha = 0.6) +
  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 Vermont.[^urban] Specifically, the average temperature of the rural areas of Vermont is substracted 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 Vermont are most pronounced in the inner core of urbanized areas, (see Figure \@ref(fig:mapUHIavgVT)), particularly Burlington. Other areas with elevated temperatures include Addison, St. Albans and Swanton. Temperatures in these areas are up to `r round(max(vt_blkgrpLST_sf$UHI24avg, na.rm = TRUE),1)`° warmer than the rural background average. By contrast, most of the central part of the state is closer to the rural average. 

```{r mapUHIavgVT, fig.align = "center", fig.cap="Map of July-Aug 2019 day-night average Urban Heat Island (UHI) effect across Vermont 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 Vermont
m <- tm_shape(vt_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(lwd = 0.2, alpha = 0.8) + 
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ny_state_sf) + tm_borders(lwd = 0.2, alpha = 0.8) +
  tm_shape(vt_highways) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(vt_highways2nd) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.15) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), 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/vt_24hrUHI.png",
          height = 8, width = 8, units = "in", dpi = 600)

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

```{r LST24UHIoutliers, include=FALSE}
vt_LST24UHIoutliers <- vt_blkgrpLST_sf %>% 
  arrange(desc(UHI24avg)) %>% 
  slice(1:10) %>% 
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  group_by(NAME) %>%
  summarize(unique(NAME)) %>% 
  pull(NAME) %>% 
  paste(.,sep = " ", collapse = ", ")
```

Figure \@ref(fig:boxplotUHIavgVT) is a boxplot of day-night average Urban Heat Island (UHI) effect by Block Group . The box represents UHI values ranging between the 25th and 75th percentiles. The line that divides the box into 2 parts represents the median UHI for all Block Groups, which in this case is `r round(median(vt_blkgrpLST_sf$UHI24avg, na.rm = TRUE),2)`. Half of the state's Block Groups are below the median and half are above the median. Each dot represents an individual Block Group. The top ten highest outlier Block Groups (far right) were found in the following municipalities: `r vt_LST24UHIoutliers`. 

```{r boxplotUHIavgVT, fig.align = "center", fig.cap="Boxplot of July-Aug 2019 day-night average Urban Heat Island (UHI) across Vermont 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
vt_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  drop_na(UHI24avg) %>% 
  ggplot(aes(x = STATE, y = UHI24avg, fill = STATE)) + 
  geom_boxplot() + 
  geom_jitter(color = "black", size = 0.4, alpha = 0.6) + 
  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 Priority Populations

```{r LST24hrTempPop, include=FALSE}
# Pop Weighted avg of 24hr LST exposure for all Groups in Vermont relative to VT average
LST24hrPop <- vt_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  # filter(STATE == "Vermont") %>%
  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 Vermont relative to VT average
LST24hrPopTab <- vt_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  # filter(STATE == "Vermont") %>%
  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:plotLST24hrTempPopAvgVT) and Table \@ref(tab:tabLST24hrPopAvgVT) show population-weighted exposures for priority populations relative to the day-night average LST for the state of `r round(mean(vt_blkgrpLST_sf$meanAvgLST,na.rm=TRUE),1)`°. For example, limited English speaking households in Vermont 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 state day-night average LST. Similarly, People of Color 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 state 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 state average.

```{r plotLST24hrTempPopAvgVT, echo=FALSE, message=FALSE, warning=FALSE, fig.align = "center", fig.cap="Population-weighted average exposures to day-night average temperatures for priority populations in Vermont relative to the state average."}
LST24hrPop %>% 
  mutate(Group = recode(Group, "Minority" = "People of Color",
                        "No HS Dip" = "No HS Diploma")) %>%
  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 = 71.5, label = "Above state\naverage"),
            color = "gray48") +
  geom_text(aes(x = "Low Income", y = 69, label = "Below state\naverage"),
            color = "gray48") +
  # geom_segment(aes(x = "Under 5", xend = "Under 5", y = 81, yend = 82), 
  #              arrow = arrow(length = unit(0.3,"cm"))) +
  expand_limits(y = c(68.5,72.5))

ggsave("images/VT_LST24_graph.png")
```


```{r tabLST24hrPopAvgVT, fig.align = "center", fig.cap="Population-weighted average exposures to day-night average LST for priority populations in Vermont relative to the state average."}
LST24hrPopTab %>% 
  mutate(Group = recode(Group, "Minority" = "People of Color",
                        "No HS Dip" = "No HS Diploma")) %>%
  arrange(Group) %>% 
  kableExtra::kable(longtable = T, booktabs = T, digits = 2, 
                    col.names = c("Group","Avg 24hr LST (°F)",
                                  "Difference from State Avg* (°F)",
                                  "Pct Above/Below State 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("State average day-night LST is ", round(mean(vt_blkgrpLST_sf$meanAvgLST,na.rm=TRUE),2), "°"))
```

```{r UHI24hrTempPop, include=FALSE}
# Pop Weighted avg of 24hr LST exposure for all Groups in Vermont relative to VT average
UHI24hrPop <- vt_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  # filter(STATE == "Vermont") %>%
  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 Vermont relative to VT average
UHI24hrPopTab <- vt_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  # filter(STATE == "Vermont") %>%
  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:plotUHI24hrTempPopAvgVT) and Table \@ref(tab:tabUHI24hrPopAvgVT) show population-weighted exposures for priority populations relative to the day-night average Urban Heat Island (UHI) effect for Vermont. All priority populations are exposed to the UHI effect on average. Indeed, the average UHI effect across Vermont is approximately `r round(mean(vt_blkgrpLST_sf$UHI24avg,na.rm=TRUE),1)`°. However, some populations experience significantly higher UHI effects. For example, limited English speaking households in Vermont 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 state day-night average UHI. Similarly, People of Color 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 state 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 state average.

```{r plotUHI24hrTempPopAvgVT, echo=FALSE, message=FALSE, warning=FALSE, fig.align = "center", fig.cap="Population-weighted average exposures to day-night average Urban Heat Island effect for priority populations in Vermont."}
UHI24hrPop %>% 
  mutate(Group = recode(Group, "Minority" = "People of Color",
                        "No HS Dip" = "No HS Diploma")) %>%
  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 = 10, label = "Above\nrural average"),
  #           color = "gray48") +
  # geom_segment(aes(x = "Over 64", xend = "Over 64", y = 8, yend = 11), 
  #              arrow = arrow(length = unit(0.3,"cm"))) +
  expand_limits(y = c(-1,5))

ggsave("images/VT_UHI24_graph.png")
```


```{r tabUHI24hrPopAvgVT, fig.align = "center", fig.cap="Population-weighted average exposures to day-night average Urban Heat Island effect for priority populations in Vermont relative to the rural average."}
UHI24hrPopTab %>% 
  mutate(Group = recode(Group, "Minority" = "People of Color",
                        "No HS Dip" = "No HS Diploma")) %>%
  arrange(Group) %>% 
  kableExtra::kable(longtable = T, booktabs = T, digits = 2, 
                    col.names = c("Group","Day-Night Avg UHI (°F)",
                                  "Difference from State UHI Avg* (°F)",
                                  "Pct Above/Below State 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("State average day-night UHI effect is ", round(mean(vt_blkgrpLST_sf$UHI24avg,na.rm=TRUE),2), "°"))
```

```{r LST24hrHighRiskStats, include=FALSE}
Pop24hr <- vt_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  filter(LST24EngLangScore == 5) %>% 
  summarize(sum(totalpopE,na.rm = TRUE))

Pop24hrPct <- round(Pop24hr/statepopsdf$TotalPop*100,1) %>% 
  pull()

HH24hr <- vt_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  filter(LST24EngLangScore == 5) %>% 
  summarize(sum(householdsE,na.rm = TRUE))

HH24hrPct <- round(HH24hr/statepopsdf$TotalHH*100,1) %>% 
  pull()

LEH24hr <- vt_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  filter(LST24EngLangScore == 5) %>% 
  summarize(sum(eng_limitE,na.rm = TRUE))

LEH24hrPct <- round(LEH24hr/statepopsdf$TotalLEH*100,1) %>% 
  pull()

MIN24hr <- vt_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  filter(LST24MinorScore == 5) %>% 
  summarize(sum(minorityE,na.rm = TRUE))

MIN24hrPct <- round(MIN24hr/statepopsdf$TotalMin*100,1) %>% 
  pull()
```


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 Vermont 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 People of Color in the top quintile for the state.[^LST24quintile] Approximately `r LEH24hr %>% format(., big.mark = ",") %>% pull()` limited English speaking households and `r MIN24hr %>% format(., big.mark = ",") %>% pull()` People of Color were resident in these highest quintile Block Groups. These numbers represent `r LEH24hrPct`% and `r MIN24hrPct`%, respectively, of those populations in the state. By comparison, `r HH24hrPct`% of all households, and `r Pop24hrPct`% of the general population were resident in these same Block Groups. 

Tables \@ref(tab:tabLST24hrHighRiskLEH) and \@ref(tab:tabLST24hrHighRiskMIN) show breakdowns by municipality of how many limited English speaking households and People of Color these represent, and the day-night average LST values for those same Census Block Groups. Maps and tables for other priority populations 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 <- vt_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.1 * xrange) # xmax - right
bbox_new[2] <- bbox_new[2] - (0.1 * 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_borders(lwd = 0.2, alpha = 0.8) + 
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ny_state_sf) + tm_borders(lwd = 0.2, alpha = 0.8) +
  tm_shape(vt_highways) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(vt_highways2nd) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.15) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.1,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/vt_LST24EngLang5.png",
          height = 8, width = 8, units = "in", dpi = 600)

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

\pagebreak
```{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."}
LST24EngLang5 %>% 
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Limited English Households` = sum(eng_limitE,na.rm = TRUE),
            `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)) %>% 
  rename("City/Town" = NAME) %>% 
  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),"°")) %>% 
  arrange(desc(`Limited English Households`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Limited English Households in Highest Quintile of Day-Night Avg LST", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header","striped")) %>% 
  kableExtra::column_spec(2:4, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 4, "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 People of Color by Census Block Group, July-Aug 2019."}
LST24Minority5 <- vt_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.1 * xrange) # xmax - right
bbox_new[2] <- bbox_new[2] - (0.1 * 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_borders(lwd = 0.2, alpha = 0.8) + 
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ny_state_sf) + tm_borders(lwd = 0.2, alpha = 0.8) +
  tm_shape(vt_highways) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(vt_highways2nd) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.15) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.1,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\nPeople of Color\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/vt_LST24Minority5.png",
          height = 8, width = 8, units = "in", dpi = 600)

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

\pagebreak
```{r tabLST24hrHighRiskMIN, fig.align = "center", fig.cap="People of Color in Census Block Groups at or above 80th percentile of day-night average Land Surface Temperature."}
# Identify municipalities and values by municipality
LST24Minority5 %>% 
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            Minorities = sum(minorityE,na.rm = TRUE),
            `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)) %>% 
  rename("City/Town" = NAME) %>% 
  mutate(`Pct of Minorities` = paste0(`Pct of Minorities`,"%"),
         Min = paste0(round(Min,2),"°"),
         Mean = paste0(round(Mean,2),"°"),
         Max = paste0(round(Max,2),"°")) %>% 
  arrange(desc(Minorities)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "People of Color in Highest Quintile of Day-Night Avg LST by Municipality", align = "r", col.names = c(names(.)[1:2], "People of Color", "Pct of POC", names(.)[5:7])) %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header","striped")) %>% 
  kableExtra::column_spec(2:4, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 4, "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. Daytime average LST for the state was `r round(mean(vt_blkgrpLST_sf$meanDayLST,na.rm=TRUE),1)`°. Daytime Land Surface Temperatures (LST) varied significantly across Vermont, with highest temperatures apparent in the most urbanized areas of the state, particularly in and around Burlington and south along the shore of Lake Champlain. Other areas of high temperature appear around Brattleboro, Bennington, Berlin, and Manchester (see Figure \@ref(fig:mapLSTdayVT) below). Across Vermont, approximately `r vt_blkgrpLST_sf %>% as.data.frame() %>% filter(pctileDayLST >= 0.8) %>% group_by() %>% summarize(sum(totalpopE, na.rm = TRUE)/10^3) %>% pull() %>% round(.,0)` thousand people live in Block Groups that exceeded the 80th percentile (temperatures over `r round(quantile(vt_blkgrpLST_sf$meanDayLST,0.8,na.rm=TRUE),1)`°) for daytime average temperatures.  

```{r mapLSTdayVT, fig.align = "center", fig.cap="Map of July-August 2019 average daytime Land Surface Temperatures (LST) across Vermont at Census Block Group level.", strip.white=TRUE}
# Map of average LST across Vermont
m <- tm_shape(vt_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(lwd = 0.2, alpha = 0.8) + 
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ny_state_sf) + tm_borders(lwd = 0.2, alpha = 0.8) +
  tm_shape(vt_highways) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(vt_highways2nd) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.15) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), 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/vt_DayLST.png",
          height = 8, width = 8, units = "in", dpi = 600)

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

```{r LSTDayoutliers, include=FALSE}
# Identify municipalities with high heat outliers
vt_LSTDayoutliers <- vt_blkgrpLST_sf %>% 
  arrange(desc(meanAvgLST)) %>% 
  slice(1:10) %>% 
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  group_by(NAME) %>%
  summarize(unique(NAME)) %>% 
  pull(NAME) %>% 
  paste(.,sep = " ", collapse = ", ")
```

Figure \@ref(fig:boxplotLSTdayVT) is a boxplot of average daytime temperatures by Block Group . The box represents temperatures ranging between the 25th and 75th percentiles. The line that divides the box into 2 parts represents the median temperature for all Block Groups, which in this case is `r round(median(vt_blkgrpLST_sf$meanDayLST, na.rm = TRUE),2)`°. Half of the state's Block Groups are below the median and half are above the median. Each dot represents an individual Block Group. The top ten highest outlier Block Groups (far right) were found in the following municipalities: `r vt_LSTDayoutliers`.

```{r boxplotLSTdayVT, fig.align = "center", fig.cap="Boxplot of July-Aug 2019 daytime average Land Surface Temperatures (LST) across Vermont at Census Block Group level."}
# boxplot of PM2.5 by state for 2016
vt_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  drop_na(meanDayLST) %>% 
  ggplot(aes(x = STATE, y = meanDayLST, fill = STATE)) + 
  geom_boxplot() + 
  geom_jitter(color = "black", size = 0.4, alpha = 0.6) +
  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 Vermont are most pronounced in the more densely urbanized areas (see Figure \@ref(fig:mapUHIdayVT)), particularly in the Burlington area, and in the cores of Berlin, Bennington, Rutland, and St. Albans. Temperatures in these areas are up to `r round(max(vt_blkgrpLST_sf$UHIDay,na.rm = TRUE),1)`° warmer than the rural background average.  

```{r mapUHIdayVT, fig.align = "center", fig.cap="Map of July-Aug 2019 daytime Urban Heat Island (UHI) effect across Vermont 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 Vermont
m <- tm_shape(vt_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(lwd = 0.2, alpha = 0.8) + 
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ny_state_sf) + tm_borders(lwd = 0.2, alpha = 0.8) +
  tm_shape(vt_highways) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(vt_highways2nd) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.15) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), 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/vt_DayUHI.png",
          height = 8, width = 8, units = "in", dpi = 600)

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

```{r LSTDayUHIoutliers, include=FALSE}
vt_LSTDayUHIoutliers <- vt_blkgrpLST_sf %>% 
  arrange(desc(UHIDay)) %>% 
  slice(1:10) %>% 
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  group_by(NAME) %>%
  summarize(unique(NAME)) %>% 
  pull(NAME) %>% 
  paste(.,sep = " ", collapse = ", ")
```

Figure \@ref(fig:boxplotUHIDayVT) is a boxplot of average daytime Urban Heat Island (UHI) effect values by Block Group . The box represents UHI values ranging between the 25th and 75th percentiles. The line that divides the box into 2 parts represents the median value for all Block Groups, which in this case is `r round(median(vt_blkgrpLST_sf$UHIDay, na.rm = TRUE),2)`°. Half of the state's Block Groups are below the median and half are above the median. Each dot represents an individual Block Group. The top ten highest outlier Block Groups (far right) were found in the following municipalities: `r vt_LSTDayUHIoutliers`. 

```{r boxplotUHIDayVT, fig.align = "center", fig.cap="Boxplot of July-Aug 2019 daytime Urban Heat Island (UHI) across Vermont 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
vt_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  drop_na(UHIDay) %>% 
  ggplot(aes(x = STATE, y = UHIDay, fill = STATE)) + 
  geom_boxplot() + 
  geom_jitter(color = "black", size = 0.4, alpha = 0.6) +
  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 Priority Populations

```{r LSTDayTempPop, include=FALSE}
# Pop Weighted avg of 24hr LST exposure for all Groups in Vermont relative to VT average
LSTDayPop <- vt_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  # filter(STATE == "Vermont") %>%
  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 Vermont relative to VT average
LSTDayPopTab <- vt_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  # filter(STATE == "Vermont") %>%
  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:plotLSTDayTempPopAvgVT) and Table \@ref(tab:tabLSTDayPopAvgVT) show population-weighted exposures for priority populations relative to the average daytime average LST for the state. For example, limited English speaking households in Vermont 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 state daytime average LST. Similarly, People of Color 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 state 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 state average.

```{r plotLSTDayTempPopAvgVT, fig.align = "center", fig.cap="Population-weighted average exposures to average daytime Land Surface Temperatures for priority populations in Vermont relative to the state average."}
LSTDayPop %>% 
  mutate(Group = recode(Group, "Minority" = "People of Color",
                        "No HS Dip" = "No HS Diploma")) %>%
  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 = 82, label = "Above state\naverage"),
            color = "gray48") +
  geom_text(aes(x = "Low Income", y = 80.3, label = "Below state\naverage"),
            color = "gray48") +
  # geom_segment(aes(x = "Under 5", xend = "Under 5", y = 94, yend = 97), 
  #              arrow = arrow(length = unit(0.3,"cm"))) +
  expand_limits(y = c(80,83.5))

ggsave("images/VT_LSTDay_graph.png")
```


```{r tabLSTDayPopAvgVT, fig.align = "center", fig.cap="Population-weighted average exposures to Daytime LST for priority populations in Vermont relative to the state average."}
LSTDayPopTab %>% 
  mutate(Group = recode(Group, "Minority" = "People of Color",
                        "No HS Dip" = "No HS Diploma")) %>%
  arrange(Group) %>% 
  kableExtra::kable(longtable = T, booktabs = T, digits = 2, 
                    col.names = c("Group","Avg Day LST (°F)",
                                  "Difference from State Avg* (°F)",
                                  "Pct Above/Below State 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("State average daytime LST is ", round(mean(vt_blkgrpLST_sf$meanDayLST,na.rm=TRUE),2), "°"))
```

```{r UHIDayPop, include=FALSE}
# Pop Weighted avg of 24hr LST exposure for all Groups in Vermont relative to VT average
UHIDayPop <- vt_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  # filter(STATE == "Vermont") %>%
  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 Vermont relative to VT average
UHIDayPopTab <- vt_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  # filter(STATE == "Vermont") %>%
  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:plotUHIDayTempPopAvgVT) and Table \@ref(tab:tabUHIDayPopAvgVT) show population-weighted exposures for priority populations relative to the average daytime average Urban Heat Island (UHI) effect for the state. All priority populations are exposed to the UHI effect on average. Indeed, the average UHI effect across the state is approximately `r round(mean(vt_blkgrpLST_sf$UHIDay,na.rm=TRUE),1)`°. However, some populations experience significantly higher UHI effects. For example, limited English speaking households in Vermont 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 state daytime average UHI. Similarly, People of Color 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 state 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 state average.

```{r plotUHIDayTempPopAvgVT, echo=FALSE, message=FALSE, warning=FALSE, fig.align = "center", fig.cap="Population-weighted average exposures to daytime Urban Heat Island effect for priority populations in Vermont."}
UHIDayPop %>% 
  mutate(Group = recode(Group, "Minority" = "People of Color",
                        "No HS Dip" = "No HS Diploma")) %>%
  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 = 15, 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,6))

ggsave("images/VT_UHIDay_graph.png")
```


```{r tabUHIDayPopAvgVT, fig.align = "center", fig.cap="Population-weighted average exposures to daytime Urban Heat Island effect for priority populations in Vermont relative to the rural average."}
UHIDayPopTab %>% 
  mutate(Group = recode(Group, "Minority" = "People of Color",
                        "No HS Dip" = "No HS Diploma")) %>%
  arrange(Group) %>% 
  kableExtra::kable(longtable = T, booktabs = T, digits = 2, 
                    col.names = c("Group","Avg Daytime UHI (°F)",
                                  "Difference from State UHI Avg* (°F)",
                                  "Pct Above/Below State 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("State average daytime UHI effect is ", round(mean(vt_blkgrpLST_sf$UHIDay,na.rm=TRUE),2), "°"))
```

```{r LSTDayHighRiskStats, include=FALSE}
PopDay <- vt_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  filter(LSTDayEngLangScore == 5) %>% 
  summarize(sum(totalpopE,na.rm = TRUE))

PopDayPct <- round(PopDay/statepopsdf$TotalPop*100,1) %>% 
  pull()

HHDay <- vt_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  filter(LSTDayEngLangScore == 5) %>% 
  summarize(sum(householdsE,na.rm = TRUE))

HHDayPct <- round(HHDay/statepopsdf$TotalHH*100,1) %>% 
  pull()

LEHDay <- vt_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  filter(LSTDayEngLangScore == 5) %>% 
  summarize(sum(eng_limitE,na.rm = TRUE))

LEHDayPct <- round(LEHDay/statepopsdf$TotalLEH*100,1) %>% 
  pull()

MINDay <- vt_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  filter(LSTDayMinorScore == 5) %>% 
  summarize(sum(minorityE,na.rm = TRUE))

MINDayPct <- round(MINDay/statepopsdf$TotalMin*100,1) %>% 
  pull()
```

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 Vermont which exhibited daytime average LSTs in the top quintile (i.e., 80th percentile) and also contained high percentages of limited English speaking households or People of Color in the top quintile for the state.[^Dayquintile] Approximately `r LEHDay %>% format(., big.mark = ",") %>% pull()` limited English households and `r MINDay %>% format(., big.mark = ",") %>% pull()` People of Color were resident in these highest quintile Block Groups. These numbers represent `r LEHDayPct`% and `r MINDayPct`%, respectively, of those populations in the state. By comparison, `r HHDayPct`% of all households, and `r PopDayPct`% of the general population were resident in these same Block Groups.  

Tables \@ref(tab:tabLSTDayHighRiskLEH) and \@ref(tab:tabLSTDayHighRiskMIN) show breakdowns by municipality of the numbers of limited English households and People of Color resident in these top quintile Block Groups. Maps and tables for other priority populations 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 <- vt_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.2 * xrange) # xmin - left
bbox_new[3] <- bbox_new[3] + (0.3 * xrange) # xmax - right
bbox_new[2] <- bbox_new[2] - (0.1 * 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_borders(lwd = 0.2, alpha = 0.8) + 
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ny_state_sf) + tm_borders(lwd = 0.2, alpha = 0.8) +
  tm_shape(vt_highways) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(vt_highways2nd) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.15) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.1,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/vt_LSTDayEngLang5.png",
          height = 8, width = 8, units = "in", dpi = 600)

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

\pagebreak
```{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."}
LSTDayEngLang5 %>% 
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Limited English Households` = sum(eng_limitE,na.rm = TRUE),
            `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)) %>% 
  rename("City/Town" = NAME) %>% 
  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),"°")) %>% 
  arrange(desc(`Limited English Households`)) %>% 
  kableExtra::kable(longtable = T, 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","striped")) %>% 
  kableExtra::column_spec(2:4, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 4, "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 People of Color by Census Block Group, July-Aug 2019."}
LSTDayMinority5 <- vt_blkgrpLST_sf %>%
  filter(LSTDayMinorScore == 5)

bbox_new <- st_bbox(LSTDayMinority5) # 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.2 * xrange) # xmin - left
bbox_new[3] <- bbox_new[3] + (0.5 * xrange) # xmax - right
bbox_new[2] <- bbox_new[2] - (0.1 * 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 <- LSTDayMinority5 %>%
  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_borders(lwd = 0.2, alpha = 0.8) + 
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ny_state_sf) + tm_borders(lwd = 0.2, alpha = 0.8) +
  tm_shape(vt_highways) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(vt_highways2nd) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.15) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.1,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\nPeople of Color\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/vt_LSTDayMinority5.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/vt_LSTDayMinority5.png")
```
\pagebreak
```{r tabLSTDayHighRiskMIN, fig.align = "center", fig.cap="People of Color in Census Block Groups at or above 80th percentile of daytime average Land Surface Temperature."}
LSTDayMinority5 %>% 
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            Minorities = sum(minorityE,na.rm = TRUE),
            `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)) %>% 
  rename("City/Town" = NAME) %>% 
  mutate(`Pct of Minorities` = paste0(`Pct of Minorities`,"%"),
         Min = paste0(round(Min,2),"°"),
         Mean = paste0(round(Mean,2),"°"),
         Max = paste0(round(Max,2),"°")) %>% 
  arrange(desc(Minorities)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "People of Color in Highest Quintile of Daytime Avg LST by Municipality", align = "r", col.names = c(names(.)[1:2], "People of Color", "Pct of POC", names(.)[5:7])) %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header","striped")) %>% 
  kableExtra::column_spec(2:4, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 4, "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 average LST for the state was `r round(mean(vt_blkgrpLST_sf$meanNightLST,na.rm=TRUE),1)`°. Nighttime Land Surface Temperatures (LST) varied significantly across Vermont, with highest temperatures apparent in along the shore of Lake Champlain and in the south around Bennington and Brattleboro (see Figure \@ref(fig:mapLSTnightVT) below).

```{r mapLSTnightVT, fig.align = "center", fig.cap="Map of July-August 2019 average nighttime Land Surface Temperatures (LST) across Vermont at Census Block Group level.", strip.white=TRUE}
# Map of average LST across Vermont
m <- tm_shape(vt_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(lwd = 0.2, alpha = 0.8) + 
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ny_state_sf) + tm_borders(lwd = 0.2, alpha = 0.8) +
  tm_shape(vt_highways) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(vt_highways2nd) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.15) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), 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/vt_NightLST.png",
          height = 8, width = 8, units = "in", dpi = 600)

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

```{r LSTNightoutliers, include=FALSE}
vt_LSTNightoutliers <- vt_blkgrpLST_sf %>% 
  arrange(desc(meanNightLST)) %>% 
  slice(1:10) %>% 
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  group_by(NAME) %>%
  summarize(unique(NAME)) %>% 
  pull(NAME) %>% 
  paste(.,sep = " ", collapse = ", ")
```

Figure \@ref(fig:boxplotLSTnightVT) is a boxplot of average nighttime temperatures by Block Group . The box represents temperatures ranging between the 25th and 75th percentiles. The line that divides the box into 2 parts represents the median temperature for all Block Groups, which in this case is `r round(median(vt_blkgrpLST_sf$meanNightLST, na.rm = TRUE),2)`°. Half of the state's Block Groups are below the median and half are above the median. Each dot represents an individual Block Group. Large dots on either end represents outliers, or unusually high or low values. The top ten highest outlier Block Groups (far right) were found in the following municipalities: `r vt_LSTNightoutliers`. 

```{r boxplotLSTnightVT, fig.align = "center", fig.cap="Boxplot of July-Aug 2019 nighttime Land Surface Temperatures (LST) across Vermont at Census Block Group level."}
# boxplot of PM2.5 by state for 2016
vt_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  drop_na(meanNightLST) %>% 
  ggplot(aes(x = STATE, y = meanNightLST, fill = STATE)) + 
  geom_boxplot() + 
  geom_jitter(color = "black", size = 0.4, alpha = 0.6) +
  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 Vermont are most pronounced in and around Burlington and St. Albans (see Figure \@ref(fig:mapUHInightVT)). Temperatures in these areas are up to `r round(max(vt_blkgrpLST_sf$UHINight,na.rm = TRUE),1)`° warmer than the rural background average. By contrast, most of the rest of the state is close to or cooler than rural background average. 

```{r mapUHInightVT, fig.align = "center", fig.cap="Map of July-Aug 2019 nighttime Urban Heat Island (UHI) effect across Vermont 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 Vermont
m <- tm_shape(vt_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(lwd = 0.2, alpha = 0.8) + 
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ny_state_sf) + tm_borders(lwd = 0.2, alpha = 0.8) +
  tm_shape(vt_highways) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(vt_highways2nd) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.15) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), 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/vt_NightUHI.png",
          height = 8, width = 8, units = "in", dpi = 600)

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

```{r LSTNightUHIoutliers, include=FALSE}
vt_LSTNightoutliers <- vt_blkgrpLST_sf %>% 
  arrange(desc(UHINight)) %>% 
  slice(1:10) %>% 
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  group_by(NAME) %>%
  summarize(unique(NAME)) %>% 
  pull(NAME) %>% 
  paste(.,sep = " ", collapse = ", ")
```


Figure \@ref(fig:boxplotUHINightVT) is a boxplot of nighttime Urban Heat Island (UHI) effects by Block Group . The box represents values ranging between the 25th and 75th percentiles. The line that divides the box into 2 parts represents the median temperature for all Block Groups, which in this case is `r round(median(vt_blkgrpLST_sf$UHINight, na.rm = TRUE),2)`. Half of the state's Block Groups are below the median and half are above the median. Each dot represents an individual Block Group. Large black dots represent outliers, or unusually high or low values. The top ten highest outlier Block Groups (far right) were found in the following municipalities: `r vt_LSTNightoutliers`.

```{r boxplotUHINightVT, fig.align = "center", fig.cap="Boxplot of July-Aug 2019 nighttime Urban Heat Island (UHI) across Vermont 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
vt_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  drop_na(UHINight) %>% 
  ggplot(aes(x = STATE, y = UHINight, fill = STATE)) + 
  geom_boxplot() + 
  geom_jitter(color = "black", size = 0.4, alpha = 0.6) +
  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 Priority Populations

```{r LSTNightPop, include=FALSE}
# Pop Weighted avg of 24hr LST exposure for all Groups in Vermont relative to VT average
LSTNightPop <- vt_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  # filter(STATE == "Vermont") %>%
  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 Vermont relative to VT average
LSTNightPopTab <- vt_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  # filter(STATE == "Vermont") %>%
  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:plotLSTNightTempPopAvgVT) and Table \@ref(tab:tabLSTNightPopAvgVT) show population-weighted exposures for priority populations relative to the average nighttime average LST of `r round(mean(vt_blkgrpLST_sf$meanNightLST,na.rm=TRUE),1)`° for the state. For example, limited English speaking households in Vermont 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 state nighttime average LST. Similarly, People of Color 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 state 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 state average.

```{r plotLSTNightTempPopAvgVT, fig.align = "center", fig.cap="Population-weighted average exposures to average nighttime Land Surfact Temperatures for priority populations in Vermont relative to the state average."}
LSTNightPop %>% 
  mutate(Group = recode(Group, "Minority" = "People of Color",
                        "No HS Dip" = "No HS Diploma")) %>%
  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 = 60.5, label = "Above state\naverage"),
            color = "gray48") +
  geom_text(aes(x = "Low Income", y = 58.5, label = "Below state\naverage"),
            color = "gray48") +
  # geom_segment(aes(x = "Under 5", xend = "Under 5", y = 66.9, yend = 67.5), 
  #              arrow = arrow(length = unit(0.3,"cm"))) +
  expand_limits(y = c(58,62))

ggsave("images/VT_LSTNight_graph.png")
```


```{r tabLSTNightPopAvgVT, fig.align = "center", fig.cap="Population-weighted average exposures to nighttime LST for priority populations in Vermont relative to the state average."}
LSTNightPopTab %>% 
  mutate(Group = recode(Group, "Minority" = "People of Color",
                        "No HS Dip" = "No HS Diploma")) %>%
  arrange(Group) %>% 
  kableExtra::kable(longtable = T, booktabs = T, digits = 2, 
                    col.names = c("Group","Avg Night LST (°F)",
                                  "Difference from State Avg* (°F)",
                                  "Pct Above/Below State 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("State average nighttime LST is ", round(mean(vt_blkgrpLST_sf$meanNightLST,na.rm=TRUE),2), "°"))
```

```{r UHINightPop, include=FALSE}
# Pop Weighted avg of 24hr LST exposure for all Groups in Vermont relative to VT average
UHINightPop <- vt_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  # filter(STATE == "Vermont") %>%
  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 Vermont relative to VT average
UHINightPopTab <- vt_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  # filter(STATE == "Vermont") %>%
  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:plotUHINightTempPopAvgVT) and Table \@ref(tab:tabUHINightPopAvgVT) show population-weighted exposures for priority populations relative to the average daytime average Urban Heat Island (UHI) effect for the state. All priority populations are exposed to the UHI effect on average. Indeed, the average UHI effect across the state is approximately `r round(mean(vt_blkgrpLST_sf$UHINight,na.rm=TRUE),1)`°. However, some populations experience significantly higher UHI effects. For example, limited English speaking households in Vermont 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 state daytime average UHI. Similarly, People of Color 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 state 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 state average.

```{r plotUHINightTempPopAvgVT, echo=FALSE, message=FALSE, warning=FALSE, fig.align = "center", fig.cap="Population-weighted average exposures to nighttime Urban Heat Island effect for priority populations in Vermont."}
UHINightPop %>% 
  mutate(Group = recode(Group, "Minority" = "People of Color",
                        "No HS Dip" = "No HS Diploma")) %>%
  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 = 4.3, label = "Above\nrural average"),
  #           color = "gray48") +
  # geom_segment(aes(x = "Over 64", xend = "Over 64", y = 3.7, yend = 4.4), 
  #              arrow = arrow(length = unit(0.3,"cm"))) +
  expand_limits(y = c(-1,4))

ggsave("images/VT_UHINight_graph.png")
```


```{r tabUHINightPopAvgVT, fig.align = "center", fig.cap="Population-weighted average exposures to nighttime Urban Heat Island effect for priority populations in Vermont relative to the rural average."}
UHINightPopTab %>% 
  mutate(Group = recode(Group, "Minority" = "People of Color",
                        "No HS Dip" = "No HS Diploma")) %>%
  arrange(Group) %>% 
  kableExtra::kable(longtable = T, booktabs = T, digits = 2, 
                    col.names = c("Group","Avg Nighttime UHI (°F)",
                                  "Difference from State UHI Avg (°F)",
                                  "Pct Above/Below State 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("State average nighttime UHI effect is ", round(mean(vt_blkgrpLST_sf$UHINight,na.rm=TRUE),2), "°"))
```

```{r LSTNightHighRiskStats, include=FALSE}
PopNight <- vt_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  filter(LSTNightEngLangScore == 5) %>% 
  summarize(sum(totalpopE,na.rm = TRUE))

PopNightPct <- round(PopNight/statepopsdf$TotalPop*100,1) %>% 
  pull()

HHNight <- vt_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  filter(LSTNightEngLangScore == 5) %>% 
  summarize(sum(householdsE,na.rm = TRUE))

HHNightPct <- round(HHNight/statepopsdf$TotalHH*100,1) %>% 
  pull()

LEHNight <- vt_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  filter(LSTNightEngLangScore == 5) %>% 
  summarize(sum(eng_limitE,na.rm = TRUE))

LEHNightPct <- round(LEHNight/statepopsdf$TotalLEH*100,1) %>% 
  pull()

MINNight <- vt_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  filter(LSTNightMinorScore == 5) %>% 
  summarize(sum(minorityE,na.rm = TRUE))

MINNightPct <- round(MINNight/statepopsdf$TotalMin*100,1) %>% 
  pull()
```

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:mapLSTNightHighRiskLEH) and \@ref(fig:mapLSTNightHighRiskMIN) highlight Census Block Groups across Vermont 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 People of Color in the top quintile for the state.[^Nightquintile] Approximately `r LEHNight %>% format(., big.mark = ",") %>% pull()` limited English households and `r MINNight %>% format(., big.mark = ",") %>% pull()` People of Color were resident in these highest quintile Block Groups. These numbers represent `r LEHNightPct`% and `r MINNightPct`%, respectively, of those populations in the state. By comparison, `r HHNightPct`% of all households, and `r PopNightPct`% of the general population were resident in these same Block Groups. 

Tables \@ref(tab:tabLSTNightHighRiskLEH) and \@ref(tab:tabLSTNightHighRiskMIN) show breakdowns by municipality of how many limited English households and People of Color these represent, and the day-night average LST values for those same Census Block Groups. Maps and tables for other priority populations 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 mapLSTNightHighRiskLEH, 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 Households by Census Block Group, July-Aug 2019."}
LSTNightEngLang5 <- vt_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.1 * 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_borders(lwd = 0.2, alpha = 0.8) + 
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ny_state_sf) + tm_borders(lwd = 0.2, alpha = 0.8) +
  tm_shape(vt_highways) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(vt_highways2nd) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.15) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.1,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/vt_LSTNightEngLang5.png",
          height = 8, width = 8, units = "in", dpi = 600)

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


```{r tabLSTNightHighRiskLEH, fig.align = "center", fig.cap="Limited English Speaking Housheolds in Census Block Groups at or above 80th percentile of nighttime average Land Surface Temperature."}
LSTNightEngLang5 %>% 
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Limited English Households` = sum(eng_limitE,na.rm = TRUE),
            `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)) %>% 
  rename("City/Town" = NAME) %>% 
  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),"°")) %>% 
  arrange(desc(`Limited English Households`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Limited English Speaking Households in Highest Quintile of Nighttime Avg LST", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header","striped")) %>% 
  kableExtra::column_spec(2:4, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 4, "Night 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 People of Color by Census Block Group, July-Aug 2019."}
LSTNightMinority5 <- vt_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.25 * xrange) # xmax - right
bbox_new[2] <- bbox_new[2] - (0.1 * 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 <- 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_borders(lwd = 0.2, alpha = 0.8) + 
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ny_state_sf) + tm_borders(lwd = 0.2, alpha = 0.8) +
  tm_shape(vt_highways) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(vt_highways2nd) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.15) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.1,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\nPeople of Color\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/vt_LSTNightMinority5.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/vt_LSTNightMinority5.png")
```
\pagebreak
```{r tabLSTNightHighRiskMIN, fig.align = "center", fig.cap="People of Color in Census Block Groups at or above 80th percentile of Nighttime average Land Surface Temperature."}
LSTNightMinority5 %>% 
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            Minorities = sum(minorityE,na.rm = TRUE),
            `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)) %>% 
  rename("City/Town" = NAME) %>% 
  mutate(`Pct of Minorities` = paste0(`Pct of Minorities`,"%"),
         Min = paste0(round(Min,2),"°"),
         Mean = paste0(round(Mean,2),"°"),
         Max = paste0(round(Max,2),"°")) %>% 
  arrange(desc(Minorities)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "People of Color in Highest Quintile of Nighttime Avg LST by Municipality", align = "r", col.names = c(names(.)[1:2], "People of Color", "Pct of POC", names(.)[5:7])) %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header","striped")) %>% 
  kableExtra::column_spec(2:4, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 4, "Night 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, VT." *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(vt_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(vt_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(vt_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(vt_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(vt_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(vt_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 Vermont. 

### 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_vt_DayNight_sf <- readRDS("DATA/LST/isd_ne_DayNight_sf.Rds") %>% 
  st_transform(., crs = st_crs(vt_state_sf_cb)) %>% 
  st_intersection(vt_state_sf_cb,.) 
```

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 Vermont. 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_vt_DayNight_sf)` weather stations across Vermont 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 Vermont."}
tm_shape(vt_state_sf_cb, unit = "mi") + tm_borders(lwd = 0.2) +
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(isd_vt_DayNight_sf) + tm_dots(shape = 4, size = 0.08) +
   tm_scale_bar(breaks = c(0, 25, 50), text.size = 0.5, 
               position = c(0.35,0.005)) +
  tm_layout(title = "NOAA ISD Weather Stations",
            frame = FALSE, main.title.size = 0.8,
            title.position = c("center", "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, although nighttime values are not statistically significant (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_vt_DayNight_sf$gridDay, y = isd_vt_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_vt_DayNight_sf$gridNight, y = isd_vt_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_vt_DayNight_sf$gridDayNight, y = isd_vt_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_vt_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_vt_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_vt_DayNight_sf)
RMSE_daynight <- sqrt(mean(lm_daynight$residuals^2))
```

Linear models are statistically significant predictors of air temperature data for day and day-night values (see Figure \@ref(fig:LSTvalidation)). The Root Mean Square Error (RMSE) of day values, `r round(RMSE_day,2)`°, and Day-Night Avg LST, `r round(RMSE_daynight,2)`°, indicate modest deviation of predicted air temperatures from observed values. 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 Vermont."}

day <- isd_vt_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("ISD 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=93,y=77) +
  expand_limits(y=c(70,95))

night <- isd_vt_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("ISD 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=57,y=72) +
  expand_limits(y=c(55,80))

daynight <- isd_vt_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("ISD 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=73,y=67) +
  expand_limits(y=c(60,85))

library(gridExtra)
# grid.arrange(day,night,daynight,ncol=1,top="LST vs Weather Station Air Temp", heights = unit(rep(3,3),rep("inches",3)),
#              widths = unit(5,"inches"))

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/vt_validate.png", 
       width = 5, height = 8, units = "in")
knitr::include_graphics("DATA/LST/vt_validate.png")
```


### Urban Heat Island effects {-}

Urban Heat Island (UHI) effects are measured as the difference between LST raster values across Vermont for each time period and the average LST for rural states 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 states of Vermont 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 state. 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 states of Vermont 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 Vermont not within the urban areas polygons were then extracted to represent rural states. 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 size is 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 Census. 

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

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

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

## 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 People of Color 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", col.names = c(names(.)[1:2], "People of Color", "Non-POC", names(.)[5])) %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header"))
```

Population numbers of People of Color, 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 People of Color occupy Block Groups with higher temperatures than the simple average would indicate?

To calculate the population-weighted average temperature exposure for People of Color, the number of People of Color in each Block Group is used as a 'weight'. The temperature of each Block Group is multiplied by its respective number of People of Color. 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,
         TempxPOC = 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",
                    col.names = c(names(.)[1:2], "People of Color", "Non-POC", names(.)[5:6])) %>% 
  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 People of Color. 

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",
                    col.names = c(names(.)[1:2], "People of Color", "Non-POC", names(.)[5:6])) %>% 
  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 People of Color experience a higher population-weighted average temperature than the general or total population. 


# Appendix B: Supplementary Figures {-}

## Day-Night Highest Exposure Priority Populations {-}
```{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 <- vt_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.15 * xrange) # xmin - left
bbox_new[3] <- bbox_new[3] + (0.2 * xrange) # xmax - right
bbox_new[2] <- bbox_new[2] - (0.1 * 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_borders(lwd = 0.2, alpha = 0.8) + 
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ny_state_sf) + tm_borders(lwd = 0.2, alpha = 0.8) +
  tm_shape(vt_highways) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(vt_highways2nd) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.15) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.1,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/vt_LST24LowInc5.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/vt_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."}
LST24LowInc5 %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Low Income Persons` = sum(num2povE,na.rm = TRUE),
            `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)) %>% 
  rename("City/Town" = NAME) %>% 
  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),"°")) %>% 
  arrange(desc(`Low Income Persons`)) %>% 
  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","striped")) %>% 
  kableExtra::column_spec(2:4, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 4, "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 <- vt_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.1 * 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_borders(lwd = 0.2, alpha = 0.8) + 
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ny_state_sf) + tm_borders(lwd = 0.2, alpha = 0.8) +
  tm_shape(vt_highways) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(vt_highways2nd) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.15) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.1,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/vt_LST24NoHSDip5.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/vt_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."}
LST24NoHSDip5 %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `No HS Dip` = sum(lthsE,na.rm = TRUE),
            `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)) %>% 
  rename("City/Town" = NAME) %>% 
  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),"°")) %>% 
  arrange(desc(`No HS Dip`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Adults without a High School Diploma in Highest Quintile of Day-Night Avg LST", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header","striped")) %>% 
  kableExtra::column_spec(2:4, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 4, "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 <- vt_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.1 * xrange) # xmin - left
# bbox_new[3] <- bbox_new[3] + (0.25 * xrange) # xmax - right
bbox_new[2] <- bbox_new[2] - (0.1 * 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_borders(lwd = 0.2, alpha = 0.8) + 
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ny_state_sf) + tm_borders(lwd = 0.2, alpha = 0.8) +
  tm_shape(vt_highways) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(vt_highways2nd) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.15) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.1,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/vt_LST24Over645.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/vt_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."}
LST24Over645 %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Over 64` = sum(over64E,na.rm = TRUE),
            `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)) %>% 
  rename("City/Town" = NAME) %>% 
  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),"°")) %>% 
  arrange(desc(`Over 64`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Persons Over 64 in Highest Quintile of Day-Night Avg LST", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header","striped")) %>% 
  kableExtra::column_spec(2:4, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 4, "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 <- vt_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.2 * xrange) # xmin - left
bbox_new[3] <- bbox_new[3] + (0.1 * xrange) # xmax - right
bbox_new[2] <- bbox_new[2] - (0.1 * 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_borders(lwd = 0.2, alpha = 0.8) + 
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ny_state_sf) + tm_borders(lwd = 0.2, alpha = 0.8) +
  tm_shape(vt_highways) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(vt_highways2nd) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.15) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.1,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/vt_LST24Under55.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/vt_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."}
LST24Under55 %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Under 5` = sum(under5E,na.rm = TRUE),
            `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)) %>% 
  rename("City/Town" = NAME) %>% 
  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),"°")) %>% 
  arrange(desc(`Under 5`)) %>% 
  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","striped")) %>% 
  kableExtra::column_spec(2:4, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 4, "Day-Night LST (°F)" = 3))
```


## Daytime Highest Exposure Priority Populations {-}

```{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 <- vt_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.2 * xrange) # xmin - left
bbox_new[3] <- bbox_new[3] + (0.1 * xrange) # xmax - right
bbox_new[2] <- bbox_new[2] - (0.1 * 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_borders(lwd = 0.2, alpha = 0.8) + 
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ny_state_sf) + tm_borders(lwd = 0.2, alpha = 0.8) +
  tm_shape(vt_highways) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(vt_highways2nd) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.15) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.1,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/vt_LSTDayLowInc5.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/vt_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."}
LSTDayLowInc5 %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Low Income Persons` = sum(num2povE,na.rm = TRUE),
            `Pct of Low Income Persons` = round(sum(num2povE)/max(TotalLowInc)*100,1),
            Min = min(meanDayLST,na.rm=TRUE),
            Mean = mean(meanDayLST,na.rm=TRUE),
            Max = max(meanDayLST,na.rm=TRUE)) %>% 
  rename("City/Town" = NAME) %>% 
  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),"°")) %>% 
  arrange(desc(`Low Income Persons`)) %>% 
  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","striped")) %>% 
  kableExtra::column_spec(2:4, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 4, "Day 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 <- vt_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.2 * xrange) # xmin - left
bbox_new[3] <- bbox_new[3] + (0.1 * xrange) # xmax - right
bbox_new[2] <- bbox_new[2] - (0.1 * 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_borders(lwd = 0.2, alpha = 0.8) + 
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ny_state_sf) + tm_borders(lwd = 0.2, alpha = 0.8) +
  tm_shape(vt_highways) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(vt_highways2nd) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.15) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.1,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/vt_LSTDayNoHSDip5.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/vt_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."}
LSTDayNoHSDip5 %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `No HS Dip` = sum(lthsE,na.rm = TRUE),
            `Pct of No HS Dip` = round(sum(lthsE)/max(TotalNoHS)*100,1),
            Min = min(meanDayLST,na.rm=TRUE),
            Mean = mean(meanDayLST,na.rm=TRUE),
            Max = max(meanDayLST,na.rm=TRUE)) %>% 
  rename("City/Town" = NAME) %>% 
  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),"°")) %>% 
  arrange(desc(`No HS Dip`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Adults without a High School Diploma in Highest Quintile of Daytime Avg LST", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header","striped")) %>% 
  kableExtra::column_spec(2:4, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 4, "Day 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 <- vt_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.2 * xrange) # xmin - left
bbox_new[3] <- bbox_new[3] + (0.1 * xrange) # xmax - right
bbox_new[2] <- bbox_new[2] - (0.1 * 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_borders(lwd = 0.2, alpha = 0.8) + 
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ny_state_sf) + tm_borders(lwd = 0.2, alpha = 0.8) +
  tm_shape(vt_highways) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(vt_highways2nd) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.15) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.1,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/vt_LSTDayOver645.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/vt_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."}
LSTDayOver645 %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Over 64` = sum(over64E,na.rm = TRUE),
            `Pct of Over 64` = round(sum(over64E)/max(TotalOver64)*100,1),
            Min = min(meanDayLST,na.rm=TRUE),
            Mean = mean(meanDayLST,na.rm=TRUE),
            Max = max(meanDayLST,na.rm=TRUE)) %>% 
  rename("City/Town" = NAME) %>% 
  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),"°")) %>% 
  arrange(desc(`Over 64`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Persons Over 64 in Highest Quintile of Daytime LST", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header","striped")) %>% 
  kableExtra::column_spec(2:4, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 4, "Day 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 <- vt_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.2 * xrange) # xmin - left
bbox_new[3] <- bbox_new[3] + (0.1 * xrange) # xmax - right
bbox_new[2] <- bbox_new[2] - (0.1 * 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_borders(lwd = 0.2, alpha = 0.8) + 
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ny_state_sf) + tm_borders(lwd = 0.2, alpha = 0.8) +
  tm_shape(vt_highways) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(vt_highways2nd) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.15) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.1,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/vt_LSTDayUnder55.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/vt_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."}
LSTDayUnder55 %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Under 5` = sum(under5E,na.rm = TRUE),
            `Pct of Under 5` = round(sum(under5E)/max(TotalUnder5)*100,1),
            Min = min(meanDayLST,na.rm=TRUE),
            Mean = mean(meanDayLST,na.rm=TRUE),
            Max = max(meanDayLST,na.rm=TRUE)) %>% 
  rename("City/Town" = NAME) %>% 
  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),"°")) %>% 
  arrange(desc(`Under 5`)) %>% 
  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","striped")) %>% 
  kableExtra::column_spec(2:4, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 4, "Day LST (°F)" = 3))
```


## Nighttime Highest Exposure Priority Populations {-}

```{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 <- vt_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.2 * xrange) # xmin - left
bbox_new[3] <- bbox_new[3] + (0.15 * xrange) # xmax - right
bbox_new[2] <- bbox_new[2] - (0.1 * 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_borders(lwd = 0.2, alpha = 0.8) + 
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ny_state_sf) + tm_borders(lwd = 0.2, alpha = 0.8) +
  tm_shape(vt_highways) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(vt_highways2nd) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.15) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.1,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/vt_LSTNightLowInc5.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/vt_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."}
LSTNightLowInc5 %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Low Income Persons` = sum(num2povE,na.rm = TRUE),
            `Pct of Low Income Persons` = round(sum(num2povE)/max(TotalLowInc)*100,1),
            Min = min(meanNightLST,na.rm=TRUE),
            Mean = mean(meanNightLST,na.rm=TRUE),
            Max = max(meanNightLST,na.rm=TRUE)) %>% 
  rename("City/Town" = NAME) %>% 
  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),"°")) %>% 
  arrange(desc(`Low Income Persons`)) %>% 
  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","striped")) %>% 
  kableExtra::column_spec(2:4, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 4, "Night 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 <- vt_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.2 * xrange) # xmin - left
bbox_new[3] <- bbox_new[3] + (0.1 * xrange) # xmax - right
bbox_new[2] <- bbox_new[2] - (0.1 * 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_borders(lwd = 0.2, alpha = 0.8) + 
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ny_state_sf) + tm_borders(lwd = 0.2, alpha = 0.8) +
  tm_shape(vt_highways) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(vt_highways2nd) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.15) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.1,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/vt_LSTNightNoHSDip5.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/vt_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."}
LSTNightNoHSDip5 %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `No HS Dip` = sum(lthsE,na.rm = TRUE),
            `Pct of No HS Dip` = round(sum(lthsE)/max(TotalNoHS)*100,1),
            Min = min(meanNightLST,na.rm=TRUE),
            Mean = mean(meanNightLST,na.rm=TRUE),
            Max = max(meanNightLST,na.rm=TRUE)) %>% 
  rename("City/Town" = NAME) %>% 
  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),"°")) %>% 
  arrange(desc(`No HS Dip`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Adults without a High School Diploma in Highest Quintile of Nighttime Avg LST", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header","striped")) %>% 
  kableExtra::column_spec(2:4, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 4, "Night 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 <- vt_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.2 * xrange) # xmin - left
bbox_new[3] <- bbox_new[3] + (0.1 * xrange) # xmax - right
bbox_new[2] <- bbox_new[2] - (0.1 * 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_borders(lwd = 0.2, alpha = 0.8) + 
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ny_state_sf) + tm_borders(lwd = 0.2, alpha = 0.8) +
  tm_shape(vt_highways) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(vt_highways2nd) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.15) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.1,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/vt_LSTNightOver645.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/vt_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."}
LSTNightOver645 %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Over 64` = sum(over64E,na.rm = TRUE),
            `Pct of Over 64` = round(sum(over64E)/max(TotalOver64)*100,1),
            Min = min(meanNightLST,na.rm=TRUE),
            Mean = mean(meanNightLST,na.rm=TRUE),
            Max = max(meanNightLST,na.rm=TRUE)) %>% 
  rename("City/Town" = NAME) %>% 
  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),"°")) %>% 
  arrange(desc(`Over 64`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Persons Over 64 in Highest Quintile of Nighttime Avg LST", align = "r") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header","striped")) %>% 
  kableExtra::column_spec(2:4, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 4, "Night 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 <- vt_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.2 * xrange) # xmin - left
bbox_new[3] <- bbox_new[3] + (0.1 * xrange) # xmax - right
bbox_new[2] <- bbox_new[2] - (0.1 * 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_borders(lwd = 0.2, alpha = 0.8) + 
  tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
  tm_shape(ny_state_sf) + tm_borders(lwd = 0.2, alpha = 0.8) +
  tm_shape(vt_highways) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(vt_highways2nd) + tm_lines(col = "seashell4", lwd = 1.5) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.15) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.15) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.1,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/vt_LSTNightUnder55.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/LST/vt_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."}
LSTNightUnder55 %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Under 5` = sum(under5E,na.rm = TRUE),
            `Pct of Under 5` = round(sum(under5E)/max(TotalUnder5)*100,1),
            Min = min(meanNightLST,na.rm=TRUE),
            Mean = mean(meanNightLST,na.rm=TRUE),
            Max = max(meanNightLST,na.rm=TRUE)) %>% 
  rename("City/Town" = NAME) %>% 
  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),"°")) %>% 
  arrange(desc(`Under 5`)) %>% 
  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","striped")) %>% 
  kableExtra::column_spec(2:4, width = "3cm") %>%
  kableExtra::add_header_above(c(" " = 4, "Night LST (°F)" = 3))
```



## Correlations {-}

```{r cormatrixLSTpopsVT, echo=FALSE, message=FALSE, warning=FALSE, fig.align = "center", fig.cap="Spearman's correlation matrix of temperature and the proportions of priority populations by Census Block Group."}
# Create a correlation matrix of PM25 and priority populations
corrLSTpops <- vt_blkgrpLST_sf %>% 
  as.data.frame() %>% 
  dplyr::transmute(`Day-Night Avg LST` = meanAvgLST,
                   `Day LST` = meanDayLST,
                   `Night LST` = meanNightLST, 
                `People of Color` = 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 Vermont", 
           legend.title = "Spearman's \nCorrelation\nCoefficient",
           lab_size = 2)
```



```{r cormatrixLSTAvgVTpollutants, 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 priority populations
corrLSTemissions <- vt_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 Vermont",
           legend.title = "Spearman's \nCorrelation\nCoefficient", 
           lab_size = 2)
```

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


```