Transport_VermontPDF.Rmd

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

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

library(tidyverse)
library(sf)
library(tmap)
library(tidytransit)
library(sp)
library(tigris)
options(tigris_use_cache = TRUE, tigris_class = "sf")
library(lubridate)
library(maptools)
library(spdep)
library(lwgeom)
```

\pagebreak

# Analysis of adequacy of public transit, walkability, and transportation cost burden in Vermont

This is an analysis of access to public transit and walkability in Vermont.

```{r data, include=FALSE}
### Analysis of Vermont transportation options
load("DATA/ne_layers.rds")

# extract census units for state
vt_blkgrps_sf <- ne_blkgrp_sf %>% 
  filter(STATE == "Vermont") %>% 
  select(GEOID, NAME,STATE, bg_area_m2, totalpopE, minorityE, under5E, under18E, over64E, householdsE, eng_hhE, eng_limitE, age25upE, lthsE, povknownE, num2povE) %>% 
  st_transform(., crs = 2852) %>% 
  filter(!st_is_empty(.)) %>% 
  mutate(PopAcre = totalpopE/(bg_area_m2*0.000247105))

vt_tracts_sf <- ne_tracts_sf %>% 
  filter(STATE == "Vermont") %>% 
  select(GEOID, NAME, totalpopE, STATE, Over18E, disabledOver18E, totalHHE, HHnoCarE) %>% 
  st_transform(., crs = 2852) %>% 
  filter(!st_is_empty(.))

# clean up
rm(ne_blkgrp_sf,ne_tracts_sf,ne_blkgrp_sf90)


## Read in VT transit data from ... GTFS since state-provided GIS files are simply drawn from GTFS but are missing variables. GTFS from https://vermont-gtfs.org/
# # Download GTFS data for Advance Transit
# download.file(url = "http://data.trilliumtransit.com/gtfs/advancetransit-vt-us/advancetransit-vt-us.zip", destfile = "DATA/transport/VT/VT_gtfs_files/AdvanceTransit.zip")
# # Download GTFS data for Green Mountain Transit
# download.file(url = "http://data.trilliumtransit.com/gtfs/ccta-vt-us/ccta-vt-us.zip", destfile = "DATA/transport/VT/VT_gtfs_files/GreenMtnTransit.zip")
# # Download GTFS data for Current Transit
# download.file(url = "http://data.trilliumtransit.com/gtfs/crtransit-vt-us/crtransit-vt-us.zip", destfile = "DATA/transport/VT/VT_gtfs_files/CurrentTransit.zip")
# # Download GTFS data for DVTA Moover
# download.file(url = "http://data.trilliumtransit.com/gtfs/dvtamoover-vt-us/dvtamoover-vt-us.zip", destfile = "DATA/transport/VT/VT_gtfs_files/DVTAMooverTransit.zip")
# # Download GTFS data for Green Mountain Community Network
# download.file(url = "http://data.trilliumtransit.com/gtfs/greenmtncn-vt-us/greenmtncn-vt-us.zip", destfile = "DATA/transport/VT/VT_gtfs_files/GreenMtnComNetworkTransit.zip")
# # Download GTFS data for Lake Champlain Ferries
# download.file(url = "http://data.trilliumtransit.com/gtfs/lctferries-vt-us/lctferries-vt-us.zip", destfile = "DATA/transport/VT/VT_gtfs_files/LakeChamplainFerries.zip")
# # Download GTFS data for Rural Community Transportation
# download.file(url = "http://data.trilliumtransit.com/gtfs/ruralcommunity-vt-us/ruralcommunity-vt-us.zip", destfile = "DATA/transport/VT/VT_gtfs_files/RuralComTransport.zip")
# # Download GTFS data for The Bus - Marble Valley Regional Transit District: Marble Valley
# download.file(url = "http://data.trilliumtransit.com/gtfs/thebus-vt-us/thebus-vt-us.zip", destfile = "DATA/transport/VT/VT_gtfs_files/MarbleValleyTransit.zip")
# # Download GTFS data for Stagecoach Transportation Services, Tri-Valley Transit
# download.file(url = "http://data.trilliumtransit.com/gtfs/trivalleytransit-vt-us/trivalleytransit-vt-us.zip", destfile = "DATA/transport/VT/VT_gtfs_files/TriValleyTransit.zip")
# # Download GTFS data for Vermont Translines
# download.file(url = "http://data.trilliumtransit.com/gtfs/vttranslines-vt-us/vttranslines-vt-us.zip", destfile = "DATA/transport/VT/VT_gtfs_files/VermontTranslines.zip")

# create list of downloaded GTFS files
fileList <- as.list(list.files("DATA/transport/VT/VT_gtfs_files/", pattern = ".zip$", full.names = TRUE))

# extract just file names without extension
fileListNames <- lapply(fileList, 
                        function(x) fs::path_ext_remove(fs::path_file(x)))

# name items in fileList
fileList <- setNames(fileList,fileListNames)

# read in each file to gtfs
VT_gtfs_list <- lapply(fileList, 
                         function(x) read_gtfs(x) %>% 
                           gtfs_as_sf(.))

# add a table to the feed that indicates which service_id runs on which date. This is later useful for linking dates and trips via service_id.
VT_gtfs_list <- lapply(VT_gtfs_list, set_date_service_table)

# To understand service patterns better we need information on weekdays and holidays. With a calendar table we know the weekday and possible holidays for each date.
# create a function to generate calendar and join to date_servicepattern_table
calendarFunction <- function(thing) {
  tibble(date = unique(thing$.$date_service_table$date)) %>% 
  mutate(
    weekday = (function(date) {
      c("Sunday", "Monday", "Tuesday", 
        "Wednesday", "Thursday", "Friday", 
        "Saturday")[as.POSIXlt(date)$wday + 1]
    })(date)
  ) %>% 
    left_join(x = thing$.$date_servicepattern_table,y = ., by = "date", copy=TRUE)
}

# To analyse on which dates trips run and to group similar services we use service patterns. Such a pattern simply lists all dates a trip runs on. To handle these patterns we create a servicepattern_id using a hash function. 
VT_gtfs_list <- lapply(VT_gtfs_list, set_servicepattern)

# Our gtfs feed now contains the data frame service_pattern which links each servicepattern_id to an existing service_id (and by extension trip_id).
# Understand patterns of service by visualising the data.
# Iterate over list to generate service pattern tables
date_servicepattern_table_list <- lapply(X = VT_gtfs_list, calendarFunction)

# This part has to be done manually for each item. sigh.
# create function to generate service pattern graph for each item
ggplotPatterns <- function (x) {
  ggplot(x) + theme_bw() +
  geom_point(aes(x = date, y = servicepattern_id, color = weekday), size = 1) +
  scale_x_date(breaks = scales::date_breaks("1 month")) +
  theme(legend.position = "bottom")
}
# iterate through list and create plot for each. look at each and identify M-F servcie pattern ids.
ggplotServicePatterns <- lapply(date_servicepattern_table_list, ggplotPatterns)

AdvanceTransit_spids <- c("s_bd3d812","s_4834b03","s_3fbb598")
CurrentTransit_spids <- c("s_b8c28a8","s_6ec375d","s_4055657")
DVTAMooverTransit_spids <- c("s_aa691be","s_273dfe9","s_232a0a5","s_2141cf9")
GreenMtnComNetworkTransit_spids <- c("s_e633768", "s_b86cbed", "s_b134a45", "s_39b322f")
GreenMtnTransit_spids <- c("s_b258906","s_8cd8e19","s_199639b","s_121b32f")
LakeChamplainFerries_spids <- c("s_c03da1d","s_6cb5291")
MarbleValleyTransit_spids <- c("s_dec72ff","s_cc4d0f9","s_c5c8bb8","s_c045849", "s_06f20d5", "s_0418370")
RuralComTransport_spids <- c("s_f30b07e","s_a5157f4","s_60b8468")
TriValleyTransit_spids <- c("s_e5732f7","s_a2fbba2","s_1c55bc1")
VermontTranslines_spids <- "s_46f7ebd"


# create the list. make sure to remove spids_list before running, otherwise it gets included in the list!
spids_list <- mget(ls(pattern = "_spids"))

# iterate over list items to extract service ids
service_ids_list <- map2(VT_gtfs_list,spids_list, 
                         function(x,y) x$.$service_pattern %>%
                           filter(servicepattern_id %in% y) %>% 
                           pull(service_id))

# now that we’ve identified the set of service_id’s that refer to all weekday trips, we can summarize service between 6 am and 9 pm for bus service on weekdays.
# iterate over list items. skip Shore Line East because it's missing `direction-id` column
daily_stop_freq_list <- map2(VT_gtfs_list,service_ids_list,
                              function(x,y) get_stop_frequency(x, start_hour = 6, end_hour = 21, service_ids = y))

# # Convert stops to points for mapping
# Convert stops to points for mapping for each item in list
vt_transit_stops_sf_list <- lapply(X = VT_gtfs_list, function(x) stops_as_sf(x$stops))

# Join headway frequencies to stops and route descriptions
# first, join headway frequencies
vt_transit_stops_sf_list <- map2(vt_transit_stops_sf_list,daily_stop_freq_list,
                                 function(x,y) x %>% 
                                   inner_join(y, by = "stop_id"))
# next, join route descriptions
vt_transit_stops_sf_list <- map2(vt_transit_stops_sf_list,VT_gtfs_list,
                                 function(x,y) x %>% 
                                   inner_join(y$routes, by = "route_id"))

# use the get_route_frequency function to summarize transit service by route, for the same time period.
# iterate over list items
daily_route_freq_list <- map2(VT_gtfs_list,service_ids_list, 
                              function(x,y) 
                                get_route_frequency(x, service_ids = y))

# Join the route frequencies to geometry for mapping
# iterate over list to extract route geometry. need to skip Hartford Line and Shoreline East because they're missing shapes.
vt_routes_sf_list <- map2(VT_gtfs_list,service_ids_list, 
                          function(x,y) get_route_geometry(x, service_ids = y))


# join calculated frequencies to geometry and route descriptions
# first, join calculated frequencies to geometry. Again, skip Hartford Line and Shoreline East.
vt_routes_sf_list <- map2(vt_routes_sf_list, daily_route_freq_list,
                          function(x,y) inner_join(x,y, by = "route_id"))
# next, join route descriptions to allow filtering
vt_routes_sf_list <- map2(vt_routes_sf_list,VT_gtfs_list, 
                          function(x,y) 
                            inner_join(x,y$routes, by = "route_id"))

# Extract bus routes and common set of variables and combine into one layer
vt_routes_sf_Bus <- lapply(vt_routes_sf_list, function(x) 
  filter(x, route_type == 3) %>%
    select(route_id:mean_headways,agency_id,route_short_name,route_long_name)) %>% 
  do.call(rbind, .) %>% 
  mutate(agency = case_when(
    agency_id == "738" ~ "Advance Transit",
    agency_id == "232" ~ "The Current",
    agency_id == "218" ~ "MOOver - DVTA",
    agency_id == "219" ~ "Green Mountain Community Network",
    agency_id == "461" ~ "Green Mountain Transit",
    agency_id == "885" ~ "Special Services Transportation",
    agency_id == "789" ~ "Lake Champlain Ferries",
    agency_id == "220" ~ "Marble Valley Regional Transit",
    agency_id == "221" ~ "Rural Community Transportation",
    agency_id == "986" ~ "Tri Valley Transit",
    agency_id == "228" ~ "Vermont Translines"
  )) %>% 
  st_transform(., crs = 2852) %>% 
  filter(!st_is_empty(.)) %>%
  st_make_valid()

# create a simplified version for faster mapping
# vt_routes_sf_Bus_simple <- vt_routes_sf_Bus %>% 
#   st_simplify(., dTolerance = 100) %>% 
#   st_make_valid()

vt_routes_sf_Bus_simple <- vt_routes_sf_Bus %>% 
  # unique(.) %>% 
  group_by(route_short_name) %>% 
  summarize(agency = first(agency),agency_id = first(agency_id),mean_headways = mean(mean_headways,na.rm = TRUE)) %>% 
  st_make_valid()

# extract bus stops and merge to one layer
vt_stops_sf_Bus <- lapply(vt_transit_stops_sf_list, function(x)
  filter(x, route_type == 3) %>% 
    select(agency_id,stop_id,route_type,headway)) %>% 
  do.call(rbind, .) %>% 
  mutate(agency = case_when(
    agency_id == "738" ~ "Advance Transit",
    agency_id == "232" ~ "The Current",
    agency_id == "218" ~ "MOOver - DVTA",
    agency_id == "219" ~ "Green Mountain Community Network",
    agency_id == "461" ~ "Green Mountain Transit",
    agency_id == "885" ~ "Special Services Transportation",
    agency_id == "789" ~ "Lake Champlain Ferries",
    agency_id == "220" ~ "Marble Valley Regional Transit",
    agency_id == "221" ~ "Rural Community Transportation",
    agency_id == "986" ~ "Tri Valley Transit",
    agency_id == "228" ~ "Vermont Translines"
  )) %>% 
  st_transform(., crs = 2852) %>% 
  filter(!st_is_empty(.)) %>%
  st_make_valid()

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

# read in park-and-ride lots from Vermont Open Geodata Portal at https://geodata.vermont.gov/ . Keep only park-and-rides that are serviced by public transit. 
vt_parknride_sf <- st_read(dsn = "DATA/transport/VT/VT_Park_and_Rides",
                           layer = "VT_Park_and_Rides") %>% 
  filter(TRANSIT != "NONE") %>% 
  st_transform(., crs = 2852) %>% 
  st_make_valid()


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

##### Return to working in R ######
# # read in processed vt_blkgrps and vt_tracts
# st_layers(dsn = "DATA/FEMA/VT")
vt_blkgrps_developed <- st_read(dsn = "DATA/FEMA/VT",
                                layer = "vt_blkgrps_developed") 

vt_blkgrps_developed <- vt_blkgrps_sf %>% 
  as.data.frame() %>% 
  dplyr::select(-geometry) %>%
  left_join(vt_blkgrps_developed, ., by = "GEOID") %>% 
  st_transform(., crs = 2852) %>% 
  mutate(NewArea = st_area(.)) %>%
  st_make_valid() %>% 
  mutate(Proportion = as.numeric(NewArea/OldArea),
         NewPop = totalpopE*Proportion,
         NewMinority = minorityE*Proportion,
         NewUnder5 = under5E*Proportion,
         NewOver64 = over64E*Proportion,
         NewUnder18 = under18E*Proportion,
         NewEng_limit = eng_limitE*Proportion,
         NewPov = num2povE*Proportion,
         NewLths = lthsE*Proportion,
         NewPopAcre = NewPop/(NewArea*0.000247105))

vt_tracts_developed <- st_read(dsn = "DATA/FEMA/VT",
                                layer = "vt_tracts_developed") %>% 
  left_join(., as.data.frame(vt_tracts_sf), by = "GEOID") %>% 
  st_transform(., crs = 2852) %>% 
  mutate(NewArea = st_area(.)) %>%
  st_make_valid() %>% 
  mutate(Proportion = as.numeric(NewArea/OldArea),
         NewDisabled = disabledOver18E*Proportion,
         NewOver18 = Over18E*Proportion,
         NewHHNoCar = HHnoCarE*Proportion)

# create a simplified version for faster mapping
vt_blkgrps_developed_simple <- st_simplify(vt_blkgrps_developed,
                                           dTolerance = 100) %>% 
  filter(!st_is_empty(.)) %>% 
  st_make_valid()
# repeat for tracts
vt_tracts_developed_simple <- st_simplify(vt_tracts_developed, 
                                          dTolerance = 100) %>% 
  filter(!st_is_empty(.)) %>% 
  st_make_valid()

# Create a 400m buffer around bus stops
vt_busBuff400m <- st_buffer(vt_stops_sf_Bus,dist = 400) %>% 
  st_union() %>% 
  st_as_sf()

# Create buffers around park-and-ride stops (4800m)
vt_prBuff4800m <- st_buffer(vt_parknride_sf, dist = 4800) %>% 
  st_union() %>% 
  st_as_sf()

vt_transitAllBuff <- rbind(vt_busBuff400m,
                           vt_prBuff4800m) %>% 
  st_union() %>% 
  st_as_sf() %>% 
  mutate(Sqm = st_area(.))

# Use areal interpolation to calculate priority populations within buffer of accessibility
vt_busBuff400m_sf <- vt_blkgrps_developed %>% 
  select(GEOID,NewArea:NewLths) %>% 
  mutate(OldArea = NewArea) %>% 
  st_intersection(vt_busBuff400m,.) %>% 
  mutate(NewArea = st_area(.),
         Proportion = NewArea/OldArea,
         NewPop = as.integer(NewPop*Proportion),
         NewMinority = as.integer(NewMinority*Proportion),
         NewUnder5 = as.integer(NewUnder5*Proportion),
         NewOver64 = as.integer(NewOver64*Proportion),
         NewUnder18 = as.integer(NewUnder18*Proportion),
         NewEng_limit = as.integer(NewEng_limit*Proportion),
         NewPov = as.integer(NewPov*Proportion),
         NewLths = as.integer(NewLths*Proportion))

vt_prBuff4800m_sf <- vt_blkgrps_developed %>% 
  select(GEOID,NewArea:NewLths) %>% 
  mutate(OldArea = NewArea) %>% 
  st_intersection(vt_prBuff4800m,.) %>% 
  mutate(NewArea = st_area(.),
         Proportion = NewArea/OldArea,
         NewPop = as.integer(NewPop*Proportion),
         NewMinority = as.integer(NewMinority*Proportion),
         NewUnder5 = as.integer(NewUnder5*Proportion),
         NewOver64 = as.integer(NewOver64*Proportion),
         NewUnder18 = as.integer(NewUnder18*Proportion),
         NewEng_limit = as.integer(NewEng_limit*Proportion),
         NewPov = as.integer(NewPov*Proportion),
         NewLths = as.integer(NewLths*Proportion))

vt_TransitAllBuff_sf <- vt_blkgrps_developed %>% 
  select(GEOID,NewArea:NewLths) %>% 
  mutate(OldArea = NewArea) %>% 
  st_intersection(vt_transitAllBuff,.) %>% 
  mutate(NewArea = st_area(.),
         Proportion = NewArea/OldArea,
         NewPop = as.integer(NewPop*Proportion),
         NewMinority = as.integer(NewMinority*Proportion),
         NewUnder5 = as.integer(NewUnder5*Proportion),
         NewOver64 = as.integer(NewOver64*Proportion),
         NewUnder18 = as.integer(NewUnder18*Proportion),
         NewEng_limit = as.integer(NewEng_limit*Proportion),
         NewPov = as.integer(NewPov*Proportion),
         NewLths = as.integer(NewLths*Proportion))

# Repeat for tracts
vt_busBuff400mTracts_sf <- vt_tracts_developed %>% 
  select(GEOID,NewDisabled:NewHHNoCar) %>% 
  mutate(OldArea = st_area(.)) %>% 
  st_intersection(vt_busBuff400m,.) %>% 
  mutate(NewArea = st_area(.),
         Proportion = NewArea/OldArea,
         NewDisabled = as.integer(NewDisabled*Proportion),
         NewNoCar = as.integer(NewHHNoCar*Proportion))

vt_prBuff4800mTracts_sf <- vt_tracts_developed %>% 
  select(GEOID,NewDisabled:NewHHNoCar) %>% 
  mutate(OldArea = st_area(.)) %>% 
  st_intersection(vt_prBuff4800m,.) %>% 
  mutate(NewArea = st_area(.),
         Proportion = NewArea/OldArea,
         NewDisabled = as.integer(NewDisabled*Proportion),
         NewNoCar = as.integer(NewHHNoCar*Proportion))

vt_TransitAllBuffTracts_sf <- vt_tracts_developed %>% 
  select(GEOID,NewDisabled:NewHHNoCar) %>% 
  mutate(OldArea = st_area(.)) %>% 
  st_intersection(vt_transitAllBuff,.) %>% 
  mutate(NewArea = st_area(.),
         Proportion = NewArea/OldArea,
         NewDisabled = as.integer(NewDisabled*Proportion),
         NewNoCar = as.integer(NewHHNoCar*Proportion))

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

vt_prBuff4800m_df <- vt_prBuff4800m_sf  %>% 
  as.data.frame() %>% 
  summarize(`Total Pop` = sum(NewPop),
            Minority = sum(NewMinority),
            `Under 5` = sum(NewUnder5),
            `Over 64` = sum(NewOver64),
            `Under 18` = sum(NewUnder18),
            `Limited English HH` = sum(NewEng_limit),
            `Low Income` = sum(NewPov),
            `No HS Dip` = sum(NewLths)) %>% 
  gather(key = Group, value = ParknRidePop)

vt_TransitAllBuff_df <- vt_TransitAllBuff_sf %>% 
  as.data.frame() %>% 
  summarize(`Total Pop` = sum(NewPop),
            Minority = sum(NewMinority),
            `Under 5` = sum(NewUnder5),
            `Over 64` = sum(NewOver64),
            `Under 18` = sum(NewUnder18),
            `Limited English HH` = sum(NewEng_limit),
            `Low Income` = sum(NewPov),
            `No HS Dip` = sum(NewLths)) %>% 
  gather(key = Group, value = TransitPop)

# Compute total tract populations within bus stop buffer
vt_busBuff400mTracts_df <- vt_busBuff400mTracts_sf %>% 
  as.data.frame() %>% 
  summarize(`Disabled` = sum(NewDisabled),
            `No Car HH` = sum(NewNoCar)) %>% 
  gather(key = Group, value = BusPop)

vt_prBuff4800mTracts_df <- vt_prBuff4800mTracts_sf %>% 
  as.data.frame() %>% 
  summarize(`Disabled` = sum(NewDisabled),
            `No Car HH` = sum(NewNoCar)) %>% 
  gather(key = Group, value = ParknRidePop)

vt_TransitAllBuffTracts_df <- vt_TransitAllBuffTracts_sf %>% 
  as.data.frame() %>% 
  summarize(`Disabled` = sum(NewDisabled),
            `No Car HH` = sum(NewNoCar)) %>% 
  gather(key = Group, value = TransitPop)

# Compute total tract populations within the state for same groups
vt_tract_pops_df <- vt_tracts_sf %>% 
  as.data.frame() %>% 
  summarize(`Disabled` = sum(disabledOver18E),
            `No Car HH` = sum(HHnoCarE)) %>% 
  gather(key = Group, value = VTPop) %>% 
  left_join(.,vt_busBuff400mTracts_df, by = "Group") %>% 
  left_join(.,vt_prBuff4800mTracts_df, by = "Group") %>% 
  left_join(.,vt_TransitAllBuffTracts_df, by = "Group")

# Compute populations for state, join with buffer pops
vt_transitAccessPops_df <- vt_blkgrps_sf %>% 
  as.data.frame() %>%  
  summarize(`Total Pop` = sum(totalpopE),
            Minority = sum(minorityE),
            `Under 5` = sum(under5E),
            `Over 64` = sum(over64E),
            `Under 18` = sum(under18E),
            `Limited English HH` = sum(eng_limitE),
            `Low Income` = sum(num2povE),
            `No HS Dip` = sum(lthsE)) %>% 
  gather(key = Group, value = VTPop) %>% 
  left_join(., vt_busBuff400m_df, by = "Group") %>% 
  left_join(., vt_prBuff4800m_df, by = "Group") %>% 
  left_join(., vt_TransitAllBuff_df, by = "Group") %>% 
  rbind(.,vt_tract_pops_df) %>% 
  mutate(PctBus = BusPop/VTPop*100,
         PctParknRide = ParknRidePop/VTPop*100,
         PctAllTransit = TransitPop/VTPop*100)

# create point layer of towns for context
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")

# # recode routes to match VTPTA categories
# vt_busroutes_sf <- vt_busroutes_sf %>% 
#   mutate(type = case_when(
#     ROUTE %in% c(3,4,6,13,17,18,19,21,22,29,30,32,33,34,35,40,49,51,55,57,58,63,64,67,72,73,75,76,78,80,87) ~ "Urban Service",
#     ROUTE == 11 ~ "Rapid Line",
#     ROUTE %in% c(8,9,10,12,14,54,59,60,61,62,65,66,95) ~ "Express Service",
#     ROUTE %in% c(1,20,27,28,31,50,56,71,92) ~ "Key Corridor",
#     ROUTE %in% c(16,203,204,210,231,242,281,282,301) ~ "Flex Service"))

# Calculate totals of priority populations for use in computing percentages in later tables
# Calculate tract populations by state
statepopsdf_tracts <- vt_tracts_sf %>% 
  as.data.frame() %>% 
  group_by(STATE) %>% 
  summarize(TotalDisabled = sum(disabledOver18E, na.rm = TRUE),
            TotalOver18 = sum(Over18E, na.rm = TRUE),
            TotalNoCar = sum(HHnoCarE, na.rm = TRUE))

statepopsdf <- vt_blkgrps_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),
            TotalUnder18 = sum(under18E, na.rm = TRUE)) %>% 
  left_join(., statepopsdf_tracts,by="STATE")

# 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_tracts <- vt_tracts_sf %>% 
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  group_by(NAME) %>% 
  summarize(TotalDisabled = sum(disabledOver18E, na.rm = TRUE),
            TotalOver18 = sum(Over18E, na.rm = TRUE),
            TotalNoCar = sum(HHnoCarE, na.rm = TRUE))

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

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


# Identify areas >= 1 mile from bus stops
# Create a 1 mile buffer around bus stops
vt_busBuff1mile <- st_buffer(vt_stops_sf_Bus,dist = 1609.34) %>% 
  st_union() %>% 
  st_as_sf()
# identify block groups that do not intersect with 1 mile from bus stop
vt_blkgrps_sf_noBus <- vt_blkgrps_developed %>% 
  filter(st_disjoint(.,vt_busBuff1mile,sparse = FALSE))
# Identify tracts that do not intersect with buffer
vt_tracts_sf_noBus <- vt_tracts_developed %>% 
  filter(st_disjoint(.,vt_busBuff1mile,sparse = FALSE))
# create simplified versions for faster mapping
vt_blkgrps_sf_noBus_simple <- vt_blkgrps_sf_noBus %>% 
  st_simplify(., dTolerance = 100) %>% 
  filter(!st_is_empty(.)) %>% 
  st_make_valid()
vt_tracts_sf_noBus_simple <- vt_tracts_sf_noBus %>% 
  st_simplify(., dTolerance = 100) %>% 
  filter(!st_is_empty(.)) %>% 
  st_make_valid()

# Identify areas >= 10 miles from park-and-ride stops
vt_prBuff10mile <- vt_parknride_sf %>% 
  st_buffer(., dist = (10*1609.34)) %>% 
  st_union() %>% 
  st_as_sf()
# identify block groups that do not intersect with 10 mile buffer
vt_blkgrps_sf_noPR <- vt_blkgrps_developed %>% 
  filter(st_disjoint(.,vt_prBuff10mile,sparse = FALSE))
# create simplified version of block groups for faster mapping
vt_blkgrps_sf_noPR_simple <- vt_blkgrps_sf_noPR %>% 
  select(NewPop:NewLths) %>% 
  st_simplify(., dTolerance = 100) %>% 
  filter(!st_is_empty(.)) %>% 
  st_make_valid()

# identify tracts that do not intersect with 10 mile buffer
vt_tracts_sf_noPR <- vt_tracts_developed %>% 
  filter(st_disjoint(.,vt_prBuff10mile,sparse = FALSE))
# create simplified version of tracts for faster mapping
vt_tracts_sf_noPR_simple <- vt_tracts_sf_noPR %>% 
  select(NewDisabled:NewHHNoCar) %>% 
  st_simplify(., dTolerance = 100) %>% 
  filter(!st_is_empty(.)) %>% 
  st_make_valid()

# Identify areas with no access to any type of transit stop. Union all buffers to create one buffer and clip to state boundary. 
vt_transitBuffAll <- rbind(vt_busBuff1mile,
                             vt_prBuff10mile) %>% 
  st_union() %>% 
  st_as_sf() %>% 
  tmaptools::crop_shape(.,vt_state_sf_cb,polygon = TRUE)

# identify block groups that do not intersect with any transit buffer
vt_blkgrps_sf_noTransit <- vt_blkgrps_developed %>% 
  filter(st_disjoint(.,vt_transitBuffAll,sparse = FALSE))
# identify tracts that do not intersect with 5 mile buffer
vt_tracts_sf_noTransit <- vt_tracts_developed %>% 
  filter(st_disjoint(.,vt_transitBuffAll,sparse = FALSE))
# create simplified versions for faster mapping
vt_blkgrps_sf_noTransit_simple <- vt_blkgrps_sf_noTransit %>% 
  st_simplify(., dTolerance = 100) %>% 
  filter(!st_is_empty(.)) %>% 
  st_make_valid()
vt_tracts_sf_noTransit_simple <- vt_tracts_sf_noTransit %>% 
  st_simplify(., dTolerance = 100) %>% 
  filter(!st_is_empty(.)) %>% 
  st_make_valid()

# write out to file for later analysis
save(vt_blkgrps_sf_noTransit, vt_tracts_sf_noTransit, 
     vt_blkgrps_sf_noBus, vt_transitAccessPops_df, vt_busBuff400m_sf,
     file = "DATA/transport/VT/noTransit.Rds")
```

## Access to Public Transit

This analysis seeks to identify communities that are underserved by access to public transit. Access to transit is measured in terms of distance to transit boarding stops and in terms of frequency of service. 

Vermont is served by 10 local bus agencies or transit districts, two inter-state rail lines operated by Amtrak, local and inter-state ferry services on Lake Champlain, and several private inter-city bus services. In addition, there are over 100 park-and-ride lots throughout the state, variously owned by the state, municipalities, and by public transit agencies. The analysis presented here only considers public transit services for which publicly available geospatial data is available, which includes the `r length(unique(vt_routes_sf_Bus$agency))`
regional public transit agencies and the `r nrow(vt_parknride_sf)` park-and-ride lots that are directly serviced by these transit agencies. Amtrak train service, private inter-city buses, ferry services, water taxis, and seasonal transportation services are not considered in this analysis. Figure \@ref(fig:transitMap) shows `r formatC(nrow(vt_routes_sf_Bus),big.mark = ",")` fixed public transit routes across the state along with population density.

```{r transitMap, fig.align = "center", fig.cap="Transit service and population density in Vermont."}
# bbox_new <- st_bbox(ct_blkgrps_sf) # 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] - (2.2 * xrange) # xmin - left
# # bbox_new[3] <- bbox_new[3] + (2 * xrange) # xmax - right
# bbox_new[2] <- bbox_new[2] - (0.2 * yrange) # ymin - bottom
# # bbox_new[4] <- bbox_new[4] + (2 * yrange) # ymax - top
# 
# bbox_new <- bbox_new %>%  # take the bounding box ...
#   st_as_sfc() # ... and make it a sf polygon

m <- tm_shape(vt_towns_sf, unit = "mi") +
  tm_borders(col = "gray", lwd = 0.2) +
  tm_shape(vt_blkgrps_developed_simple) +
  tm_fill(col = "NewPopAcre", style = "quantile",
          colorNA = "white", colorNULL = "gray", title = "Pop per Acre",
          legend.format = list(fun = function(x) formatC(x, digits = 1, format = "f")), palette = "YlGn") +
  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(vt_routes_sf_Bus_simple) + tm_lines(col = "coral", lwd = 1) +
  tm_shape(vt_parknride_sf) + tm_dots(col = "red", size = 0.1) +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5,
               position = c(0.6,0.005)) +
  tm_add_legend(type = "line", labels = "Bus Routes", col = "coral") +
  tm_add_legend(type = "symbol", labels = "Park-and-Ride Lots", 
                col = "red", size = 0.2, border.lwd = 0) +
  tm_layout(title = "Vermont\nTransit Routes\nand Population\nDensity",
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_transit.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_transit.png")
```

Approximately `r paste0(round(vt_transitAccessPops_df %>% filter(Group == "Total Pop") %>% dplyr::select(PctAllTransit) %>% pull(),1),"%")` of the state population lives within reasonable walking distance of a bus stop (defined as approximately 400 meters or 1/4 mile) or 3 miles of a park-and-ride lot that is serviced by transit. 

However, this access varies by population group. Figure \@ref(fig:loliBusDist) below compare the percentages of different population groups living within reasonable walking distance of bus stops (1/4 mile) or within 3 miles of park-and-ride lots across the state. More than half of the priority populations, including transit-dependent groups, have access similar to or better than the general population. However, children under 5, people under 18 and over 64, and adults without a high school diploma have consistently less access to transit on average compared to the general population.

```{r loliBusDist, echo=FALSE, message=FALSE, warning=FALSE, fig.align = "center", fig.cap="Percentage of population groups within 1/4 mile (400 meters) of a bus stop or 3 miles (4800m) of park-and-ride lot across Vermont"}
vt_transitAccessPops_df %>% 
  mutate(Group = recode(Group, "Minority" = "People of Color",
                        "No HS Dip" = "No HS Diploma")) %>%
  ggplot(aes(x = reorder(Group,-PctAllTransit), 
             y = PctAllTransit)) +
  geom_segment(aes(x = reorder(Group,-PctAllTransit), xend = reorder(Group,-PctAllTransit),
                   y = vt_transitAccessPops_df[1,8], yend = PctAllTransit), 
               color = "skyblue") +
  geom_point(color = "blue", size = 4, alpha = 0.8) +
  coord_flip() + xlab("") + ylab("") + ggtitle("Vermont Populations within Walking Distance (400m)\nof Bus Stops or 4800m of Park-and-Ride Lots") + theme_light() +
  theme(panel.grid.major.y = element_blank(),
        panel.border = element_blank(),
        axis.ticks.y = element_blank()) +
  geom_text(aes(x = Group, y = PctAllTransit + 0.2 * sign(PctAllTransit), 
                label = paste0(round(PctAllTransit,1),"%")), 
            hjust = 0, vjust = -.4, size = 3,
            color=rgb(100,100,100, maxColorValue=255)) +
  scale_y_continuous(labels = function(x) paste0(x, "%")) +
  geom_hline(yintercept = vt_transitAccessPops_df[1,8], linetype = "dashed") +
  geom_text(aes(x = "Total Pop", y = 43, label = "Above state\naverage"),
            color = "gray48") +
  geom_text(aes(x = "Total Pop", y = 36, label = "Below state\naverage"),
            color = "gray48") +
  expand_limits(y = c(35,51))

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

Specific values of the same can be seen in Table \@ref(tab:statsTransitDist) below.

```{r statsTransitDist, echo=FALSE, message=FALSE, warning=FALSE, fig.align = "center", fig.cap="Vermont populations living within 1/4 mile (400 meters) of a bus stop or 3 miles (4800m) of a park-and-ride lot. Based on ACS 2018 Block Group data."}
# See table of distance by group
vt_transitAccessPops_df %>% 
  mutate(Group = recode(Group, "Minority" = "People of Color",
                        "No HS Dip" = "No HS Diploma")) %>%
  mutate(PctBus = paste0(round(PctBus,1),"%"),
         PctParknRide = paste0(round(PctParknRide,1),"%"),
         PctAllTransit = paste0(round(PctAllTransit,1),"%")) %>% 
  kableExtra::kable(longtable = T, booktabs = T, 
                    format.args = list(big.mark = ','), align = "r",
                    caption = "Populations with Access to Transit", 
                    col.names = c("Group","VT Pop","Bus","ParknRide","Any Transit","Bus","ParknRide","Any Transit")) %>% 
  kableExtra::column_spec(3:8, width = "1.7cm") %>%
  kableExtra::add_header_above(c(" " = 2, "Total Pop with Access" = 3, "Pct Pop with Access" = 3)) %>% 
  kableExtra::footnote(general = "Access is 1/4 mile (400m) from bus stop or 3 miles (4800m) from park-and-ride lots.") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
  
```

The population percentages above provide a glimpse of average access to public transit across the state for different population groups. However, access also varies for specific communities. Indeed, it is important to consider vulnerable populations or transportation-limited populations that do not have reasonable access to transit. 

```{r TransitNoAccess, include=FALSE}
TransitNoAccess1 <- vt_blkgrps_sf_noTransit %>%
  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(NewOver64,na.rm = TRUE),
            `Pct of Under 5` = round(sum(NewUnder5)/max(TotalUnder5)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Pct of Under 5` > 0) %>% 
  mutate(`Pct of Under 5` = if_else(`Pct of Under 5` > 100, 100, `Pct of Under 5`)) %>% 
  mutate(`Pct of Under 5` = paste0(`Pct of Under 5`,"%")) %>% 
  arrange(desc(`Under 5`))

TransitNoAccess2 <- vt_blkgrps_sf_noTransit %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Under 18` = sum(NewOver64,na.rm = TRUE),
            `Pct of Under 18` = round(sum(NewUnder18)/max(TotalUnder18)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Pct of Under 18` > 0) %>% 
  mutate(`Pct of Under 18` = if_else(`Pct of Under 18` > 100, 100, `Pct of Under 18`)) %>% 
  mutate(`Pct of Under 18` = paste0(`Pct of Under 18`,"%")) %>% 
  arrange(desc(`Under 18`))
```

Figures \@ref(fig:mapTransitAccessU5) and \@ref(fig:mapTransitAccessU18) below highlight developed portions of Census Block Groups across Vermont without reasonable access to bus or park-and-ride lots, and the numbers of persons under 5 and under 18 who live there. Tables \@ref(tab:tabTransitAccessU5) and \@ref(tab:tabTransitAccessU18) show breakdowns of the same by municipality. For example, Table \@ref(tab:tabTransitAccessU5) shows that `r formatC(as.numeric(TransitNoAccess1[1,3]),big.mark=",")` children under 5 across `r TransitNoAccess1[1,2]` Block Groups in `r TransitNoAccess1[1,1]` resided one or more miles from the nearest bust stop. These people represented `r TransitNoAccess1[1,4]` of children under 5 in `r TransitNoAccess1[1,1]`. Similarly, Table \@ref(tab:tabTransitAccessU18) shows that `r formatC(as.numeric(TransitNoAccess2[1,3]),big.mark=",",digits = 5)` people under 18 across `r TransitNoAccess2[1,2]` Block Groups in `r TransitNoAccess2[1,1]` resided one or more miles from the nearest transit stop for any mode. These people represented `r TransitNoAccess2[1,4]` of people under 18 in `r TransitNoAccess2[1,1]`.
Maps and tables of other priority populations without access to the various modes of transit can be found in Appendix B.

```{r mapTransitAccessU5, fig.align = "center", fig.cap="Children under 5 without reasonable access to any form of transit by Census Block Group."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") +
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(vt_blkgrps_sf_noTransit_simple, unit = "mi") +
  tm_fill(col = "NewUnder5", style = "quantile",
          colorNA = "white", colorNULL = "gray", title = "Number of\nChildren Under 5",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd") +
  tm_shape(vt_tracts_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5,
               position = c(0.6,0.005)) +
  tm_layout(title = "Children Under 5\nwithout Reasonable\nAccess to any\nMode of Transit",
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_transitU5access.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_transitU5access.png")
```

```{r tabTransitAccessU5, fig.align = "center", fig.cap="Children under 5 in Census Block Groups without reasonable access to any mode of transit."}
vt_blkgrps_sf_noTransit %>%
  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(NewUnder5,na.rm = TRUE),
            `Pct of Under 5` = round(sum(NewUnder5)/max(TotalUnder5)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Pct of Under 5` > 0) %>% 
  mutate(`Pct of Under 5` = if_else(`Pct of Under 5` > 100,100,`Pct of Under 5`)) %>% 
  mutate(`Pct of Under 5` = paste0(`Pct of Under 5`,"%")) %>% 
  arrange(desc(`Under 5`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Children Under 5 without Access to any Mode of Transit", align = "r", col.names = c(names(.)[1:2],"Number of Under 5 in Census Block Groups","Pct of Under 5 in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped")) 
```

```{r mapTransitAccessU18, fig.align = "center", fig.cap="Persons under 18 without reasonable access to any form of transit by Census Block Group."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") +
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(vt_blkgrps_sf_noTransit_simple, unit = "mi") +
  tm_fill(col = "NewUnder18", style = "quantile",
          colorNA = "white", colorNULL = "gray", title = "Number of\nPersons Under 18",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd") +
  tm_shape(vt_tracts_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5,
               position = c(0.6,0.005)) +
  tm_layout(title = "Persons Under 18\nwithout Reasonable\nAccess to any\nMode of Transit",
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_transitU18access.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_transitU18access.png")
```

```{r tabTransitAccessU18, fig.align = "center", fig.cap="Persons under 18 in Census Block Groups without reasonable access to any mode of transit."}
vt_blkgrps_sf_noTransit %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Under 18` = sum(NewUnder18,na.rm = TRUE),
            `Pct of Under 18` = round(sum(NewUnder18)/max(TotalUnder18)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Pct of Under 18` > 0) %>% 
  mutate(`Pct of Under 18` = if_else(`Pct of Under 18` > 100,100,`Pct of Under 18`)) %>% 
  mutate(`Pct of Under 18` = paste0(`Pct of Under 18`,"%")) %>% 
  arrange(desc(`Under 18`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Persons Under 18 without Access to any Mode of Transit", align = "r", col.names = c(names(.)[1:2],"Number of Under 18 in Census Block Groups","Pct of Under 18 in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped")) 
```



## Frequency of Transit Service
The quality of access to public transit can also be assessed by the level or frequency of service, which is measured by "headway." Headway describes the scheduled or observed time between transit vehicle arrivals at a transit stop. For example, a bus route with a 15 minute headway would mean that the bus arrives at a given stop four times an hour (it's frequency), or once every 15 minutes. Transit headways are affected by the scheduled frequency of service, the number of vehicles on a given route, traffic delays, and dispatch management of vehicle spacing.[^mindGap] Headways are significant because they affect the desirability and useability of transit service.[^headways] Headways can affect:

* average wait times
* the amount of planning and preparation needed to use transit and stay on schedule
* the amount of time lost when transit schedules do not directly conform to work, school, or activity schedules
* the time penalty for missing a train or bus
* public use or support of transit

Headways are significant for transit dependent populations and can affect the quality of life and economic opportunities of transit riders. Figure \@ref(fig:mapGTFSbus) below shows average scheduled headways for bus routes throughout Vermont.[^headwayReal] Headways vary significantly by route and by type of service across the state. 

```{r GTFS, include=FALSE}
# Calculate average headway for block groups within walking distance of bus stops and then compute weighted average headway by population group
vt_bus_stopHeadway_df <- vt_stops_sf_Bus %>% 
  st_transform(., crs = 2852) %>% 
  st_join(., vt_busBuff400m_sf) %>% 
  as.data.frame() %>% 
  group_by(GEOID) %>% 
  summarize(AvgStopHeadway = mean(headway,na.rm = TRUE))
# Create same for tracts
vt_bus_stopHeadwayTracts_df <- vt_stops_sf_Bus %>% 
  st_transform(., crs = 2852) %>% 
  st_join(., vt_busBuff400mTracts_sf) %>% 
  as.data.frame() %>% 
  group_by(GEOID) %>% 
  summarize(AvgStopHeadway = mean(headway,na.rm = TRUE))
# Create version of vt_bus_stopHeadwayTracts_df for rbind
vt_bus_stopHeadwayTracts <- vt_bus_stopHeadwayTracts_df %>% 
  as.data.frame() %>% 
  left_join(vt_busBuff400mTracts_sf,.,by="GEOID") %>% 
  as.data.frame() %>% 
  transmute(Disabled = NewDisabled,
            `No Car HH` = NewNoCar,
            AvgStopHeadway) %>% 
  gather(key = Group, value = Pop, Disabled:`No Car HH`) %>% 
  group_by(Group) %>% 
  summarize(AvgWHeadway = weighted.mean(x = AvgStopHeadway, 
                                        w = Pop, na.rm = TRUE))

# Join stop average headway to block groups within buffer
vt_busHeadwayPops_df <- vt_busBuff400m_sf %>% 
  as.data.frame() %>% 
  left_join(.,vt_bus_stopHeadway_df, by = "GEOID") %>% 
  transmute(`Total Pop` = NewPop,
            Minority = NewMinority,
            `Under 5` = NewUnder5,
            `Over 64` = NewOver64,
            `Under 18` = NewUnder18,
            `English Limited HH` = NewEng_limit,
            `Low Income` = NewPov,
            `No HS Dip` = NewLths,
            AvgStopHeadway) %>% 
  gather(key = Group, value = Pop, `Total Pop`:`No HS Dip`) %>% 
  group_by(Group) %>% 
  summarize(AvgWHeadway = weighted.mean(x = AvgStopHeadway, 
                                        w = Pop, na.rm = TRUE)) %>% 
  rbind(.,vt_bus_stopHeadwayTracts)
```

```{r mapGTFSbus, fig.align = "center", fig.cap="Mean headways by bus route for Vermont routes running Monday through Friday, from 6am - 9pm."}
# arrange row order descending so that higher frequency lines are on top
vt_routes_sf_Bus <- vt_routes_sf_Bus %>% 
  arrange(desc(mean_headways)) %>% 
  # st_simplify(., dTolerance = 100) %>% 
  st_make_valid()

# m <- tm_shape(vt_state_sf_cb, unit = "mi") +
#   tm_fill(col = "gray",alpha = 0.2) +
#   tm_shape(vt_towns_sf, unit = "mi") +
#   tm_borders(col = "gray", lwd = 0.2) +
#   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(vt_routes_sf_Bus) +
#   tm_lines(col = "mean_headways", style = "quantile", lwd = 2,
#            title.col = "Minutes", 
#            legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = c("#1a9641","#a6d96a","#ffff00","#fdae61","#d7191c")) +
#   tm_shape(I89roadSegment) +
#   tm_symbols(shape = I89, border.lwd = NA, size = 0.1) +
#   tm_shape(I91roadSegment) +
#   tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
#   tm_shape(I91roadSegment2) +
#   tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
#   tm_shape(I93roadSegment) +
#   tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
#   tm_shape(vt_towns_sf_pts) + tm_dots() +
#   tm_text("NAME", size = 0.5, col = "black",
#           xmod = 0.7, ymod = 0.2, shadow = TRUE) +
#   tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5,
#                position = c(0.6,0.005)) +
#   tm_layout(title = "Mean\nHeadways\nby Bus\nRoute",
#             frame = FALSE, main.title.size = 0.8,
#             legend.outside = TRUE,
#             legend.title.size = 0.8,
#             legend.outside.position = c("right", "top"))
# 
# tmap_save(m, "DATA/transport/VT/figures/vt_busHeadways.png",
#           height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_busHeadways.png")
```

Table \@ref(tab:tabGTFSbusTransit) shows headways for the bus transit districts for which GTFS data was available. The average population-weighted headway for all bus routes in the state is `r round(vt_busHeadwayPops_df %>% filter(Group == "Total Pop") %>% select(AvgWHeadway) %>% pull(),1)` minutes.  

```{r tabGTFSbusTransit, fig.align = "center", fig.cap="Headway by rapid transit line based on on GTFS schedule for Monday to Friday service, 6am - 9pm."}
# Show table of mean headways by line
vt_stops_sf_Bus %>% 
  as.data.frame() %>% 
  # mutate(ROUTE = as.character(ROUTE)) %>% 
  # left_join(., as.data.frame(vt_routes_sf), by=c("ROUTE" = "route_id")) %>% 
  group_by(agency) %>% 
  summarize(Number = n(), Min = min(headway), Mean = mean(headway), Max = max(headway)) %>% 
  arrange(agency) %>% 
  kableExtra::kable(longtable = T, booktabs = T, digits = 1,
                    format.args = list(big.mark=","),
                    col.names = c("District","Number of Stops",
                                  "Min",
                                  "Mean",
                                  "Max"), 
                    caption = "Bus Stop Headways by Bus District", 
                    align = "r") %>% 
  kableExtra::add_header_above(c(" " = 2, "Headway (minutes)" = 3)) %>%
  # kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header"))
```

However headways, or service frequencies, vary significantly across population groups. Figure \@ref(fig:loliBusHeadway) below compares the population-weighted average bus headways for different groups across the state. Most groups have average headways shorter than the general population. However, persons under 18, disabled persons, and adults without a high school diploma have headways that exceed that of the general population. The shortest population-weighted average headway, experienced by limited English speaking households is `r round(vt_busHeadwayPops_df %>% arrange(AvgWHeadway) %>% slice(.,1) %>% select(AvgWHeadway) %>% pull(),0)` minutes.

```{r loliBusHeadway, echo=FALSE, message=FALSE, warning=FALSE, fig.align = "center", fig.cap="Population-weighted average bus headways for groups within 400 meters of a bus stop across Vermont."}
vt_busHeadwayPops_df %>% 
  mutate(Group = recode(Group, "Minority" = "People of Color",
                        "No HS Dip" = "No HS Diploma")) %>%
  ggplot(aes(x = reorder(Group,AvgWHeadway), y = AvgWHeadway)) +
  geom_segment(aes(x = reorder(Group,AvgWHeadway), 
                   xend = reorder(Group,AvgWHeadway),
                   y = pull(vt_busHeadwayPops_df[6,2]), yend = AvgWHeadway), 
               color = "tan1") +
  geom_point(color = "orange", size = 4, alpha = 0.8) +
  coord_flip() + xlab("") + ylab("minutes") + ggtitle("Population-Weighted Average Stop Headway for\nVermont Groups within Walking Distance (400m)\nof Bus Stops") + theme_light() +
  theme(panel.grid.major.y = element_blank(),
        panel.border = element_blank(),
        axis.ticks.y = element_blank()) +
  geom_text(aes(x = Group, y = AvgWHeadway + 0.2 * sign(AvgWHeadway), 
                label = round(AvgWHeadway,0)), 
            hjust = 1.5, vjust = -0.6, size = 3,
            color=rgb(100,100,100, maxColorValue=255)) +
  geom_hline(yintercept = pull(vt_busHeadwayPops_df[6,2]), 
             linetype = "dashed") +
  geom_text(aes(x = "Total Pop", y = 305, label = "Above state\naverage"),
            color = "gray48") +
  geom_text(aes(x = "Total Pop", y = 277, label = "Below state\naverage"),
            color = "gray48") +
  expand_limits(y = c(250,320))

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

```{r excessiveHeadway, include=FALSE}
# Subset block groups and tracts with average headways exceeding 80th percentile headway for all routes
headway80thblkgrps_bus <- vt_busBuff400m_sf %>% 
  left_join(., vt_bus_stopHeadway_df, by = "GEOID") %>% 
  filter(AvgStopHeadway >= quantile(vt_routes_sf_Bus$mean_headways,0.8))

headway80thtracts_bus <- vt_busBuff400mTracts_sf %>% 
  left_join(., vt_bus_stopHeadwayTracts_df, by = "GEOID") %>%
  filter(AvgStopHeadway >= quantile(vt_routes_sf_Bus$mean_headways,0.8))

# calculate numbers of affected people by mode by municipality
headway80_1 <- headway80thblkgrps_bus %>%
  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(NewUnder5,na.rm = TRUE),
            `Pct of Under 5` = round(sum(NewUnder5)/max(TotalUnder5)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  mutate(`Pct of Under 5` = paste0(`Pct of Under 5`,"%")) %>% 
  arrange(desc(`Under 5`))

headway80_2 <- headway80thblkgrps_bus %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Under 18` = sum(NewUnder18,na.rm = TRUE),
            `Pct of Under 18` = round(sum(NewUnder18)/max(TotalUnder18)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  mutate(`Pct of Under 18` = paste0(`Pct of Under 18`,"%")) %>% 
  arrange(desc(`Under 18`))

# write out to file for use in later analysis
save(headway80thblkgrps_bus,
     headway80thtracts_bus,
     vt_bus_stopHeadway_df, vt_busHeadwayPops_df,
     file = "DATA/transport/VT/headway80th.Rds")

st_write(headway80thblkgrps_bus,"DATA/transport/VT/headway80thblkgrps_bus.shp",
         delete_layer = TRUE)
st_write(headway80thtracts_bus, "DATA/transport/VT/headway80thtracts_bus.shp",
         delete_layer = TRUE)
```
The population-weighted average headways above provide a glimpse of average headways across the state for different population groups. However, headways also vary for specific communities. Indeed, it is important to consider priority populations or transportation-limited populations that experience excessively long headways or low service frequency. Figures \@ref(fig:mapHeadwayUnder5bus) and \@ref(fig:mapHeadwayUnder18bus) below highlight communities of persons under 5 and under 18 who are within reasonable distance of bus stops (i.e., 400m or 1/4 mile), but where the mean headway for these stops exceeds the 80th percentile of mean route headways for the state. Tables \@ref(tab:tabHeadwayUnder5bus) and \@ref(tab:tabHeadwayUnder18bus) show breakdowns by municipality. For example, Table \@ref(tab:tabHeadwayUnder5bus) shows that `r formatC(as.numeric(headway80_1[1,3]),big.mark=",")` children under 5 across `r headway80_1[1,2]` Block Groups in `r headway80_1[1,1]` resided in areas where the average headway of accessible bus stops exceeded `r quantile(vt_routes_sf_Bus$mean_headways,0.8)` minutes. These people represented `r headway80_1[1,4]` of children under 5 in `r headway80_1[1,1]`. Similarly, Table \@ref(tab:tabHeadwayUnder18bus) shows that `r formatC(as.numeric(headway80_2[1,3]),big.mark=",")` persons under 18 across `r headway80_2[1,2]` Block Groups in `r headway80_2[1,1]` resided in areas where the average headway of accessible bus stops exceeded `r quantile(vt_stops_sf_Bus$headway,.8)` minutes. These people represented `r headway80_2[1,4]` of persons under 18 in `r headway80_2[1,1]`.
Maps and tables of other priority populations with excessive headways can be found in Appendix B. 

```{r mapHeadwayUnder5bus, fig.align = "center", fig.cap="Children under 5 within 400m or 1/4 mile of transit stops in which the mean headway for these stops is above the 80th percentile for the state."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") +
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(vt_towns_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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(headway80thblkgrps_bus) +
  tm_fill(col = "NewUnder5", style = "quantile",
          colorNA = "white", colorNULL = "gray",
          title = "Number of\nChildren Under 5",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd") +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5,
               position = c(0.6,0.005)) +
  # tm_add_legend(type = "line", lwd=1, alpha = 0.3, labels = "Bus Route", col = "blue") +
  tm_layout(title = paste("Children Under 5","with Bus Headways",paste0("Over ",round(quantile(vt_routes_sf_Bus$mean_headways,0.8),0)," Minutes"),sep = "\n"),
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_headwayUnder5bus.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_headwayUnder5bus.png")
```

```{r tabHeadwayUnder5bus, fig.align = "center", fig.cap="Children under 5 within 400m or 1/4 mile of transit stops in which the mean headway for these stops is above the 80th percentile for the state."}
headway80thblkgrps_bus %>%
  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(NewUnder5,na.rm = TRUE),
            `Pct of Under 5` = round(sum(NewUnder5)/max(TotalUnder5)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Pct of Under 5` > 0) %>% 
  mutate(`Pct of Under 5` = if_else(`Pct of Under 5` > 100, 100, `Pct of Under 5`)) %>% 
  mutate(`Pct of Under 5` = paste0(`Pct of Under 5`,"%")) %>% 
  arrange(desc(`Under 5`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = paste("Children Under 5","with Average Headways",paste0("Over ",round(quantile(vt_routes_sf_Bus$mean_headways,0.8),0)," Minutes"),sep = "\n"), align = "r", col.names = c(names(.)[1:2],"Number of Under 5 in Block Groups","Pct of Under 5 in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))

```

```{r mapHeadwayUnder18bus, fig.align = "center", fig.cap="Persons under 18 within 400m or 1/4 mile of transit stops in which the mean headway for these stops is above the 80th percentile for the state."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") +
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(vt_towns_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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(headway80thblkgrps_bus) +
  tm_fill(col = "NewUnder18", style = "quantile",
          colorNA = "white", colorNULL = "gray",
          title = "Number of\nPersons Under 18",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd") +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5,
               position = c(0.6,0.005)) +
  # tm_add_legend(type = "line", lwd=1, alpha = 0.3, labels = "Bus Route", col = "blue") +
  tm_layout(title = paste("Persons Under 18","with Bus Headways",paste0("Over ",round(quantile(vt_routes_sf_Bus$mean_headways,0.8),0)," Minutes"),sep = "\n"),
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_headwayUnder18bus.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_headwayUnder18bus.png")
```

```{r tabHeadwayUnder18bus, fig.align = "center", fig.cap="Persons under 18 within 400m or 1/4 mile of transit stops in which the mean headway for these stops is above the 80th percentile for the state."}
headway80thblkgrps_bus %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Under 18` = sum(NewUnder18,na.rm = TRUE),
            `Pct of Under 18` = round(sum(NewUnder18)/max(TotalUnder18)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Pct of Under 18` > 0) %>% 
  mutate(`Pct of Under 18` = if_else(`Pct of Under 18` > 100, 100, `Pct of Under 18`)) %>% 
  mutate(`Pct of Under 18` = paste0(`Pct of Under 18`,"%")) %>% 
  arrange(desc(`Under 18`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = paste("Persons Under 18","with Average Headways",paste0("Over ",round(quantile(vt_routes_sf_Bus$mean_headways,0.8),0)," Minutes"),sep = "\n"), align = "r", col.names = c(names(.)[1:2],"Number of Under 18 in Block Groups","Pct of Under 18 in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))

```



## Walkability

More walkable communities are positively associated with a variety of quality of life and public health benefits. Conversely, less walkable communities are associated with a range of public health challenges, including higher rates of cardiovascular disease and diabetes.[^WalkCardio] The walkability of neighborhoods is increasingly regarded as an important component of healthier and more sustainable communities, particularly with regard to reducing motor vehicle travel and promoting alternative modes of transport that are more affordable, accessible, and less polluting.  

The National Walkability Index is a nationwide geographic data resource produced by the U.S. Environmental Protection Agency that ranks Census Block Groups according to their relative walkability.[^Walkability] Walkability was modeled based on characteristics of the built environment that influence the likelihood of walking being used as a mode of travel. These characteristics of Block Groups include:

* mix of employment types (greater mix = more walkable)
* amount or density of occupied housing (higher density = more walkable)
* street intersection density (higher density = more walkable)
* proportion of workers who carpool (more carpooling = more walkable)

```{r walkabilityData, include=FALSE}
# Analyze Walkability Index for VT. See https://edg.epa.gov/metadata/catalog/search/resource/details.page?uuid=%7B251AFDD9-23A7-4068-9B27-A3048A7E6012%7D
# Import walkability layer
walkability <- st_read(dsn = "DATA/Walkability/Natl_WI_SHP",
        layer = "Natl_WI") 

# Join NatWalkIndex values to other block group variables
vt_blkgrp_walkability_sf <- walkability %>% 
  as.data.frame() %>% 
  select(GEOID10,SFIPS,CFIPS,TRFIPS,NatWalkInd) %>% 
  left_join(vt_blkgrps_sf,., by = c("GEOID" = "GEOID10"))

# Aggregate NatWalkIndex values to tract level sf
vt_tract_walkability_sf <- vt_blkgrp_walkability_sf %>% 
  as.data.frame() %>% 
  select(SFIPS,CFIPS,TRFIPS,NatWalkInd) %>% 
  mutate(TRAVTID = paste0(SFIPS,CFIPS,TRFIPS)) %>% 
  group_by(TRAVTID) %>% 
  summarize(NatWalkInd = mean(NatWalkInd)) %>% 
  left_join(vt_tracts_sf,., by = c("GEOID" = "TRAVTID"))

# Create summarized df of population weighted avg walkability by tract pop
vt_tract_walkability_df <- vt_tract_walkability_sf %>% 
  as.data.frame() %>% 
  transmute(Disabled = disabledOver18E,
            `No Car HH` = HHnoCarE,
            NatWalkInd) %>% 
  gather(key = Group, value = Pop, Disabled:`No Car HH`) %>% 
  group_by(Group) %>% 
  summarize(wNatWalkIndex = weighted.mean(NatWalkInd, Pop, na.rm = TRUE))

# Populated-Weighted Average Walkability Index for Priority Populations
vt_wAvgWalkability_df <- vt_blkgrp_walkability_sf %>% 
  as.data.frame() %>% 
  transmute(`Total Pop` = totalpopE,
            Minority = minorityE,
            `Low Income` = num2povE,
            `Under 5` = under5E,
            `Under 18` = under18E,
            `Over 64` = over64E,
            `Limited English HH` = eng_limitE,
            `No HS Dip` = lthsE,
            NatWalkInd) %>% 
  gather(key = Group, value = Pop, `Total Pop`:`No HS Dip`) %>% 
  group_by(Group) %>% 
  summarize(wNatWalkIndex = weighted.mean(NatWalkInd, Pop, na.rm = TRUE)) %>% 
  rbind(., vt_tract_walkability_df)

# save data for later analysis
save(vt_blkgrp_walkability_sf,
     vt_tract_walkability_sf,
     vt_wAvgWalkability_df,
     file = "DATA/transport/VT/walkability.Rds")
```

Figure \@ref(fig:mapWalkability) shows walkability scores by Census Block Group across Vermont.  The population-weighted average walkability score across the state is `r round(vt_wAvgWalkability_df %>% filter(Group == "Total Pop") %>% select(wNatWalkIndex) %>% pull(),1)`. Walkability varies across the state, although most of the state is classified as below average to least walkable. Walkability is highest in portions of Burlington, Montpelier, Fairhaven, and Rutland.

```{r mapWalkability, echo=FALSE, message=FALSE, warning=FALSE, fig.align = "center", fig.cap="Walkability Index scores for Census Block Groups across Vermont."}
m <- tm_shape(vt_blkgrp_walkability_sf) +
  tm_fill(col = "NatWalkInd", alpha = 0.8, style = "equal", n = 4,
          colorNA = "white", colorNULL = "gray", title = "Walkability Index",
          legend.format = list(fun = function(x) formatC(x, digits = 1, format = "f")), palette = c("#FECC61","#FEF6BA","#B9E293","#238443"), 
          labels = c("1 - 5.8 (Least Walkable)",
                     "5.8 - 10.5 (Below Avg)",
                     "10.5 - 15.2 (Above Avg)",
                     "15.2 - 20 (Most Walkable)")) +
  tm_shape(vt_towns_sf, unit = "mi") + 
  tm_borders(col = "gray", lwd = 0.2) +
  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.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "Walkability", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_walk.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_walk.png")
```

Walkability varies for different population groups across the state. Figure \@ref(fig:loliWalkability) below compares population-weighted average walkability index scores for different populations. Most groups live in communities with below average walkability index scores that are also at or above that of the general population. However, persons over age 64 tend to live in communities with walkability scores below the state average. 

```{r loliWalkability, echo=FALSE, message=FALSE, warning=FALSE, fig.align = "center", fig.cap="population-weighted average walkability index scores for groups across Vermont."}
# Create lollipop graph of pop-weighted average walk index
vt_wAvgWalkability_df %>% 
  mutate(Group = recode(Group, "Minority" = "People of Color",
                        "No HS Dip" = "No HS Diploma")) %>%
  ggplot(aes(x = reorder(Group,-wNatWalkIndex), y = wNatWalkIndex)) +
  geom_segment(aes(x = reorder(Group,-wNatWalkIndex), 
                   xend = reorder(Group,-wNatWalkIndex),
                   y = pull(vt_wAvgWalkability_df[6,2]), yend = wNatWalkIndex), 
               color = "green4") +
  geom_point(color = "green", size = 4, alpha = 0.8) +
  coord_flip() + xlab("") + ylab("Walk Index") + ggtitle("Population-Weighted Walkability Index") + theme_light() +
  theme(panel.grid.major.y = element_blank(),
        panel.border = element_blank(),
        axis.ticks.y = element_blank()) +
  geom_text(aes(x = Group, y = wNatWalkIndex + 0.2 * sign(wNatWalkIndex), 
                label = round(wNatWalkIndex,1)), 
            hjust = 2.5, vjust=-0.8, size = 3,
            color=rgb(100,100,100, maxColorValue=255)) +
  geom_hline(yintercept = pull(vt_wAvgWalkability_df[6,2]), linetype = "dashed") +
  geom_text(aes(x = "Total Pop", y = 7.4, label = "Above state\naverage"),
            color = "gray48") +
  geom_text(aes(x = "Total Pop", y = 6.3, label = "Below state\naverage"),
            color = "gray48") +
  expand_limits(y = c(6,8))

ggsave("images/VT_Walk_graph.png")
```
```{r lowWalk, include=FALSE}
# Isolate least walkable block groups and tracts and add area variable
leastWalkable <- vt_blkgrp_walkability_sf %>% 
  filter(NatWalkInd <= 5.8) %>% 
  mutate(SqKm = bg_area_m2/10^6)

leastWalkableTract <- vt_tract_walkability_sf %>% 
  filter(NatWalkInd <= 5.8) %>% 
  mutate(SqKm = as.numeric(st_area(.))/10^6)

leastWalkable1 <- leastWalkable %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Under 18` = sum(under18E,na.rm = TRUE),
            `Pct of Under 18` = round(sum(under18E)/max(TotalUnder18)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  mutate(`Pct of Under 18` = paste0(`Pct of Under 18`,"%")) %>% 
  arrange(desc(`Under 18`))

leastWalkable2 <- leastWalkable %>%
  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)) %>% 
  rename("City/Town" = NAME) %>% 
  mutate(`Pct of Over 64` = paste0(`Pct of Over 64`,"%")) %>% 
  arrange(desc(`Over 64`))

```
The population-weighted average walkability scores provide a glimpse of average walkability across the state for different population groups. However, walkability also varies for specific communities. Indeed, it is important to consider priority populations that reside in communities with low walkability. Figures \@ref(fig:mapLowWalkU18) and  \@ref(fig:mapLowWalkO64) below highlight Census Block Groups across Vermont that are least walkable and the numbers of persons under 18 or over 64 that live there. Tables \@ref(tab:tabLowWalkU18) and \@ref(tab:tabLowWalkO64) show breakdowns by municipality. For example, Table \@ref(tab:tabLowWalkU18) shows that `r formatC(as.numeric(leastWalkable1[1,3]),big.mark=",")` people under 18 across `r leastWalkable1[1,2]` Block Groups in `r leastWalkable1[1,1]` resided in the least walkable areas. These people represented `r leastWalkable1[1,4]` of people under 18 in `r leastWalkable1[1,1]`. Similarly, Table \@ref(tab:tabLowWalkO64) shows that `r formatC(as.numeric(leastWalkable2[1,3]),big.mark=",")` people over 64 across `r TransitNoAccess1[1,2]` Block Groups in `r leastWalkable2[1,1]` resided the least walkable areas. These people represented `r leastWalkable2[1,4]` of people over 64 in `r leastWalkable2[1,1]`.
Maps and tables of other priority populations in the least walkable Block Groups can be found in Appendix B.

```{r mapLowWalkU18, fig.align = "center", fig.cap="Persons under 18 per square kilometer in least walkable Census Block Groups (i.e., Walkability Index Score less than or equal to 5.8)."}
leastWalkable18 <- leastWalkable %>% 
  select(under18E,SqKm) %>% 
  filter(under18E > 0)

m <- tm_shape(vt_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(leastWalkable, unit = "mi") +
  tm_fill(col = "under18E", alpha = 0.9, style = "fisher",
          colorNA = "white", colorNULL = "gray", title = "Persons Under 18\nPer Square Kilometer",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd", convert2density = TRUE, area = "SqKm") +
  tm_shape(vt_tracts_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "Persons Under 18\nin Least Walkable\nBlock Groups", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_walkU18.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_walkU18.png")
```

```{r tabLowWalkU18, fig.align = "center", fig.cap="Persons under 18 by municipality in least walkable Census Block Groups (i.e., Walkability Index Score less than or equal to 5.8)."}
leastWalkable %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Under 18` = sum(under18E,na.rm = TRUE),
            `Pct of Under 18` = round(sum(under18E)/max(TotalUnder18)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Under 18` > 0 & `Pct of Under 18` > 0) %>% 
  mutate(`Pct of Under 18` = if_else(`Pct of Under 18` > 100, 100, `Pct of Under 18`)) %>% 
  mutate(`Pct of Under 18` = paste0(`Pct of Under 18`,"%")) %>% 
  arrange(desc(`Under 18`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Persons Under 18 in Least Walkable Block Groups", align = "r", col.names = c(names(.)[1:2],"Number of Under 18 in Block Groups","Pct of Under 18 in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapLowWalkO64, fig.align = "center", fig.cap="Persons over 64 per square kilometer in least walkable Census Block Groups (i.e., Walkability Index Score less than or equal to 5.8)."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(leastWalkable, unit = "mi") +
  tm_fill(col = "over64E", alpha = 0.9, style = "fisher",
          colorNA = "white", colorNULL = "gray", title = "Persons Over 64\nPer Square Kilometer",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd", convert2density = TRUE, area = "SqKm") +
  tm_shape(vt_tracts_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "Persons Over 64\nin Least Walkable\nBlock Groups", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_walkO64.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_walkO64.png")
```

```{r tabLowWalkO64, fig.align = "center", fig.cap="Persons over 64 by municipality in least walkable Census Block Groups (i.e., Walkability Index Score less than or equal to 5.8)."}
leastWalkable %>%
  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)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Over 64` > 0 & `Pct of Over 64` > 0) %>% 
  mutate(`Pct of Over 64` = if_else(`Pct of Over 64` > 100, 100, `Pct of Over 64`)) %>% 
  mutate(`Pct of Over 64` = paste0(`Pct of Over 64`,"%")) %>% 
  arrange(desc(`Over 64`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Persons Over 64 in Least Walkable Block Groups", align = "r", col.names = c(names(.)[1:2],"Number of Over 64 in Block Groups","Pct of Over 64 in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```



## Transportation Cost Burden

Transportation is the second-largest expenditure category for American households, approximately 16% of annual expenditures on average between 2015 and 2018.[^Expenditures] Only housing costs (~32%) exceed those of transportation. The high cost of transportation is due in large part to the costs of dependence on private vehicle ownership, including purchase and financing, insurance, maintenance, and fuel costs. Households without reasonable access to other modes of transportation typically experience higher transportation costs.[^CostBurden] Transportation costs can be significant economic burdens for moderate and lower income households and for communities where destinations are distant or limited. For many working-class and rural households, transportation costs exceed housing costs. 

The Location Affordability Index (LAI) is a nationwide geographic data resource developed by the U.S. Department of Housing and Urban Development (HUD) in collaboration with the U.S. Department of Transportation under the federal Partnership for Sustainable Communities.[^LAIabout] The LAI was developed as a way of integrating both housing and transportation costs into estimates of the affordability of specific neighborhoods and cities. The LAI provides estimated total expenses from transportation, whether transit or motor vehicle, as well as transportation cost as a percentage of household income, for eight different household types — which vary by household income, size, and number of commuters — at the Census Tract level. For the purposes of this analysis, transportation costs as a percentage of income for moderate-income households (households with one commuter making 80% or less of area median household income) was used as a proxy to represent transportation cost burden for all areas. LAI data were also downscaled to the Block Group level for the purpose of demographic analyses and comparisons (see Appendix B Data and Methodology). 

```{r TcostBurden, include=FALSE}
# Read in Location Affordability Index data. Using Version 3 based on ACS 2012-2016 data at Tract level. Downloaded from https://hudgis-hud.opendata.arcgis.com/datasets/location-affordability-index-v-3?selectedAttribute=avg_hh_size_renters
lai_ri <- st_read("DATA/transport/Location_Affordability_Index_v.3", "Location_Affordability_Index_v.3") %>% 
  filter(STUSAB == "VT") %>% 
  st_make_valid() %>% 
  st_transform(., crs = 2852)

# Downscale LAI to Block Group by assigning tract value to all overlapping block groups (i.e. block groups with same tract GEOID prefix)
lai_vt_blkgrp <- vt_blkgrps_sf %>% 
  mutate(GEOIDT = str_sub(GEOID,1,11)) %>% 
  left_join(., as.data.frame(lai_ri),by=c("GEOIDT"="GEOID"))

# Join LAI to tracts
lai_vt_tract <- vt_tracts_sf %>% 
  left_join(., as.data.frame(lai_ri),by="GEOID")

# create summarized df of population weighted average transportation cost burden by tract pop
lai_vt_tract_df <- lai_vt_tract %>% 
  as.data.frame() %>% 
  transmute(Disabled = disabledOver18E,
            `No Car HH` = HHnoCarE,
            hh7_t) %>% 
  gather(key = Group, value = Pop, Disabled:`No Car HH`) %>% 
  group_by(Group) %>% 
  summarize(wTBurden = weighted.mean(hh7_t, Pop, na.rm = TRUE))

# Population-weighted average transportation cost burden for priority populations
lai_vt_df <- lai_vt_blkgrp %>% 
  as.data.frame() %>% 
  transmute(`Total Pop` = totalpopE,
            Minority = minorityE,
            `Low Income` = num2povE,
            `Under 5` = under5E,
            `Under 18` = under18E,
            `Over 64` = over64E,
            `Limited English HH` = eng_limitE,
            `No HS Dip` = lthsE,
            hh7_t) %>% 
  gather(key = Group, value = Pop, `Total Pop`:`No HS Dip`) %>% 
  group_by(Group) %>% 
  summarize(wTBurden = weighted.mean(hh7_t, Pop, na.rm = TRUE)) %>% 
  rbind(., lai_vt_tract_df)

# save data for later analysis
save(lai_vt_blkgrp,
     lai_vt_tract,
     lai_vt_df,
     file = "DATA/transport/VT/costBurden.Rds")

# isolate places where transportation cost as a percentage of income exceeds housing as a percentage of income
lai_ToverH <- lai_ri %>%
  filter(hh7_t > hh7_h)
```

Figure \@ref(fig:mapLAI) below shows transportation cost burden (i.e., transportation cost as a percentage of household income) by Census Tract across Vermont. Transportation cost burden varies across the state. It is lowest around Burlington, Middlebury and St. Albans. In `r nrow(lai_ToverH)` Census Tracts, the transportation cost burden exceeds the cost burden of housing. 

```{r mapLAI, fig.align = "center", fig.cap="Transportation cost burden for Census Tracts across Vermont."}
m <- tm_shape(lai_ri) +
  tm_fill(col = "hh7_t", alpha = 0.8, style = "quantile", title = "Transportation as\nPercent of Income", showNA = FALSE,
          legend.format = list(fun = function(x) formatC(x, digits = 1, format = "f"))) +
  tm_shape(lai_ToverH) + tm_dots(col = "gray", shape = 8, size = 0.5) +
  tm_shape(vt_towns_sf, unit = "mi") + 
  tm_borders(col = "gray", lwd = 0.2) +
  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.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_add_legend(type = "symbol",labels = "Transport > Housing burden", shape = 8, size = 0.5, col = "gray") +
  tm_layout(title = "Transportation\nCost Burden", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_tburden.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_tburden.png")
````

Transportation cost burden varies for different population groups across the state. Figure \@ref(fig:loliLAI) below compares population-weighted average transportation cost burden for different populations. Most groups live in communities with transportation cost burdens as a percentage of household income that are greater than the general population. Persons over age 64 tend to live in communities with the highest transportation cost burdens as a percentage of household income. 

```{r loliLAI, fig.align = "center", fig.cap="Population-weighted average transportation cost burden as a percentagae of household income for groups across Vermont."}
# Create lollipop graph of pop-weighted average LAI
lai_vt_df %>% 
  mutate(Group = recode(Group, "Minority" = "People of Color",
                        "No HS Dip" = "No HS Diploma")) %>%
  ggplot(aes(x = reorder(Group,wTBurden), y = wTBurden)) +
  geom_segment(aes(x = reorder(Group,wTBurden), 
                   xend = reorder(Group,wTBurden),
                   y = pull(lai_vt_df[6,2]), yend = wTBurden), 
               color = "coral4") +
  geom_point(color = "coral3", size = 4, alpha = 0.8) +
  coord_flip() + xlab("") + ylab("Transportation as Percent of Income") + ggtitle("Population-Weighted Transportation Cost Burden") + theme_light() +
  theme(panel.grid.major.y = element_blank(),
        panel.border = element_blank(),
        axis.ticks.y = element_blank()) +
  geom_text(aes(x = Group, y = wTBurden + 0.2 * sign(wTBurden), 
                label = round(wTBurden,1)), 
            hjust = 0.6, vjust=-0.8, size = 3,
            color=rgb(100,100,100, maxColorValue=255)) +
  scale_y_continuous(labels = function(x) paste0(x, "%")) +
  geom_hline(yintercept = pull(lai_vt_df[6,2]), linetype = "dashed") +
  geom_text(aes(x = "Total Pop", y = 29, label = "Above state\naverage"),
            color = "gray48") +
  geom_text(aes(x = "Total Pop", y = 27, label = "Below state\naverage"),
            color = "gray48") +
  expand_limits(y = c(25,30))

ggsave("images/VT_Cost_graph.png")
```
```{r HiBurden, include=FALSE}
# Isolate highest transport burden block groups and tracts and add area variable
HiBurden <- lai_vt_blkgrp %>% 
  filter(percent_rank(hh7_t) > 0.8) %>% 
  mutate(SqKm = bg_area_m2/10^6)

HiBurdenTract <- lai_vt_tract %>% 
  filter(percent_rank(hh7_t) > 0.8) %>% 
  mutate(SqKm = as.numeric(st_area(.))/10^6)

HiBurden1 <- HiBurden %>%
  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)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Pct of Over 64` > 0) %>% 
  mutate(`Pct of Over 64` = if_else(`Pct of Over 64` > 100, 100, `Pct of Over 64`)) %>% 
  mutate(`Pct of Over 64` = paste0(`Pct of Over 64`,"%")) %>% 
  arrange(desc(`Over 64`))

HiBurden2 <- HiBurden %>%
  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)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Pct of No HS Dip` > 0) %>% 
  mutate(`Pct of No HS Dip` = if_else(`Pct of No HS Dip` > 100, 100, `Pct of No HS Dip`)) %>%
  mutate(`Pct of No HS Dip` = paste0(`Pct of No HS Dip`,"%")) %>% 
  arrange(desc(`No HS Dip`))
```
The population-weighted average transportation cost burden values provide a glimpse of average transportation cost burden across the state for different population groups. However, transportation cost burden also varies for specific communities. Indeed, it is important to consider priority populations that reside in communities with a high transportation cost burden. Figures \@ref(fig:mapHiBurdenO64) and \@ref(fig:mapHiBurdenLTHS) below highlight Census Block Groups across Vermont where the transportation cost burden is in the highest quintile (i.e., 80th percentile) and the numbers of persons over 64 or adults without a high school diploma that live there. Tables \@ref(tab:tabHiBurdenO64) and \@ref(tab:tabHiBurdenLTHS) show breakdowns by municipality. For example, Table \@ref(tab:tabHiBurdenO64) shows that `r formatC(as.numeric(HiBurden1[1,3]),big.mark=",")` people over 64 across `r HiBurden1[1,2]` Block Groups in `r HiBurden1[1,1]` resided in areas with the highest transportation cost burden. These people represented `r HiBurden1[1,4]` of people over 64 in `r HiBurden1[1,1]`. Similarly, Table \@ref(tab:tabHiBurdenLTHS) shows that `r formatC(as.numeric(HiBurden2[1,3]),big.mark=",")` adults without a high school diploma across `r HiBurden2[1,2]` Block Groups in `r HiBurden2[1,1]` resided in areas with the highest transportation cost burden. These people represented `r HiBurden2[1,4]` of adults without a high school diploma in `r HiBurden2[1,1]`.
Maps and tables of other priority populations in areas with the highest transportation cost burden can be found in Appendix B.

```{r mapHiBurdenO64, fig.align = "center", fig.cap="Persons over 64 per square kilometer in Census Block Groups with the highest transportation cost burden (i.e., transportation cost burden in the 80th percentile for the state)."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(HiBurden, unit = "mi") +
  tm_fill(col = "over64E", alpha = 0.9, style = "fisher",
          colorNA = "white", colorNULL = "gray", title = "Persons Over 64\nPer Square Kilometer",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd", convert2density = TRUE, area = "SqKm") +
  tm_shape(vt_tracts_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "Persons Over 64\nin Highest\nTransportation\nCost Burden\nBlock Groups", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_HiBurdenO64.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_HiBurdenO64.png")
```

```{r tabHiBurdenO64, fig.align = "center", fig.cap="Persons over 64 by municipality in Census Block Groups with the highest transportation cost burden (i.e., transportation cost burden in the 80th percentile for the state)."}
HiBurden %>%
  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)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Over 64` > 0 & `Pct of Over 64` > 0) %>%
  mutate(`Pct of Over 64` = if_else(`Pct of Over 64` > 100, 100, `Pct of Over 64`)) %>% 
  mutate(`Pct of Over 64` = paste0(`Pct of Over 64`,"%")) %>% 
  arrange(desc(`Over 64`)) %>% 
  filter(`Over 64` > 0) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Persons Over 64 in Highest Transportation Cost Burden Block Groups", align = "r", col.names = c(names(.)[1:2],"Number of Over 64 in Block Groups","Pct of Over 64 in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapHiBurdenLTHS, fig.align = "center", fig.cap="Adults without a high school diploma per square kilometer in Census Block Groups with the highest transportation cost burden (i.e., transportation cost burden in the 80th percentile for the state)."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(HiBurden, unit = "mi") +
  tm_fill(col = "lthsE", alpha = 0.9, style = "fisher",
          colorNA = "white", colorNULL = "gray", title = "Adults without a\nHigh School\nDiploma\nPer Square Kilometer",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd", convert2density = TRUE, area = "SqKm") +
  tm_shape(vt_tracts_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "Adults without a\nHigh School\nDiploma\nin Highest\nTransportation\nCost Burden\nBlock Groups", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_HiBurdenLTHS.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_HiBurdenLTHS.png")
```

```{r tabHiBurdenLTHS, fig.align = "center", fig.cap="Adults without a high school diploma by municipality in Census Block Groups with the highest transportation cost burden (i.e., transportation cost burden in the 80th percentile for the state)."}
HiBurden %>%
  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)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`No HS Dip` > 0 & `Pct of No HS Dip` > 0) %>%
  mutate(`Pct of No HS Dip` = if_else(`Pct of No HS Dip` > 100, 100, `Pct of No HS Dip`)) %>% 
  mutate(`Pct of No HS Dip` = paste0(`Pct of No HS Dip`,"%")) %>% 
  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 Transportation Cost Burden Block Groups", align = "r", col.names = c(names(.)[1:2],"Number of No HS Dip in Block Groups","Pct of No HS Dip in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```



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

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

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

[^headways]: TransitWiki. Headway. https://www.transitwiki.org/TransitWiki/index.php/Headway

[^mindGap]: Mind the Gap: Best practices in bus headway management. https://www.metro-magazine.com/10002870/mind-the-gap-best-practices-in-bus-headway-management

[^headwayReal]: Note that headways analyzed here are based on static schedules provided by the transit agency through the General Transit Feed System (GTFS). Scheduled headways likely differ from actual vehicle frequency due to traffic, dispatch management, accidents or breakdowns, and other issues. For an analysis of the difference between headways based on static schedules and actually observed headways, see Nate Wessel and Steven Farber. 2019. "On the accuracy of schedule-based GTFS for measuring accessibility." *The Journal of Transport and Land Use* 12(1):475-500. 

[^GTFSStaticOverview]: See Google Transit APIs. "GTFS Static Overview." https://developers.google.com/transit/gtfs

[^WalkCardio]: See Nicholas A. Howell, Jack V. Tu, Rahim Moineddin, Anna Chu, and Gillian L. Booth. 2019. "Association Between Neighborhood Walkability and Predicted 10‐Year Cardiovascular Disease Risk: The CANHEART (Cardiovascular Health in Ambulatory Care Research Team) Cohort." *Journal of the American Heart Association* 8(21). https://doi.org/10.1161/JAHA.119.013146; also see  Rachel C. Colley, Tanya Christidis, Isabelle Michaud, Michael Tjepkema and Nancy A. Ross. 2019. "The association between walkable neighbourhoods and physical activity across the lifespan." *Health Reports* 30(9).  http://dx.doi.org/10.25318/82-003-x201900900001-eng; see also Hao Wang and Yuqi Yang. 2019. "Neighbourhood walkability: A review and bibliometric analysis." *Cities* 93:43-61. https://doi.org/10.1016/j.cities.2019.04.015

[^Walkability]: See US EPA National Walkability Index. https://www.epa.gov/smartgrowth/smart-location-mapping#walkability; Also see Kathleen B. Watson, Geoffrey P. Whitfield, John V. Thomas, David Berrigan, Janet E. Fulton, and Susan A. Carlson. 2020. "Associations between the National Walkability Index and walking among US Adults — National Health Interview Survey, 2015." *Preventive Medicine* 137. https://doi.org/10.1016/j.ypmed.2020.106122.

[^Expenditures]: See U.S. Bureau of Labor Statistics. 2020. "Consumer Expenditures in 2018." BLS Reports May 2020.  https://www.bls.gov/opub/reports/consumer-expenditures/2018/home.htm; see also Thomas W. Sanchez, Carrie Makarewicz, Peter Haas, and Casey J. Dawkins. 2006. "Transportation costs, inequities, and tradeoffs." https://www.researchgate.net/profile/Thomas_Sanchez/publication/278017364_Transportation_Costs_Inequities_and_Tradeoffs/links/5578577e08aeb6d8c01f1382/Transportation-Costs-Inequities-and-Tradeoffs.pdf

[^CostBurden]: See Thomas W. Sanchez, Carrie Makarewicz, Peter Haas, and Casey J. Dawkins. 2006. "Transportation costs, inequities, and tradeoffs." https://www.researchgate.net/profile/Thomas_Sanchez/publication/278017364_Transportation_Costs_Inequities_and_Tradeoffs/links/5578577e08aeb6d8c01f1382/Transportation-Costs-Inequities-and-Tradeoffs.pdf

[^LAIabout]: See U.S. Department of Housing and Urban Development - HUD Exchange. "About the Location Affordability Index."  https://www.hudexchange.info/programs/location-affordability-index/about/

[^TransitDistance]: See C.K. Ayvalik and C.J. Khisty. 2002. "Heuristic analysis of impacts of commuter rail station consolidation on pedestrian access." *Transportation Research Record* 1793:47–54; E.A. Beimborn, M.J. Greenwald, and X. Jin. 2003. "Accessibility, connectivity, and captivity: impacts on transit choice." *Transportation Research Record* 1835: 1–9; S. Biba. 2010. "A new method for determining the population with walking access to transit." *International Journal of Geographical Information Science* 24(3): 347-364; J.E. Burkhardt. 2003. "Critical measures of transit service quality in the eyes of older travelers." *Transportation Research Record* 1835: 84–92; J. Dill. 2003. "Transit use and proximity to rail: results from large employment sites in the San Francisco, California, bay area." *Transportation Research Record* 1835:19–24; J.E. Evans, V. Perincherry, and G.B.I. Douglas. 1997. "Transit friendliness factor: approach to quantifying transit access environment in a transportation planning model." *Transportation Research Record* 1604: 32–39; S. Hsiao et al., 1997. "Use of geographic information system for analysis of transit pedestrian access." *Transportation Research Record* 1604: 50–59; D.R. Loutzenheiser. 1997. "Pedestrian access to transit: model of walk trips and their design and urban form determinants around bar area rapid transit stations." *Transportation Research Record* 1604: 40–49; W.A. O’Neill, R.D. Ramsey, and J. Chou. 1992. "Analysis of transit service areas using geographic information systems." *Transportation Research Record* 1364: 131–138; S.  O’Sullivan and J. Morrall. 1996. "Walking distances to and from light-rail transit stations." *Transportation Research Record* 1538: 19–26; C.G. Phillips  and H.R. Edwards. 2002. "Socioeconomic, community-based approach for developing integrated mass transit systems application to city of Baltimore, Maryland." *Transportation Research Record* 1797: 71–79; F. Zhao et al., 2003. "Forecasting transit walk accessibility – regression model alternative to buffer method." *Transportation Research Record* 1835: 34–41.


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

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

* General Transit Feed Specification (GTFS)
* MRLC National Land Cover Database (NLCD)
* US EPA Walkability Index
* US Department of Housing and Urban Development Location Affordability Index
* American Community Survey 5-year Estimates
  + Population demographics

## General Transit Feed Specification (GTFS) {-}

The transit route linework and transit stop locations for Vermont were derived from General Transit Feed Specification (GTFS) data.  

The General Transit Feed Specification (GTFS) is a a common format for public transportation schedules and associated geographic information. Public transit agencies increasingly publish their transit data electronically via GTFS "feeds" primarily to allow software developers to create applications which can tap into these feeds and provide users with detailed or timely transit schedules or route planning.[^GTFSStaticOverview] GTFS is also increasingly used by transportation researchers because it provides an efficient way to acquire detailed scheduling and routing information for transit networks. While GTFS feeds can provide real-time data on vehicle movements and locations, most data is provided as static schedules, and vehicle frequencies and headways are inferred from the schedule. 

GTFS data for the Vermont transit agencies were acquired from two sources: 

* Vermont GTFS data at https://vermont-gtfs.org/
* Open Mobility Data at https://transitfeeds.com/l/415-vermont-usa 

GTFS data was processed using the `tidytransit` package in R. Schedules for all routes running Monday through Friday, 6am to 9pm, were extracted. Headways for each transit stop on these routes was then calculated based on route schedules, and these headways were also averaged for each route.

To calculate headways experienced by different population groups, transit stop points with calculated headways were intersected with the developed portions of Census Block Groups and that were within 400m, and 4800m of these stops for buses and commuter rail, respectively. Headways for these intersecting stops were then aggregated to provide mean headways for their intersecting Block Groups. Populations within those intersecting Block Groups were assumed to be subject to those mean headways. 

A resident is considered to have access if she is within one‐quarter mile of a bus stop, one‐half mile of a rapid transit (i.e. subway or light rail) station, and up to three miles of a commuter rail station or park-and-ride lot. This geographic definition of access is consistent with FTA Circular 4702.1B Title VI Requirements and Guidelines for Federal Transit Administration Recipients, Chapter IV‐14, as well as decades of research on transit behavior, and the judgment of most planners and transit researchers.[^TransitDistance] 

To identify the populations with access to transit, 400 meter, 800 meter, and 4800 meter buffers were generated around all bus, rapid transit, and commuter rail stops and park-and-ride lots, respectively. The transit buffer was intersected with the developed portions of Census Block Groups. Developed portions of Census Block Groups were identified based upon the National Land Cover Database (NLCD), using a process of areal apportionment (see below). The population with access to transit was calculated as the product of the  areal proportion of the intersecting buffer and developed Block Group polygons:

*Population with transit access = Proportion of developed Block Group Intersection x Population of developed Block Group*

For example, if 10% of the developed area of a Census Block Group intersected/overlapped with the 400m buffer around transit stops, it was assumed that 10% of the population has access to those transit stops. Assuming a population of 100 people in the developed portion of the Block Group, this would mean *100 x .10 = 10* people would have access to transit via those stops. This transit access analysis was done in R.


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

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

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

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

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

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

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


## US EPA Walkability Index {-}

The National Walkability Index is a nationwide geographic data resource produced by the U.S. Environmental Protection Agency that ranks Census Block Groups according to their relative walkability.[^Walkability] Walkability was modeled based on characteristics of the built environment that influence the likelihood of walking being used as a mode of travel. These characteristics of Block Groups include:

* mix of employment types (greater mix = more walkable)
* amount or density of occupied housing (higher density = more walkable)
* street intersection density (higher density = more walkable)
* proportion of workers who carpool (more carpooling = more walkable)

The Walkability Index was downloaded as a shapefile from the US EPA's Smart Location Mapping website at https://www.epa.gov/smartgrowth/smart-location-mapping#walkability. To compare walkability index values with disabled populations and households without a car, Block Group walkability index values were aggregated to the overlying Census Tract as simple averages. Population-weighted averages for the walkability index were also computed for Block Groups and Tracts to compare population groups. 


## US Department of Housing and Urban Development Location Affordability Index {-}

The Location Affordability Index (LAI) is a nationwide geographic data resource developed by the U.S. Department of Housing and Urban Development (HUD) in collaboration with the U.S. Department of Transportation under the federal Partnership for Sustainable Communities.[^LAIabout] The LAI was developed as a way of integrating both housing and transportation costs into estimates of the affordability of specific neighborhoods and cities. The LAI provides estimated total expenses from transportation, whether transit or motor vehicle, as well as transportation cost as a percentage of household income, for eight different household types — which vary by household income, size, and number of commuters — at the Census Tract level. For this analysis Location Affordability Index Model Version 3 (LAIM Version 3 or LAIM3), released in March 2019, was used. LAIM Version 3 estimates housing and transportation costs for eight different household profiles, to focus on the impact of the built environment on these costs by holding demographic characteristics constant. These household profiles are described in Table \@ref(tab:tabLAIprofiles). The LAI can be interpreted as the percentage of the specified household profile’s income that would be needed to cover housing, transportation, or combined costs if that household type were to live in the specified Census Tract. So, for example, a household profile (such as single-parent family) with a housing LAI of 40% would be expected to pay 40% of their income for housing if they lived in the specified location. Details on how housing and transportation costs were modeled are descrbed in the Location Affordability Index Version 3 Data and Methodology document. 

```{r tabLAIprofiles}
data.frame(hhProfile = c("Median-Income Family", "Very Low-Income Individual", "Working Individual", "Single Professional", "Retired Couple", "Single-Parent Family", "Moderate-Income Family", "Dual-Professional Family"),
hhhIncome = c("MHHI", "National poverty line", "50% of MHHI", "135% of MHHI", "80% of MHHI", "50% of MHHI", "80% of MHHI", "150% of MHHI"),
hhSize = c(4,1,1,1,2,3,3,4),
Commuters = c(2,1,1,1,0,1,1,2)) %>% 
  kableExtra::kable(longtable = T, booktabs = T, 
                    caption = "Location Affordability Index Household Profiles", 
                    align = "l",
                    col.names = c("Household Profile","Income","Household Size","Commuters")) %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header","striped")) %>% 
  kableExtra::footnote(., general = "MHHI = Median household income for a given area (CBSA or County)")
```

For the purposes of this analysis, transportation costs as a percentage of income for moderate-income households (three-member households with one commuter making 80% or less of area median household income) was used as a proxy to represent transportation cost burden for all areas. LAI data were also downscaled to the Block Group level for the purpose of demographic analyses and comparisons. Census Block Groups belonging to a given Census Tract were assigned the same value as the overlying Tract. 


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

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

## Population-weighted averages {-}

Wherever feasible, averages are reported as population-weighted averages. 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 experience longer headways than the general or Total Population. Consider the table below.

```{r wmean1}
data.frame(BG = c("BG1","BG2","BG3","BG4"),
           Headway = 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 Headways", 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 headways, are each reported by Block Group. Since each Block Group is associated with one headway value, we might assume (incorrectly) that everyone is equally exposed to the average headways 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 headway value is associated with a different number of people. Do People of Color occupy Block Groups with higher headways than the simple average would indicate?

To calculate the population-weighted average headway for People of Color, the number of People of Color in each Block Group is used as a 'weight'. The headway 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"),
           Headway = c(65,75,85,70),
           Minority = c(5,10,20,10),
           NonMinority = c(20,12,5,10)) %>% 
  mutate(TotalPop = Minority + NonMinority,
         HeadwayxPOC = Headway*Minority) %>% 
  bind_rows(summarize_all(., funs(if(is.numeric(.)) sum(.) else "Total"))) %>% 
  kableExtra::kable(longtable = T, booktabs = T, 
                    caption = "Block Group Populations and Headway", 
                    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., POC x Headway) is then divided by the sum of the weights (i.e., total People of Color), so that 3475/45 = `r round(3475/45,2)`. The result is a weighted average headway 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 headways 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"),
           Headway = c(65,75,85,70),
           Minority = c(5,10,20,10),
           NonMinority = c(20,12,5,10)) %>% 
  mutate(TotalPop = Minority + NonMinority,
         TempxTotalPop = Headway*TotalPop) %>% 
  bind_rows(summarize_all(., funs(if(is.numeric(.)) sum(.) else "Total"))) %>% 
  kableExtra::kable(longtable = T, booktabs = T, 
                    caption = "Block Group Populations and Headway", 
                    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 headway is 6800/92 = `r round(6800/92,2)`. Thus we can see that People of Color experience a higher population-weighted average headway than the general or total population. 


\pagebreak
# Appendix B: Supplementary Figures {-}

## No Transit Access and Priority Populations {-}
```{r mapTransitAccessO64, fig.align = "center", fig.cap="Persons over 64 without reasonable access to any form of transit by Census Block Group."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") +
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(vt_blkgrps_sf_noTransit_simple, unit = "mi") +
  tm_fill(col = "NewOver64", style = "fisher",
          colorNA = "white", colorNULL = "gray", title = "Number of\nPersons Over 64",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd") +
  tm_shape(vt_tracts_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5,
               position = c(0.6,0.005)) +
  tm_layout(title = "Persons Over 64\nwithout Reasonable\nAccess to any\nMode of Transit",
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_transitO64access.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_transitO64access.png")
```

```{r tabTransitAccessO64, fig.align = "center", fig.cap="Persons over 64 in Census Block Groups without reasonable access to any mode of transit."}
vt_blkgrps_sf_noTransit %>%
  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(NewOver64,na.rm = TRUE),
            `Pct of Over 64` = round(sum(NewOver64)/max(TotalOver64)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Pct of Over 64` > 0) %>% 
  mutate(`Pct of Over 64` = if_else(`Pct of Over 64` > 100,100,`Pct of Over 64`)) %>% 
  mutate(`Pct of Over 64` = paste0(`Pct of Over 64`,"%")) %>% 
  arrange(desc(`Over 64`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Persons Over 64 without Access to any Mode of Transit", align = "r", col.names = c(names(.)[1:2],"Number of Over 64 in Census Block Groups","Pct of Over 64 in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped")) 
```

```{r mapTransitAccessLTHS, fig.align = "center", fig.cap="Adults without a high school diploma without reasonable access to any form of transit by Census Block Group."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") +
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(vt_blkgrps_sf_noTransit_simple, unit = "mi") +
  tm_fill(col = "NewLths", style = "fisher",
          colorNA = "white", colorNULL = "gray", title = "Number of\nNo HS Dip",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd") +
  tm_shape(vt_tracts_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5,
               position = c(0.6,0.005)) +
  tm_layout(title = "Adults without a\nHigh School\nDiploma\nwithout Reasonable\nAccess to any\nMode of Transit",
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_transitLTHSaccess.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_transitLTHSaccess.png")
```

```{r tabTransitAccessLTHS, fig.align = "center", fig.cap="Adults without a high school diploma in Census Block Groups without reasonable access to any mode of transit."}
vt_blkgrps_sf_noTransit %>%
  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(NewLths,na.rm = TRUE),
            `Pct of No HS Dip` = round(sum(NewLths)/max(TotalNoHS)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Pct of No HS Dip` > 0) %>% 
  mutate(`Pct of No HS Dip` = if_else(`Pct of No HS Dip` > 100,100,`Pct of No HS Dip`)) %>% 
  mutate(`Pct of No HS Dip` = paste0(`Pct of No HS Dip`,"%")) %>% 
  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 without Access to any Mode of Transit", align = "r", col.names = c(names(.)[1:2],"Number of No HS Dip in Census Block Groups","Pct of No HS Dip in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped")) 
```

```{r mapTransitAccessDisabled, fig.align = "center", fig.cap="Disabled persons without reasonable access to any form of transit by Census Block Group."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") +
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(vt_tracts_sf_noTransit_simple, unit = "mi") +
  tm_fill(col = "NewDisabled", style = "fisher",
          colorNA = "white", colorNULL = "gray", title = "Number of\nDisabled Persons",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd") +
  tm_shape(vt_tracts_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5,
               position = c(0.6,0.005)) +
  tm_layout(title = "Disabled Persons\nwithout Reasonable\nAccess to any\nMode of Transit",
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_transitDisabledaccess.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_transitDisabledaccess.png")
```

```{r tabTransitAccessDisabled, fig.align = "center", fig.cap="Disabled persons in Census Tracts without reasonable access to any mode of transit."}
vt_tracts_sf_noTransit %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Tracts` = n(),
            `Disabled` = sum(NewDisabled,na.rm = TRUE),
            `Pct of Disabled` = round(sum(NewDisabled)/max(TotalDisabled)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Pct of Disabled` > 0) %>% 
  mutate(`Pct of Disabled` = if_else(`Pct of Disabled` > 100,100,`Pct of Disabled`)) %>% 
  mutate(`Pct of Disabled` = paste0(`Pct of Disabled`,"%")) %>% 
  arrange(desc(`Disabled`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Disabled Persons without Access to any Mode of Transit", align = "r", col.names = c(names(.)[1:2],"Number of Disabled in Census Tracts","Pct of Disabled in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped")) 
```

```{r mapTransitAccessINC, fig.align = "center", fig.cap="Low income persons without reasonable access to any form of transit by Census Block Group."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") +
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(vt_blkgrps_sf_noTransit_simple, unit = "mi") +
  tm_fill(col = "NewPov", style = "fisher",
          colorNA = "white", colorNULL = "gray", title = "Number of\nLow Income Persons",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd") +
  tm_shape(vt_tracts_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5,
               position = c(0.6,0.005)) +
  tm_layout(title = "Low Income\nPersons\nwithout Reasonable\nAccess to any\nMode of Transit",
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_transitINCaccess.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_transitINCaccess.png")
```

```{r tabTransitAccessINC, fig.align = "center", fig.cap="Low income persons in Census Block Groups without reasonable access to any mode of transit."}
vt_blkgrps_sf_noTransit %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Low Income` = sum(NewPov,na.rm = TRUE),
            `Pct of Low Income` = round(sum(NewPov)/max(TotalLowInc)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Pct of Low Income` > 0) %>% 
  mutate(`Pct of Low Income` = if_else(`Pct of Low Income` > 100,100,`Pct of Low Income`)) %>% 
  mutate(`Pct of Low Income` = paste0(`Pct of Low Income`,"%")) %>% 
  arrange(desc(`Low Income`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Low Income Persons without Access to any Mode of Transit", align = "r", col.names = c(names(.)[1:2],"Number of Low Income in Census Block Groups","Pct of Low Income in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped")) 
```

```{r mapTransitAccessMIN, fig.align = "center", fig.cap="People of Color without reasonable access to any form of transit by Census Block Group."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") +
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(vt_blkgrps_sf_noTransit_simple, unit = "mi") +
  tm_fill(col = "NewMinority", style = "fisher",
          colorNA = "white", colorNULL = "gray", title = "Number of\nPeople of Color",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd") +
  tm_shape(vt_tracts_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5,
               position = c(0.6,0.005)) +
  tm_layout(title = "People of Color\nwithout Reasonable\nAccess to any\nMode of Transit",
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_transitMINaccess.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_transitMINaccess.png")
```

```{r tabTransitAccessMIN, fig.align = "center", fig.cap="People of Color in Census Block Groups without reasonable access to any mode of transit."}
vt_blkgrps_sf_noTransit %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Minority` = sum(NewMinority,na.rm = TRUE),
            `Pct of Minority` = round(sum(NewMinority)/max(TotalMin)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Pct of Minority` > 0) %>% 
  mutate(`Pct of Minority` = if_else(`Pct of Minority` > 100,100,`Pct of Minority`)) %>% 
  mutate(`Pct of Minority` = paste0(`Pct of Minority`,"%")) %>% 
  arrange(desc(`Minority`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "People of Color without Access to any Mode of Transit", align = "r", col.names = c(names(.)[1:2],"Number of POC in Census Block Groups","Pct of POC in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped")) 
```

```{r mapTransitAccessLEH, fig.align = "center", fig.cap="Limited English speaking households without reasonable access to any form of transit by Census Block Group."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") +
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(vt_blkgrps_sf_noTransit_simple, unit = "mi") +
  tm_fill(col = "NewLths", style = "fisher",
          colorNA = "white", colorNULL = "gray", title = "Number of\nLimited English\nHouseholds",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd") +
  tm_shape(vt_tracts_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5,
               position = c(0.6,0.005)) +
  tm_layout(title = "Limited English\nHouseholds\nwithout Reasonable\nAccess to any\nMode of Transit",
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_transitLEHaccess.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_transitLEHaccess.png")
```

```{r tabTransitAccessLEH, fig.align = "center", fig.cap="Limited English speaking households in Census Block Groups without reasonable access to any mode of transit."}
vt_blkgrps_sf_noTransit %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Limited English HH` = sum(NewLths,na.rm = TRUE),
            `Pct of Limited English HH` = round(sum(NewLths)/max(TotalLEH)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Pct of Limited English HH` > 0) %>% 
  mutate(`Pct of Limited English HH` = if_else(`Pct of Limited English HH` > 100,100,`Pct of Limited English HH`)) %>% 
  mutate(`Pct of Limited English HH` = paste0(`Pct of Limited English HH`,"%")) %>% 
  arrange(desc(`Limited English HH`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Limited English Speaking Households without Access to any Mode of Transit", align = "r", col.names = c(names(.)[1:2],"Number of Limited English HH in Census Block Groups","Pct of Limited English HH in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped")) 
```

```{r mapTransitAccessNoCarHH, fig.align = "center", fig.cap="Households without a car without reasonable access to any form of transit by Census Block Group."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") +
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(vt_tracts_sf_noTransit_simple, unit = "mi") +
  tm_fill(col = "NewHHNoCar", style = "fisher",
          colorNA = "white", colorNULL = "gray", title = "Number of\nNo Car Households",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd") +
  tm_shape(vt_tracts_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5,
               position = c(0.6,0.005)) +
  tm_layout(title = "No Car Households\nwithout Reasonable\nAccess to any\nMode of Transit",
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_transitNoCarHHaccess.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_transitNoCarHHaccess.png")
```

```{r tabTransitAccessNoCarHH, fig.align = "center", fig.cap="Households without a car in Census Block Groups without reasonable access to any mode of transit."}
vt_tracts_sf_noTransit %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Tracts` = n(),
            `No Car Households` = sum(NewHHNoCar,na.rm = TRUE),
            `Pct of No Car Households` = round(sum(NewHHNoCar)/max(TotalNoCar)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Pct of No Car Households` > 0) %>% 
  mutate(`Pct of No Car Households` = if_else(`Pct of No Car Households` > 100,100,`Pct of No Car Households`)) %>% 
  mutate(`Pct of No Car Households` = paste0(`Pct of No Car Households`,"%")) %>% 
  arrange(desc(`No Car Households`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "No Car Households without Access to any Mode of Transit", align = "r", col.names = c(names(.)[1:2],"Number of No Car Households in Census Tracts","Pct of No Car Households in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped")) 
```


## Excessive Headways and Priority Populations {-}
```{r mapHeadwayDisabledbus, fig.align = "center", fig.cap="Disabled persons within 400m or 1/4 mile of transit stops in which the mean headway for these stops is above the 80th percentile for the state."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") +
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(vt_towns_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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(headway80thtracts_bus) +
  tm_fill(col = "NewDisabled", style = "fisher",
          colorNA = "white", colorNULL = "gray",
          title = "Number of\nDisabled Persons",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd") +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5,
               position = c(0.6,0.005)) +
  # tm_add_legend(type = "line", lwd=1, alpha = 0.3, labels = "Bus Route", col = "blue") +
  tm_layout(title = paste("Disabled Persons","with Bus Headways",paste0("Over ",round(quantile(vt_routes_sf_Bus$mean_headways,0.8),0)," Minutes"),sep = "\n"),
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_headwayDisabledbus.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_headwayDisabledbus.png")
```

```{r tabHeadwayDisabledbus, fig.align = "center", fig.cap="Disabled persons within 400m or 1/4 mile of transit stops in which the mean headway for these stops is above the 80th percentile for the state."}
headway80thtracts_bus %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Tracts` = n(),
            `Disabled` = sum(NewDisabled,na.rm = TRUE),
            `Pct of Disabled` = round(sum(NewDisabled)/max(TotalDisabled)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Pct of Disabled` > 0) %>% 
  mutate(`Pct of Disabled` = if_else(`Pct of Disabled` > 100, 100, `Pct of Disabled`)) %>% 
  mutate(`Pct of Disabled` = paste0(`Pct of Disabled`,"%")) %>% 
  arrange(desc(`Disabled`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = paste("Disabled Persons","with Average Headways",paste0("Over ",round(quantile(vt_routes_sf_Bus$mean_headways,0.8),0)," Minutes"),sep = "\n"), align = "r", col.names = c(names(.)[1:2],"Number of Disabled in Census Tracts","Pct of Disabled in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapHeadwayNoHSDipbus, fig.align = "center", fig.cap="Adults without a high school diploma within 400m or 1/4 mile of transit stops in which the mean headway for these stops is above the 80th percentile for the state."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") +
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(vt_towns_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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(headway80thblkgrps_bus) +
  tm_fill(col = "NewLths", style = "fisher",
          colorNA = "white", colorNULL = "gray",
          title = "Number of\nNo HS Dip",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd") +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5,
               position = c(0.6,0.005)) +
  # tm_add_legend(type = "line", lwd=1, alpha = 0.3, labels = "Bus Route", col = "blue") +
  tm_layout(title = paste("Adults without a","High School","Diploma","with Bus Headways",paste0("Over ",round(quantile(vt_routes_sf_Bus$mean_headways,0.8),0)," Minutes"),sep = "\n"),
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_headwayNoHSDipbus.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_headwayNoHSDipbus.png")
```

```{r tabHeadwayNoHSDipbus, fig.align = "center", fig.cap="Adults without a high school diploma within 400m or 1/4 mile of transit stops in which the mean headway for these stops is above the 80th percentile for the state."}
headway80thblkgrps_bus %>%
  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(NewLths,na.rm = TRUE),
            `Pct of No HS Dip` = round(sum(NewLths)/max(TotalNoHS)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Pct of No HS Dip` > 0) %>% 
  mutate(`Pct of No HS Dip` = if_else(`Pct of No HS Dip` > 100, 100, `Pct of No HS Dip`)) %>% 
  mutate(`Pct of No HS Dip` = paste0(`Pct of No HS Dip`,"%")) %>% 
  arrange(desc(`No HS Dip`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = paste("Adults without a","High School","Diploma","with Bus Headways",paste0("Over ", round(quantile(vt_routes_sf_Bus$mean_headways,0.8),0), " Minutes"),sep = "\n"), align = "r", col.names = c(names(.)[1:2],"Number of No HS Dip in Block Groups","Pct of No HS Dip in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapHeadwayNoCarHHbus, fig.align = "center", fig.cap="Households without a car within 400m or 1/4 mile of transit stops in which the mean headway for these stops is above the 80th percentile for the state."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") +
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(vt_towns_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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(headway80thtracts_bus) +
  tm_fill(col = "NewHHNoCar", style = "fisher",
          colorNA = "white", colorNULL = "gray",
          title = "Number of\nNo Car Households",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd") +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5,
               position = c(0.6,0.005)) +
  # tm_add_legend(type = "line", lwd=1, alpha = 0.3, labels = "Bus Route", col = "blue") +
  tm_layout(title = paste("No Car Households","with Bus Headways",paste0("Over ",round(quantile(vt_routes_sf_Bus$mean_headways,0.8),0)," Minutes"),sep = "\n"),
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_headwayNoCarHHbus.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_headwayNoCarHHbus.png")
```

```{r tabHeadwayNoCarHHbus, fig.align = "center", fig.cap="Households without a car within 400m or 1/4 mile of transit stops in which the mean headway for these stops is above the 80th percentile for the state."}
headway80thtracts_bus %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `No Car Households` = sum(NewHHNoCar,na.rm = TRUE),
            `Pct of No Car Households` = round(sum(NewHHNoCar)/max(TotalNoCar)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Pct of No Car Households` > 0) %>% 
  mutate(`Pct of No Car Households` = if_else(`Pct of No Car Households` > 100, 100, `Pct of No Car Households`)) %>% 
  mutate(`Pct of No Car Households` = paste0(`Pct of No Car Households`,"%")) %>% 
  arrange(desc(`No Car Households`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = paste("No Car Households","with Average Headways",paste0("Over ",round(quantile(vt_routes_sf_Bus$mean_headways,0.8),0)," Minutes"),sep = "\n"), align = "r", col.names = c(names(.)[1:2],"Number of No Car Households in Block Groups","Pct of No Car Households in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapHeadwayLowIncomebus, fig.align = "center", fig.cap="Low income persons within 400m or 1/4 mile of transit stops in which the mean headway for these stops is above the 80th percentile for the state."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") +
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(vt_towns_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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(headway80thblkgrps_bus) +
  tm_fill(col = "NewPov", style = "fisher",
          colorNA = "white", colorNULL = "gray",
          title = "Number of\nLow Income Persons",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd") +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5,
               position = c(0.6,0.005)) +
  # tm_add_legend(type = "line", lwd=1, alpha = 0.3, labels = "Bus Route", col = "blue") +
  tm_layout(title = paste("Low Income","Persons","with Bus Headways",paste0("Over ",round(quantile(vt_routes_sf_Bus$mean_headways,0.8),0)," Minutes"),sep = "\n"),
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_headwayLowIncomebus.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_headwayLowIncomebus.png")
```

```{r tabHeadwayLowIncomebus, fig.align = "center", fig.cap="Low income persons within 400m or 1/4 mile of transit stops in which the mean headway for these stops is above the 80th percentile for the state."}
headway80thblkgrps_bus %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Low Income` = sum(NewPov,na.rm = TRUE),
            `Pct of Low Income` = round(sum(NewPov)/max(TotalLowInc)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Pct of Low Income` > 0) %>% 
  mutate(`Pct of Low Income` = if_else(`Pct of Low Income` > 100, 100, `Pct of Low Income`)) %>% 
  mutate(`Pct of Low Income` = paste0(`Pct of Low Income`,"%")) %>% 
  arrange(desc(`Low Income`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = paste("Low Income Persons","with Average Headways",paste0("Over ",quantile(vt_routes_sf_Bus$mean_headways,0.8)," Minutes"),sep = "\n"), align = "r", col.names = c(names(.)[1:2],"Number of Low Income in Block Groups","Pct of Low Income in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapHeadwayMinoritybus, fig.align = "center", fig.cap="People of Color within 400m or 1/4 mile of transit stops in which the mean headway for these stops is above the 80th percentile for the state."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") +
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(vt_towns_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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(headway80thblkgrps_bus) +
  tm_fill(col = "NewMinority", style = "fisher",
          colorNA = "white", colorNULL = "gray",
          title = "Number of\nPeople of Color",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd") +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5,
               position = c(0.6,0.005)) +
  # tm_add_legend(type = "line", lwd=1, alpha = 0.3, labels = "Bus Route", col = "blue") +
  tm_layout(title = paste("People of Color","with Bus Headways",paste0("Over ",quantile(vt_routes_sf_Bus$mean_headways,0.8)," Minutes"),sep = "\n"),
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_headwayMinoritybus.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_headwayMinoritybus.png")
```

```{r tabHeadwayMinoritybus, fig.align = "center", fig.cap="People of Color within 400m or 1/4 mile of transit stops in which the mean headway for these stops is above the 80th percentile for the state."}
headway80thblkgrps_bus %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Minority` = sum(NewMinority,na.rm = TRUE),
            `Pct of Minority` = round(sum(NewMinority)/max(TotalMin)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Pct of Minority` > 0) %>% 
  mutate(`Pct of Minority` = if_else(`Pct of Minority` > 100, 100, `Pct of Minority`)) %>% 
  mutate(`Pct of Minority` = paste0(`Pct of Minority`,"%")) %>% 
  arrange(desc(`Minority`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = paste("People of Color","with Average Headways",paste0("Over ",quantile(vt_routes_sf_Bus$mean_headways,0.8)," Minutes"),sep = "\n"), align = "r", col.names = c(names(.)[1:2],"Number of POC in Block Groups","Pct of POC in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapHeadwayOver64bus, fig.align = "center", fig.cap="Persons over 64 within 400m or 1/4 mile of transit stops in which the mean headway for these stops is above the 80th percentile for the state."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") +
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(vt_towns_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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(headway80thblkgrps_bus) +
  tm_fill(col = "NewOver64", style = "fisher",
          colorNA = "white", colorNULL = "gray",
          title = "Number of\nPersons Over 64",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd") +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5,
               position = c(0.6,0.005)) +
  # tm_add_legend(type = "line", lwd=1, alpha = 0.3, labels = "Bus Route", col = "blue") +
  tm_layout(title = paste("Persons Over 64","with Bus Headways",paste0("Over ",round(quantile(vt_routes_sf_Bus$mean_headways,0.8),0)," Minutes"),sep = "\n"),
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_headwayOver64bus.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_headwayOver64bus.png")
```

```{r tabHeadwayOver64bus, fig.align = "center", fig.cap="Persons over 64 within 400m or 1/4 mile of transit stops in which the mean headway for these stops is above the 80th percentile for the state."}
headway80thblkgrps_bus %>%
  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(NewOver64,na.rm = TRUE),
            `Pct of Over 64` = round(sum(NewOver64)/max(TotalOver64)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Pct of Over 64` > 0) %>% 
  mutate(`Pct of Over 64` = if_else(`Pct of Over 64` > 100, 100, `Pct of Over 64`)) %>% 
  mutate(`Pct of Over 64` = paste0(`Pct of Over 64`,"%")) %>% 
  arrange(desc(`Over 64`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = paste("Persons Over 64","with Average Headways",paste0("Over ",quantile(vt_routes_sf_Bus$mean_headways,0.8)," Minutes"),sep = "\n"), align = "r", col.names = c(names(.)[1:2],"Number of Over 64 in Block Groups","Pct of Over 64 in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapHeadwayLEHbus, fig.align = "center", fig.cap="Limited English speaking households within 400m or 1/4 mile of transit stops in which the mean headway for these stops is above the 80th percentile for the state."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") +
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(vt_towns_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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(headway80thblkgrps_bus) +
  tm_fill(col = "NewEng_limit", style = "fisher",
          colorNA = "white", colorNULL = "gray",
          title = "Number of\nLimited English\nHouseholds",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd") +
  tm_shape(I89roadSegment) +
  tm_symbols(shape = I89, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black",
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5,
               position = c(0.6,0.005)) +
  # tm_add_legend(type = "line", lwd=1, alpha = 0.3, labels = "Bus Route", col = "blue") +
  tm_layout(title = paste("Limited English", "Households","with Bus Headways",paste0("Over ",round(quantile(vt_routes_sf_Bus$mean_headways,0.8),0)," Minutes"),sep = "\n"),
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_headwayLEHbus.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_headwayLEHbus.png")
```

```{r tabHeadwayLEHbus, fig.align = "center", fig.cap="Limited English speaking households within 400m or 1/4 mile of transit stops in which the mean headway for these stops is above the 80th percentile for the state."}
headway80thblkgrps_bus %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Under 18` = sum(NewEng_limit,na.rm = TRUE),
            `Pct of Under 18` = round(sum(NewEng_limit)/max(TotalLEH)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Pct of Under 18` > 0) %>% 
  mutate(`Pct of Under 18` = if_else(`Pct of Under 18` > 100, 100, `Pct of Under 18`)) %>% 
  mutate(`Pct of Under 18` = paste0(`Pct of Under 18`,"%")) %>% 
  arrange(desc(`Under 18`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = paste("Limited English Households","with Average Headways",paste0("Over ",quantile(vt_routes_sf_Bus$mean_headways,0.8)," Minutes"),sep = "\n"), align = "r", col.names = c(names(.)[1:2],"Number of Under 18 in Block Groups","Pct of Under 18 in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```


## Low Walkability and Priority Populations {-}
```{r mapLowWalkU5, fig.align = "center", fig.cap="Children under 5 per square kilometer in least walkable Census Block Groups (i.e., Walkability Index Score less than or equal to 5.8)."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(leastWalkable, unit = "mi") +
  tm_fill(col = "under5E", alpha = 0.9, style = "fisher",
          colorNA = "white", colorNULL = "gray", title = "Children Under 5\nPer Square Kilometer",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd", convert2density = TRUE, area = "SqKm") +
  tm_shape(vt_tracts_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "Children Under 5\nin Least Walkable\nBlock Groups", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_walkU5.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_walkU5.png")
```

```{r tabLowWalkU5, fig.align = "center", fig.cap="Children under 5 by municipality in least walkable Census Block Groups (i.e., Walkability Index Score less than or equal to 5.8)."}
leastWalkable %>%
  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)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Under 5` > 0 & `Pct of Under 5` > 0) %>% 
  mutate(`Pct of Under 5` = if_else(`Pct of Under 5` > 100, 100, `Pct of Under 5`)) %>% 
  mutate(`Pct of Under 5` = paste0(`Pct of Under 5`,"%")) %>% 
  arrange(desc(`Under 5`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Children Under 5 in Least Walkable Block Groups", align = "r", col.names = c(names(.)[1:2],"Number of Under 5 in Block Groups","Pct of Under 5 in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapLowWalkLTHS, fig.align = "center", fig.cap="Adults without a high school diploma per square kilometer in least walkable Census Block Groups (i.e., Walkability Index Score less than or equal to 5.8)."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(leastWalkable, unit = "mi") +
  tm_fill(col = "lthsE", alpha = 0.9, style = "fisher",
          colorNA = "white", colorNULL = "gray", title = "Adults without a\nHigh School\nDiploma\nPer Square Kilometer",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd", convert2density = TRUE, area = "SqKm") +
  tm_shape(vt_tracts_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "Adults without a\nHigh School\nDiploma\nin Least Walkable\nBlock Groups", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_walkLTHS.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_walkLTHS.png")
```

```{r tabLowWalkLTHS, fig.align = "center", fig.cap="Adults without a high school diploma by municipality in least walkable Census Block Groups (i.e., Walkability Index Score less than or equal to 5.8)."}
leastWalkable %>%
  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)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`No HS Dip` > 0 & `Pct of No HS Dip` > 0) %>% 
  mutate(`Pct of No HS Dip` = if_else(`Pct of No HS Dip` > 100, 100, `Pct of No HS Dip`)) %>% 
  mutate(`Pct of No HS Dip` = paste0(`Pct of No HS Dip`,"%")) %>% 
  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 Least Walkable Block Groups", align = "r", col.names = c(names(.)[1:2],"Number of No HS Dip in Block Groups","Pct of No HS Dip in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapLowWalkDisabled, fig.align = "center", fig.cap="Disabled persons per square kilometer in least walkable Census Block Groups (i.e., Walkability Index Score less than or equal to 5.8)."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(leastWalkableTract, unit = "mi") +
  tm_fill(col = "disabledOver18E", alpha = 0.9, style = "fisher",
          colorNA = "white", colorNULL = "gray", title = "Disabled Persons\nPer Square Kilometer",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd", convert2density = TRUE, area = "SqKm") +
  tm_shape(vt_tracts_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "Disabled Persons\nin Least Walkable\nBlock Groups", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_walkDisabled.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_walkDisabled.png")
```

```{r tabLowWalkDisabled, fig.align = "center", fig.cap="Disabled persons by municipality in least walkable Census Block Groups (i.e., Walkability Index Score less than or equal to 5.8)."}
leastWalkableTract %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Disabled` = sum(disabledOver18E,na.rm = TRUE),
            `Pct of Disabled` = round(sum(disabledOver18E)/max(TotalDisabled)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Disabled` > 0 & `Pct of Disabled` > 0) %>% 
  mutate(`Pct of Disabled` = if_else(`Pct of Disabled` > 100, 100, `Pct of Disabled`)) %>% 
  mutate(`Pct of Disabled` = paste0(`Pct of Disabled`,"%")) %>% 
  arrange(desc(`Disabled`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Disabled Persons in Least Walkable Block Groups", align = "r", col.names = c(names(.)[1:2],"Number of Disabled in Block Groups","Pct of Disabled in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapLowWalkINC, fig.align = "center", fig.cap="Low income persons per square kilometer in least walkable Census Block Groups (i.e., Walkability Index Score less than or equal to 5.8)."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(leastWalkable, unit = "mi") +
  tm_fill(col = "num2povE", alpha = 0.9, style = "fisher",
          colorNA = "white", colorNULL = "gray", title = "Low Income\nPersons\nPer Square Kilometer",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd", convert2density = TRUE, area = "SqKm") +
  tm_shape(vt_tracts_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "Low Income\nPersons\nin Least Walkable\nBlock Groups", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_walkINC.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_walkINC.png")
```

```{r tabLowWalkINC, fig.align = "center", fig.cap="Low income persons by municipality in least walkable Census Block Groups (i.e., Walkability Index Score less than or equal to 5.8)."}
leastWalkable %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Low Income` = sum(num2povE,na.rm = TRUE),
            `Pct of Low Income` = round(sum(num2povE)/max(TotalLowInc)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Low Income` > 0 & `Pct of Low Income` > 0) %>% 
  mutate(`Pct of Low Income` = if_else(`Pct of Low Income` > 100, 100, `Pct of Low Income`)) %>% 
  mutate(`Pct of Low Income` = paste0(`Pct of Low Income`,"%")) %>% 
  arrange(desc(`Low Income`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Low Income Persons in Least Walkable Block Groups", align = "r", col.names = c(names(.)[1:2],"Number of Low Income in Block Groups","Pct of Low Income in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapLowWalkMIN, fig.align = "center", fig.cap="People of Color per square kilometer in least walkable Census Block Groups (i.e., Walkability Index Score less than or equal to 5.8)."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(leastWalkable, unit = "mi") +
  tm_fill(col = "minorityE", alpha = 0.9, style = "fisher",
          colorNA = "white", colorNULL = "gray", title = "People of Color\nPer Square Kilometer",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd", convert2density = TRUE, area = "SqKm") +
  tm_shape(vt_tracts_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "People of Color\nin Least Walkable\nBlock Groups", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_walkMIN.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_walkMIN.png")
```

```{r tabLowWalkMIN, fig.align = "center", fig.cap="People of Color by municipality in least walkable Census Block Groups (i.e., Walkability Index Score less than or equal to 5.8)."}
leastWalkable %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Minority` = sum(minorityE,na.rm = TRUE),
            `Pct of Minority` = round(sum(minorityE)/max(TotalMin)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Minority` > 0 & `Pct of Minority` > 0) %>% 
  mutate(`Pct of Minority` = if_else(`Pct of Minority` > 100, 100, `Pct of Minority`)) %>% 
  mutate(`Pct of Minority` = paste0(`Pct of Minority`,"%")) %>% 
  arrange(desc(`Minority`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "People of Color in Least Walkable Block Groups", align = "r", col.names = c(names(.)[1:2],"Number of POC in Block Groups","Pct of POC in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapLowWalkLEH, fig.align = "center", fig.cap="Limited English speaking households per square kilometer in least walkable Census Block Groups (i.e., Walkability Index Score less than or equal to 5.8)."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(leastWalkable, unit = "mi") +
  tm_fill(col = "eng_limitE", alpha = 0.9, style = "fisher",
          colorNA = "white", colorNULL = "gray", title = "Limited English\nSpeaking Households\nPer Square Kilometer",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd", convert2density = TRUE, area = "SqKm") +
  tm_shape(vt_tracts_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "Limited English\nSpeaking\nHouseholds\nin Least Walkable\nBlock Groups", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_walkLEH.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_walkLEH.png")
```

```{r tabLowWalkLEH, fig.align = "center", fig.cap="Limited English speaking households by municipality in least walkable Census Block Groups (i.e., Walkability Index Score less than or equal to 5.8)."}
leastWalkable %>%
  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)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Limited English Households` > 0 & `Pct of Limited English Households` > 0) %>% 
  mutate(`Pct of Limited English Households` = if_else(`Pct of Limited English Households` > 100, 100, `Pct of Limited English Households`)) %>% 
  mutate(`Pct of Limited English Households` = paste0(`Pct of Limited English Households`,"%")) %>% 
  arrange(desc(`Limited English Households`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Limited English Speaking Households in Least Walkable Block Groups", align = "r", col.names = c(names(.)[1:2],"Number of Limited English Households in Block Groups","Pct of Limited English Households in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapLowWalkNoCarHH, fig.align = "center", fig.cap="Households without a car per square kilometer in least walkable Census Block Groups (i.e., Walkability Index Score less than or equal to 5.8)."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(leastWalkableTract, unit = "mi") +
  tm_fill(col = "HHnoCarE", alpha = 0.9, style = "fisher",
          colorNA = "white", colorNULL = "gray", title = "No Car Households\nPer Square Kilometer",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd", convert2density = TRUE, area = "SqKm") +
  tm_shape(vt_tracts_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "No Car Households\nin Least Walkable\nBlock Groups", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_walkNoCarHH.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_walkNoCarHH.png")
```

```{r tabLowWalkNoCarHH, fig.align = "center", fig.cap="Households without a car by municipality in least walkable Census Block Groups (i.e., Walkability Index Score less than or equal to 5.8)."}
leastWalkableTract %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `No Car Households` = sum(HHnoCarE,na.rm = TRUE),
            `Pct of No Car Households` = round(sum(HHnoCarE)/max(TotalNoCar)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`No Car Households` > 0 & `Pct of No Car Households` > 0) %>% 
  mutate(`Pct of No Car Households` = if_else(`Pct of No Car Households` > 100, 100, `Pct of No Car Households`)) %>% 
  mutate(`Pct of No Car Households` = paste0(`Pct of No Car Households`,"%")) %>% 
  arrange(desc(`No Car Households`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "No Car Households in Least Walkable Block Groups", align = "r", col.names = c(names(.)[1:2],"Number of No Car Households in Block Groups","Pct of No Car Households in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```


## Highest Transportation Cost Burden and Priority Populations {-}
```{r mapHiBurdenU18, fig.align = "center", fig.cap="Persons under 18 per square kilometer in Census Block Groups with the highest transportation cost burden (i.e., transportation cost burden in the 80th percentile for the state)."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(HiBurden, unit = "mi") +
  tm_fill(col = "under18E", alpha = 0.9, style = "fisher",
          colorNA = "white", colorNULL = "gray", title = "Persons Under 18\nPer Square Kilometer",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd", convert2density = TRUE, area = "SqKm") +
  tm_shape(vt_tracts_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "Persons Under 18\nin Highest\nTransportation\nCost Burden\nBlock Groups", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_HiBurdenU18.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_HiBurdenU18.png")
```

```{r tabHiBurdenU18, fig.align = "center", fig.cap="Persons under 18 by municipality in Census Block Groups with the highest transportation cost burden (i.e., transportation cost burden in the 80th percentile for the state)."}
HiBurden %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Under 18` = sum(under18E,na.rm = TRUE),
            `Pct of Under 18` = round(sum(under18E)/max(TotalUnder18)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Under 18` > 0 & `Pct of Under 18` > 0) %>%
  mutate(`Pct of Under 18` = if_else(`Pct of Under 18` > 100, 100, `Pct of Under 18`)) %>% 
  mutate(`Pct of Under 18` = paste0(`Pct of Under 18`,"%")) %>% 
  arrange(desc(`Under 18`)) %>%  
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Persons Under 18 in Highest Transportation Cost Burden Block Groups", align = "r", col.names = c(names(.)[1:2],"Number of Under 18 in Block Groups","Pct of Under 18 in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapHiBurdenDisabled, fig.align = "center", fig.cap="Disabled persons per square kilometer in Census Block Groups with the highest transportation cost burden (i.e., transportation cost burden in the 80th percentile for the state)."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(HiBurdenTract, unit = "mi") +
  tm_fill(col = "disabledOver18E", alpha = 0.9, style = "fisher",
          colorNA = "white", colorNULL = "gray", title = "Disabled Persons\nPer Square Kilometer",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd", convert2density = TRUE, area = "SqKm") +
  tm_shape(vt_tracts_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "Disabled Persons\nin Highest\nTransportation\nCost Burden\nBlock Groups", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_HiBurdenDisabled.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_HiBurdenDisabled.png")
```

```{r tabHiBurdenDisabled, fig.align = "center", fig.cap="Disabled persons by municipality in Census Tracts with the highest transportation cost burden (i.e., transportation cost burden in the 80th percentile for the state)."}
HiBurdenTract %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Tracts` = n(),
            `Disabled` = sum(disabledOver18E,na.rm = TRUE),
            `Pct of Disabled` = round(sum(disabledOver18E)/max(TotalDisabled)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Disabled` > 0 & `Pct of Disabled` > 0) %>%
  mutate(`Pct of Disabled` = if_else(`Pct of Disabled` > 100, 100, `Pct of Disabled`)) %>% 
  mutate(`Pct of Disabled` = paste0(`Pct of Disabled`,"%")) %>% 
  arrange(desc(`Disabled`)) %>%  
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Disabled Persons in Highest Transportation Cost Burden Census Tracts", align = "r", col.names = c(names(.)[1:2],"Number of Disabled in Census Tracts","Pct of Disabled in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapHiBurdenU5, fig.align = "center", fig.cap="Children under 5 per square kilometer in Census Block Groups with the highest transportation cost burden (i.e., transportation cost burden in the 80th percentile for the state)."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(HiBurden, unit = "mi") +
  tm_fill(col = "under5E", alpha = 0.9, style = "fisher",
          colorNA = "white", colorNULL = "gray", title = "Children Under 5\nPer Square Kilometer",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd", convert2density = TRUE, area = "SqKm") +
  tm_shape(vt_tracts_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "Children Under 5\nin Highest\nTransportation\nCost Burden\nBlock Groups", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_HiBurdenU5.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_HiBurdenU5.png")
```

```{r tabHiBurdenU5, fig.align = "center", fig.cap="Children under 5 by municipality in Census Block Groups with the highest transportation cost burden (i.e., transportation cost burden in the 80th percentile for the state)."}
HiBurden %>%
  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)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Under 5` > 0 & `Pct of Under 5` > 0) %>%
  mutate(`Pct of Under 5` = if_else(`Pct of Under 5` > 100, 100, `Pct of Under 5`)) %>% 
  mutate(`Pct of Under 5` = paste0(`Pct of Under 5`,"%")) %>% 
  arrange(desc(`Under 5`)) %>%  
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Children Under 5 in Highest Transportation Cost Burden Block Groups", align = "r", col.names = c(names(.)[1:2],"Number of Under 5 in Block Groups","Pct of Under 5 in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapHiBurdenINC, fig.align = "center", fig.cap="Low income persons per square kilometer in Census Block Groups with the highest transportation cost burden (i.e., transportation cost burden in the 80th percentile for the state)."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(HiBurden, unit = "mi") +
  tm_fill(col = "num2povE", alpha = 0.9, style = "fisher",
          colorNA = "white", colorNULL = "gray", title = "Low Income\nPersons\nPer Square Kilometer",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd", convert2density = TRUE, area = "SqKm") +
  tm_shape(vt_tracts_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "Low Income\nPersons\nin Highest\nTransportation\nCost Burden\nBlock Groups", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_HiBurdenINC.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_HiBurdenINC.png")
```

```{r tabHiBurdenINC, fig.align = "center", fig.cap="Low income persons by municipality in Census Block Groups with the highest transportation cost burden (i.e., transportation cost burden in the 80th percentile for the state)."}
HiBurden %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Low Income` = sum(num2povE,na.rm = TRUE),
            `Pct of Low Income` = round(sum(num2povE)/max(TotalLowInc)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Low Income` > 0 & `Pct of Low Income` > 0) %>%
  mutate(`Pct of Low Income` = if_else(`Pct of Low Income` > 100, 100, `Pct of Low Income`)) %>% 
  mutate(`Pct of Low Income` = paste0(`Pct of Low Income`,"%")) %>% 
  arrange(desc(`Low Income`)) %>%  
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Low Income Persons in Highest Transportation Cost Burden Block Groups", align = "r", col.names = c(names(.)[1:2],"Number of Low Income in Block Groups","Pct of Low Income in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapHiBurdenNoCarHH, fig.align = "center", fig.cap="Households without a car per square kilometer in Census Block Groups with the highest transportation cost burden (i.e., transportation cost burden in the 80th percentile for the state)."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(HiBurdenTract, unit = "mi") +
  tm_fill(col = "HHnoCarE", alpha = 0.9, style = "fisher",
          colorNA = "white", colorNULL = "gray", title = "No Car Households\nPer Square Kilometer",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd", convert2density = TRUE, area = "SqKm") +
  tm_shape(vt_tracts_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "No Car Households\nin Highest\nTransportation\nCost Burden\nBlock Groups", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_HiBurdenNoCarHH.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_HiBurdenNoCarHH.png")
```

```{r tabHiBurdenNoCarHH, fig.align = "center", fig.cap="Households without a car by municipality in Census Block Groups with the highest transportation cost burden (i.e., transportation cost burden in the 80th percentile for the state)."}
HiBurdenTract %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `No Car Households` = sum(HHnoCarE,na.rm = TRUE),
            `Pct of No Car Households` = round(sum(HHnoCarE)/max(TotalNoCar)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`No Car Households` > 0 & `Pct of No Car Households` > 0) %>%
  mutate(`Pct of No Car Households` = if_else(`Pct of No Car Households` > 100, 100, `Pct of No Car Households`)) %>% 
  mutate(`Pct of No Car Households` = paste0(`Pct of No Car Households`,"%")) %>% 
  arrange(desc(`No Car Households`)) %>%  
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "No Car Households in Highest Transportation Cost Burden Block Groups", align = "r", col.names = c(names(.)[1:2],"Number of No Car Households in Block Groups","Pct of No Car Households in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapHiBurdenMIN, fig.align = "center", fig.cap="People of Color per square kilometer in Census Block Groups with the highest transportation cost burden (i.e., transportation cost burden in the 80th percentile for the state)."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(HiBurden, unit = "mi") +
  tm_fill(col = "minorityE", alpha = 0.9, style = "fisher",
          colorNA = "white", colorNULL = "gray", title = "People of Color\nPer Square Kilometer",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd", convert2density = TRUE, area = "SqKm") +
  tm_shape(vt_tracts_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "People of Color\nin Highest\nTransportation\nCost Burden\nBlock Groups", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_HiBurdenMIN.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_HiBurdenMIN.png")
```

```{r tabHiBurdenMIN, fig.align = "center", fig.cap="People of Color by municipality in Census Block Groups with the highest transportation cost burden (i.e., transportation cost burden in the 80th percentile for the state)."}
HiBurden %>%
  st_centroid(.) %>% 
  st_intersection(vt_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Minority` = sum(minorityE,na.rm = TRUE),
            `Pct of Minority` = round(sum(minorityE)/max(TotalMin)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Minority` > 0 & `Pct of Minority` > 0) %>%
  mutate(`Pct of Minority` = if_else(`Pct of Minority` > 100, 100, `Pct of Minority`)) %>% 
  mutate(`Pct of Minority` = paste0(`Pct of Minority`,"%")) %>% 
  arrange(desc(`Minority`)) %>%  
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "People of Color in Highest Transportation Cost Burden Block Groups", align = "r", col.names = c(names(.)[1:2],"Number of POC in Block Groups","Pct of POC in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

```{r mapHiBurdenLEH, fig.align = "center", fig.cap="Limited English speaking households per square kilometer in Census Block Groups with the highest transportation cost burden (i.e., transportation cost burden in the 80th percentile for the state)."}
m <- tm_shape(vt_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(HiBurden, unit = "mi") +
  tm_fill(col = "eng_limitE", alpha = 0.9, style = "fisher",
          colorNA = "white", colorNULL = "gray", title = "Limited English\nHouseholds\nPer Square Kilometer",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd", convert2density = TRUE, area = "SqKm") +
  tm_shape(vt_tracts_sf) + tm_borders(alpha = 0.5, lwd = 0.2) +
  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.1) +
  tm_shape(I91roadSegment) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I91roadSegment2) +
  tm_symbols(shape = I91, border.lwd = NA, size = 0.1) +
  tm_shape(I93roadSegment) +
  tm_symbols(shape = I93, border.lwd = NA, size = 0.1) +
  tm_shape(vt_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "Limited English\nHouseholds\nin Highest\nTransportation\nCost Burden\nBlock Groups", 
            frame = FALSE, main.title.size = 0.8,
            legend.outside = TRUE,
            legend.title.size = 0.8,
            legend.outside.position = c("right", "top"))

tmap_save(m, "DATA/transport/VT/figures/vt_HiBurdenLEH.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/VT/figures/vt_HiBurdenLEH.png")
```

```{r tabHiBurdenLEH, fig.align = "center", fig.cap="Limited English speaking households by municipality in Census Block Groups with the highest transportation cost burden (i.e., transportation cost burden in the 80th percentile for the state)."}
HiBurden %>%
  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)) %>% 
  rename("City/Town" = NAME) %>% 
  filter(`Limited English Households` > 0 & `Pct of Limited English Households` > 0) %>%
  mutate(`Pct of Limited English Households` = if_else(`Pct of Limited English Households` > 100, 100, `Pct of Limited English Households`)) %>% 
  mutate(`Pct of Limited English Households` = paste0(`Pct of Limited English Households`,"%")) %>% 
  arrange(desc(`Limited English Households`)) %>%  
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Limited English Households in Highest Transportation Cost Burden Block Groups", align = "r", col.names = c(names(.)[1:2],"Number of Limited English Households in Block Groups","Pct of Limited English Households in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```