Transport_RhodeIslandPDF.Rmd

---
title: "Transportation Options in Rhode Island"
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 Rhode Island

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

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

rm(ne_blkgrp_sf90)

### Read in RI transit data. See http://www.rigis.org/
# read in RIPTA bus routes
ri_busroutes_sf <- st_read(dsn = "DATA/transport/RI/RIPTA_Bus_Routes",
                           layer = "RIPTA_Bus_Routes")
# read in RIPTA bus stops
ri_busstops_sf <- st_read(dsn = "DATA/transport/RI/RIPTA_Bus_Stops",
                          layer = "RIPTA_Bus_Stops")
# read in RI ferry routes
ri_ferryroutes_sf <- st_read(dsn = "DATA/transport/RI/RI_Ferry_Routes",
                             layer = "Ferry_Routes")
# read in RI railroad rights of way
ri_railrow_sf <- st_read(dsn = "DATA/transport/RI/RI_Railroad_Rights_of_Way",
                         layer = "Railroad_Rights_of_Way") %>% 
  filter(STATUS == "Active" & RAILUSE == "Passenger")
# read in RI park and ride stops
ri_parknride_sf <- st_read(dsn = "DATA/transport/RI/RIPTA_Park_and_Ride_Stops",
                           layer = "RIPTA_Park_and_Ride_Stops")

# Convert to projected local CRS EPSG:2840: NAD83(HARN) / Rhode Island
ri_busstops_sf <- st_transform(ri_busstops_sf, crs = 2840) 

ri_blkgrps_sf <- ne_blkgrp_sf %>% 
  filter(STATE == "Rhode Island") %>% 
  st_transform(., crs = 2840) %>% 
  mutate(PopAcre = totalpopE/(bg_area_m2*0.000247105))

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

# Get rid of empty geometries
empty_geo <- st_is_empty(ri_blkgrps_sf)
ri_blkgrps_sf <- ri_blkgrps_sf[!empty_geo,]
empty_geo <- st_is_empty(ri_tracts_sf)
ri_tracts_sf <- ri_tracts_sf[!empty_geo,]

# 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 ma_blkgrps and ma_tracts
# st_layers(dsn = "DATA/FEMA/MA")
ri_blkgrps_developed <- st_read(dsn = "DATA/FEMA/RI",
                                layer = "ri_blkgrps_developed") 

ri_blkgrps_developed <- ri_blkgrps_sf %>% 
  mutate(RI_LOWINC = if_else(RI_INCOME == "I",totalpopE,0)) %>% 
  mutate(RI_LOWINC = replace_na(RI_LOWINC,0)) %>% 
  mutate(RI_MINORITIES = if_else(RI_MINORITY == "M",totalpopE,0)) %>% 
  mutate(RI_MINORITIES = replace_na(RI_MINORITIES,0)) %>% 
  as.data.frame() %>% 
  dplyr::select(-geometry) %>%
  left_join(ri_blkgrps_developed, ., by = "GEOID") %>% 
  st_transform(., crs = 2840) %>% 
  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,
         NewRI_LOWINC = RI_LOWINC*Proportion,
         NewRI_MINORITIES = RI_MINORITIES*Proportion,
         NewPopAcre = NewPop/(NewArea*0.000247105))

ri_tracts_developed <- st_read(dsn = "DATA/FEMA/RI",
                                layer = "ri_tract_developed") %>% 
  left_join(., as.data.frame(ri_tracts_sf), by = "GEOID") %>% 
  st_transform(., crs = 2840) %>% 
  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
ri_blkgrps_developed_simple <- st_simplify(ri_blkgrps_developed,
                                           dTolerance = 100) %>% 
  filter(!st_is_empty(.)) %>% 
  st_make_valid()
# repeat for tracts
ri_tracts_developed_simple <- st_simplify(ri_tracts_developed, 
                                          dTolerance = 100) %>% 
  filter(!st_is_empty(.)) %>% 
  st_make_valid()

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

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

# Use areal interpolation to calculate priority populations within buffer of accessibility
ri_busBuff400m_sf <- ri_blkgrps_developed %>% 
  select(GEOID,NewArea:NewRI_MINORITIES) %>% 
  mutate(OldArea = NewArea) %>% 
  st_intersection(ri_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),
         NewRI_LOWINC = as.integer(NewRI_LOWINC*Proportion),
         NewRI_MINORITIES = as.integer(NewRI_MINORITIES*Proportion))

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

# # Use areal interpolation to calculate priority populations within buffer of accessibility
# ri_busBuff400m_sf <- ri_blkgrps_sf %>% 
#   mutate(RI_LOWINC = if_else(RI_INCOME == "I",totalpopE,0)) %>% 
#   mutate(RI_LOWINC = replace_na(RI_LOWINC,0)) %>% 
#   mutate(RI_MINORITIES = if_else(RI_MINORITY == "M",totalpopE,0)) %>% 
#   mutate(RI_MINORITIES = replace_na(RI_MINORITIES,0)) %>% 
#   select(GEOID,
#          totalpopE,
#          minorityE,
#          under5E,
#          over64E,
#          under18E,
#          eng_limitE,
#          num2povE,
#          lthsE,
#          RI_LOWINC,
#          RI_MINORITIES) %>% 
#   mutate(OldArea = st_area(.)) %>% 
#   st_intersection(ri_busBuff400m,.) %>% 
#   mutate(NewArea = st_area(.),
#          Proportion = NewArea/OldArea,
#          NewPop = as.integer(totalpopE*Proportion),
#          NewMinority = as.integer(minorityE*Proportion),
#          NewUnder5 = as.integer(under5E*Proportion),
#          NewOver64 = as.integer(over64E*Proportion),
#          NewUnder18 = as.integer(under18E*Proportion),
#          NewEng_limit = as.integer(eng_limitE*Proportion),
#          NewPov = as.integer(num2povE*Proportion),
#          NewLths = as.integer(lthsE*Proportion),
#          NewRI_LOWINC = as.integer(RI_LOWINC*Proportion),
#          NewRI_MINORITIES = as.integer(RI_MINORITIES*Proportion))
# 
# # Repeat for tracts
# ri_busBuff400mTracts_sf <- ri_tracts_sf %>% 
#   select(GEOID,
#          disabledOver18E,
#          HHnoCarE) %>% 
#   mutate(OldArea = st_area(.)) %>% 
#   st_intersection(ri_busBuff400m,.) %>% 
#   mutate(NewArea = st_area(.),
#          Proportion = NewArea/OldArea,
#          NewDisabled = as.integer(disabledOver18E*Proportion),
#          NewNoCar = as.integer(HHnoCarE*Proportion))

# Compute total block group populations within bus stop buffer
ri_busBuff400m_df <- ri_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),
            `RI Low Income` = sum(NewRI_LOWINC),
            `RI Minority` = sum(NewRI_MINORITIES)) %>% 
  gather(key = Group, value = BusPop)

# Compute total tract populations within bus stop buffer
ri_busBuff400mTracts_df <- ri_busBuff400mTracts_sf %>% 
  as.data.frame() %>% 
  summarize(`Disabled` = sum(NewDisabled),
            `No Car HH` = sum(NewNoCar)) %>% 
  gather(key = Group, value = BusPop)
# Compute total tract populations within the state for same groups
ri_tract_pops_df <- ri_tracts_sf %>% 
  as.data.frame() %>% 
  summarize(`Disabled` = sum(disabledOver18E),
            `No Car HH` = sum(HHnoCarE)) %>% 
  gather(key = Group, value = RIPop) %>% 
  left_join(.,ri_busBuff400mTracts_df, by = "Group")
  

# Compute populations for state, join with buffer pops
ri_busAccessPops_df <- ri_blkgrps_sf %>% 
  as.data.frame() %>% 
  mutate(RI_LOWINC = if_else(RI_INCOME == "I",totalpopE,0)) %>% 
  mutate(RI_LOWINC = replace_na(RI_LOWINC,0)) %>% 
  mutate(RI_MINORITIES = if_else(RI_MINORITY == "M",totalpopE,0)) %>% 
  mutate(RI_MINORITIES = replace_na(RI_MINORITIES,0)) %>% 
  summarize(`Total Pop` = sum(totalpopE),
            Minority = sum(minorityE),
            `Under 5` = sum(under5E),
            `Over 64` = sum(over64E),
            `Under 18` = sum(under18E),
            `Limited English HH` = sum(eng_limitE),
            `Low Income` = sum(num2povE),
            `No HS Dip` = sum(lthsE),
            `RI Low Income` = sum(RI_LOWINC, na.rm = TRUE),
            `RI Minority` = sum(RI_MINORITIES, na.rm = TRUE)) %>% 
  gather(key = Group, value = RIPop) %>% 
  left_join(., ri_busBuff400m_df, by = "Group") %>% 
  rbind(.,ri_tract_pops_df) %>% 
  mutate(PctBus = BusPop/RIPop*100)

# create point layer of towns for context
ri_towns_sf_pts <- county_subdivisions(state = "RI", cb = TRUE) %>% 
  filter(NAME %in% c("Providence",
                     "Woonsocket",
                     "Pawtucket",
                     "Warwick",
                     "Bristol",
                     "Portsmouth",
                     "Newport",
                     "Narragansett",
                     "North Kingstown",
                     "Charlestown",
                     "Situate",
                     "Glocester",
                     "Burrillville",
                     "Scituate",
                     "Coventry",
                     "Hopkinton",
                     "South Kingstown",
                     "Tiverton",
                     "New Shoreham",
                     "Westerly")) %>% 
  st_transform(., crs = 2840) %>% 
  st_centroid(of_largest_polygon = TRUE)

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

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

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

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

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

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

# recode routes to match RIPTA categories
ri_busroutes_sf <- ri_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 <- ri_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 <- ri_blkgrps_sf %>% 
  as.data.frame() %>% 
  mutate(RI_LOWINC = if_else(RI_INCOME == "I", totalpopE, 0)) %>% 
  mutate(RI_LOWINC = replace_na(RI_LOWINC,0)) %>% 
  mutate(RI_MINORITY = if_else(RI_MINORITY == "M", totalpopE,0)) %>%
  mutate(RI_MINORITY = replace_na(RI_MINORITY,0)) %>% 
  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),
            TotalRI_LOWINC = sum(RI_LOWINC, na.rm = TRUE),
            TotalRI_MINORITY = sum(RI_MINORITY, na.rm = TRUE)) %>% 
  left_join(., statepopsdf_tracts,by="STATE")

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

ri_towns_sf <- county_subdivisions(state = "RI", cb = TRUE) %>% 
  st_transform(., crs = 2840)

townpopsdf_tracts <- ri_tracts_sf %>% 
  st_centroid(.) %>% 
  st_intersection(ri_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 <- ri_blkgrps_sf %>% 
  dplyr::select(-GEOID) %>% 
  st_centroid(.) %>% 
  st_intersection(ri_towns_sf,.) %>% 
  as.data.frame() %>% 
  mutate(RI_LOWINC = if_else(RI_INCOME == "I", totalpopE, 0)) %>% 
  mutate(RI_LOWINC = replace_na(RI_LOWINC,0)) %>% 
  mutate(RI_MINORITY = if_else(RI_MINORITY == "M", totalpopE,0)) %>%
  mutate(RI_MINORITY = replace_na(RI_MINORITY,0)) %>% 
  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),
            TotalRI_LOWINC = sum(RI_LOWINC, na.rm = TRUE),
            TotalRI_MINORITY = sum(RI_MINORITY, na.rm = TRUE)) %>% 
  left_join(.,townpopsdf_tracts,by="NAME")

ri_state_sf_cb <- ne_states_sf_cb %>% 
  filter(NAME == "Rhode Island") %>% 
  st_transform(., crs = st_crs(ri_blkgrps_sf))

# Identify block groups >= 1 mile from transit stops
# Create a 1 mile buffer around bus stops
ri_busBuff1mile <- st_buffer(ri_busstops_sf,dist = 1609.34) %>% 
  st_union() %>% 
  st_as_sf()

# Identify block groups that do not intersect with buffer
# ri_blkgrps_sf_noBus <- ri_blkgrps_sf %>% 
#   filter(st_disjoint(.,ri_busBuff1mile,sparse = FALSE))

ri_blkgrps_sf_noBus <- ri_blkgrps_developed %>% 
  filter(st_disjoint(.,ri_busBuff1mile,sparse = FALSE))

# Identify tracts that do not intersect with buffer
# ri_tracts_sf_noBus <- ri_tracts_sf %>% 
#   filter(st_disjoint(.,ri_busBuff1mile,sparse = FALSE))

ri_tracts_sf_noBus <- ri_tracts_developed %>% 
  filter(st_disjoint(.,ri_busBuff1mile,sparse = FALSE))

# write out to file for later analysis
save(ri_blkgrps_sf_noBus, ri_tracts_sf_noBus, 
     ri_busAccessPops_df, ri_blkgrps_sf_noBus, ri_busBuff400m_sf, 
     file = "DATA/transport/RI/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. 

Rhode Island is served primarily by the Rhode Island Public Transit Authority (RIPTA) which is a quasi-public, independent authority. This public transit service system is made up almost entirely by fixed-route bus service, which passes through 36 of Rhode Island's 39 municipalities. Amtrak train service passes through Rhode Island as well, but it is not considered in this analysis because it is not managed by the state of Rhode Island. This analysis looks at bus service access and frequency in Rhode Island. Figure \@ref(fig:transitMap) shows `r nrow(ri_busroutes_sf)` fixed public transit routes across the state along with population density.

```{r transitMap, fig.align = "center", fig.cap="Transit service and population density in Rhode Island."}
# bbox_new <- st_bbox(ri_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

# # create simplified version to speed up mapping
# ri_blkgrps_developed_simple <- st_simplify(ri_blkgrps_developed, dTolerance = 100)

m <- tm_shape(ri_towns_sf, unit = "mi") + 
  tm_borders(col = "gray", lwd = 0.2) +
  tm_shape(ri_blkgrps_developed_simple) +
  tm_fill(col = "NewPopAcre", alpha = 0.8, 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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(ri_busroutes_sf) + 
  tm_lines(col = "type", palette = c("purple","red","blue"), 
           colorNA = NULL, title.col = "Route Type", lwd = 2, 
           alpha = 0.5) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "Rhode Island\nPublic Transit\nAuthority\n(RIPTA)\nBus 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/RI/figures/ri_transit.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_transit.png")
```

Approximately `r paste0(round(ri_busAccessPops_df %>% filter(Group == "Total Pop") %>% dplyr::select(PctBus) %>% pull(),1),"%")` of the state population lives within reasonable walking distance of a bus stop (defined as approximately 400 meters or 1/4 mile). 

However, this access varies by population group. Figure \@ref(fig:loliBusDist) below compares the percentages of different population groups living within reasonable walking distance of bus stops across the state. Most priority populations, including transit-dependent groups, have access similar to or better than the general population. However, note that people over age 64 have less access 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 across Rhode Island."}
ri_busAccessPops_df %>% 
  mutate(Group = recode(Group, "Minority" = "People of Color",
                        "RI Minority" = "RI POC",
                        "No HS Dip" = "No HS Diploma")) %>%
  ggplot(aes(x = reorder(Group,-PctBus), 
             y = PctBus)) +
  geom_segment(aes(x = reorder(Group,-PctBus), xend = reorder(Group,-PctBus),
                   y = ri_busAccessPops_df[1,4], yend = PctBus), 
               color = "skyblue") +
  geom_point(color = "blue", size = 4, alpha = 0.8) +
  coord_flip() + xlab("") + ylab("") + ggtitle("Rhode Island Populations 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 = PctBus + 0.2 * sign(PctBus), 
                label = paste0(round(PctBus,0),"%")), 
            hjust = -0.1, vjust = -.8, size = 3,
            color=rgb(100,100,100, maxColorValue=255)) +
  scale_y_continuous(labels = function(x) paste0(x, "%")) +
  geom_hline(yintercept = ri_busAccessPops_df[1,4], linetype = "dashed") +
  geom_text(aes(x = "Under 18", y = 69, label = "Above state\naverage"),
            color = "gray48") +
  geom_text(aes(x = "Under 18", y = 45, label = "Below state\naverage"),
            color = "gray48") +
  expand_limits(y = c(40,96))

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

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

```{r statsBusDist, echo=FALSE, message=FALSE, warning=FALSE, fig.align = "center", fig.cap="Rhode Island populations living within 1/4 mile (400 meters) of a bus stop. Based on ACS 2018 Block Group data."}
# See table of distance by group
ri_busAccessPops_df %>% 
  mutate(Group = recode(Group, "Minority" = "People of Color",
                        "RI Minority" = "RI POC",
                        "No HS Dip" = "No HS Diploma")) %>%
  mutate(`% Bus Access` = paste0(round(PctBus,1),"%")) %>% 
  select(-PctBus) %>% 
  kableExtra::kable(longtable = T, booktabs = T, 
                    format.args = list(big.mark = ','), align = "r",
                    caption = "Populations Living within 400m of Bus Stops", col.names = c("Group","RI Pop","Pop with Bus Access","Pct Pop with Bus Access")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header"))
  
```

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 <- ri_blkgrps_sf_noBus %>%
  st_centroid(.) %>% 
  st_intersection(ri_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) %>% 
  mutate(`Pct of Over 64` = paste0(`Pct of Over 64`,"%")) %>% 
  arrange(desc(`Over 64`))

TransitNoAccess2 <- ri_tracts_sf_noBus %>%
  st_centroid(.) %>% 
  st_intersection(ri_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) %>% 
  mutate(`Pct of Disabled` = paste0(`Pct of Disabled`,"%")) %>% 
  arrange(desc(Disabled))
```
Figures \@ref(fig:mapTransitAccessO64) and \@ref(fig:mapTransitAccessDisabled) below highlight developed portions of Census Block Groups and Tracts across Rhode Island which are more than 1 mile from a transit stop and the numbers of persons over 64 or disabled who live there. Tables \@ref(tab:tabTransitAccessO64) and \@ref(tab:tabTransitAccessDisabled) show breakdowns by municipality. For example, Table \@ref(tab:tabTransitAccessO64) shows that `r formatC(as.numeric(TransitNoAccess1[1,3]),big.mark=",")` people over 64 across `r TransitNoAccess1[1,2]` Block Groups in `r TransitNoAccess1[1,1]` resided one or more miles from the nearest transit stop. These people represented `r TransitNoAccess1[1,4]` of people over 64 in `r TransitNoAccess1[1,1]`. Similarly, Table \@ref(tab:tabTransitAccessDisabled) shows that `r formatC(as.numeric(TransitNoAccess2[1,3]),big.mark=",")` people with disabilities across `r TransitNoAccess1[1,2]` Census Tracts in `r TransitNoAccess2[1,1]` resided one or more miles from the nearest transit stop. These people represented `r TransitNoAccess2[1,4]` of people with disabilities in `r TransitNoAccess2[1,1]`.
Maps and tables of other priority populations without access to transit can be found in Appendix B.

```{r mapTransitAccessO64, fig.align = "center", fig.cap="Persons over 64 who are one or more miles from the nearest transit stop by Census Block Group."}
m <- tm_shape(ri_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(ri_blkgrps_sf_noBus, unit = "mi") +
  tm_fill(col = "NewOver64", alpha = 0.9, style = "quantile",
          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(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  # tm_shape(ri_busroutes_sf) +
  # tm_lines(col = "type", palette = c("purple","red","blue"),
  #          colorNA = NULL, title.col = "Route Type", lwd = 2,
  #          alpha = 0.5) +
  tm_shape(ri_busroutes_sf) +
  tm_lines(lwd = 2, alpha = 0.5) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_add_legend(type = "line", lwd=2, alpha = 0.5, labels = "Transit Route") +
  tm_layout(title = "Persons Over 64\n1+ Miles from\nNearest Transit\nStop", 
            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/RI/figures/ri_transitO64access.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_transitO64access.png")
```

```{r tabTransitAccessO64, fig.align = "center", fig.cap="Persons over 64 in Census Block Groups one or more miles from the nearest bus stop."}
ri_blkgrps_sf_noBus %>%
  st_centroid(.) %>% 
  st_intersection(ri_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) %>% 
  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 One or More Miles from Nearest Transit Stop", 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 mapTransitAccessDisabled, fig.align = "center", fig.cap="Disabled persons who are one or more miles from the nearest transit stop by Census Tract."}
m <- tm_shape(ri_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(ri_tracts_sf_noBus, unit = "mi") +
  tm_fill(col = "NewDisabled", alpha = 0.9, style = "quantile",
          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(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  # tm_shape(ri_busroutes_sf) +
  # tm_lines(col = "type", palette = c("purple","red","blue"),
  #          colorNA = NULL, title.col = "Route Type", lwd = 2,
  #          alpha = 0.5) +
  tm_shape(ri_busroutes_sf) +
  tm_lines(lwd = 2, alpha = 0.5) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_add_legend(type = "line", lwd=2, alpha = 0.5, labels = "Transit Route") +
  tm_layout(title = "Disabled Persons\n1+ Miles from\nNearest Transit\nStop", 
            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/RI/figures/ri_transitDisabledaccess.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_transitDisabledaccess.png")
```

```{r tabTransitAccessDisabled, fig.align = "center", fig.cap="Persons over 64 in Census Block Groups one or more miles from the nearest bus stop."}
ri_tracts_sf_noBus %>%
  st_centroid(.) %>% 
  st_intersection(ri_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) %>% 
  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 One or More Miles from Nearest Transit Stop", 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")) 
```


## 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 is scheduled or actually 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:mapGTFS) below shows average scheduled headways for bus routes throughout Rhode Island.[^headwayReal] Table \@ref(tab:tabGTFS) shows a breakdown of headway times by service type. Headways vary significantly by route and by type of service across the state. 

```{r GTFS, include=FALSE}
# Read in RI GTFS data
# # Note that Official Download URL is erratic in availability
# ri_transit <- read_gtfs("https://www.ripta.com/stuff/contentmgr/files/0/3cda81dfa140edbe9aae214b26245b4a/files/google_transit.zip") %>% 
#   gtfs_as_sf()
# Use previously downloaded GTFS for June 2019
ri_transit <- read_gtfs("DATA/transport/RI/RI_gtfs_files/RIPTA062019gtfs.zip") %>% 
  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.
ri_transit <- set_date_service_table(ri_transit)

# 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.
holidays = tribble(~date, ~holiday,
                   ymd("2019-07-04"), "Independence Day",
                   ymd("2019-09-02"), "Labor Day",
                   ymd("2019-12-25"), "Christmas")

calendar = tibble(date = unique(ri_transit$.$date_service_table$date)) %>% 
  mutate(
    weekday = (function(date) {
      c("Sunday", "Monday", "Tuesday", 
        "Wednesday", "Thursday", "Friday", 
        "Saturday")[as.POSIXlt(date)$wday + 1]
    })(date)
  )

calendar <- calendar %>% left_join(holidays, by = "date")

# 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. 
ri_transit <- set_servicepattern(ri_transit)
# 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.
# date_servicepattern_table <- ri_transit$.$date_servicepattern_table %>%
#   left_join(calendar, by = "date")
# 
# ggplot(date_servicepattern_table) + 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")

# Use service pattern id for daily M-F service selected based on graphic. CHECK GGPLOT ABOVE TO GET CORRECT SERVICE PATTERN ID!
service_ids <- ri_transit$.$service_pattern %>% 
  filter(servicepattern_id == "s_fb9e6e7") %>% 
  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.
daily_stop_freq <- get_stop_frequency(ri_transit, start_hour = 6, end_hour = 21,
                                      service_ids = service_ids)

# Convert stops to points for mapping
ri_transit_stops_sf <- stops_as_sf(ri_transit$stops)

# Join headway frequencies to stops
ri_transit_stops_sf <- ri_transit_stops_sf %>% 
  inner_join(daily_stop_freq, by = "stop_id")

# map it out
# tmap_mode("view")
# tm_shape(ri_transit_stops_sf) + tm_dots(col = "headway", alpha = 0.6)

# use the get_route_frequency function to summarise transit service by route, for the same time period.
daily_route_freq <- get_route_frequency(ri_transit, service_ids = service_ids,
                                        start_hour = 6, end_hour = 21)

# Join the route frequencies to geometry for mapping
# get_route_geometry needs a gtfs object that includes shapes as simple feature data frames
ri_routes_sf <- get_route_geometry(ri_transit, service_ids = service_ids)

# join calculated frequencies to geometry
ri_routes_sf <- ri_routes_sf %>% 
  inner_join(daily_route_freq, by = "route_id")

# Calculate average headway for block groups within walking distance of bus stops and then compute weighted average headway by population group
ri_transit_stopHeadway_df <- ri_transit_stops_sf %>% 
  st_transform(., crs = 2840) %>% 
  st_join(., ri_busBuff400m_sf) %>% 
  as.data.frame() %>% 
  group_by(GEOID) %>% 
  summarize(AvgStopHeadway = mean(headway,na.rm = TRUE))
# Create same for tracts
ri_transit_stopHeadwayTracts_df <- ri_transit_stops_sf %>% 
  st_transform(., crs = 2840) %>% 
  st_join(., ri_busBuff400mTracts_sf) %>% 
  as.data.frame() %>% 
  group_by(GEOID) %>% 
  summarize(AvgStopHeadway = mean(headway,na.rm = TRUE))
# Create version of ri_transit_stopHeadwayTracts_df for rbind
ri_transit_stopHeadwayTracts <- ri_transit_stopHeadwayTracts_df %>% 
  as.data.frame() %>% 
  left_join(ri_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
ri_busHeadwayPops_df <- ri_busBuff400m_sf %>% 
  as.data.frame() %>% 
  left_join(.,ri_transit_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,
            `RI Low Income` = NewRI_LOWINC,
            `RI Minority` = NewRI_MINORITIES,
            AvgStopHeadway) %>% 
  gather(key = Group, value = Pop, `Total Pop`:`RI Minority`) %>% 
  group_by(Group) %>% 
  summarize(AvgWHeadway = weighted.mean(x = AvgStopHeadway, 
                                        w = Pop, na.rm = TRUE)) %>% 
  rbind(.,ri_transit_stopHeadwayTracts)
```


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

m <- tm_shape(ri_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(ri_routes_sf) + 
  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(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "Mean\nHeadways\nby Transit\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/RI/figures/ri_transitHeadways.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_transitHeadways.png")
```

```{r tabGTFS, fig.align = "center", fig.cap="Headway by service type for Rhode Island transit based on on GTFS schedule for Monday to Friday service, 6am - 9pm."}
# Show table of mean headways of service type
ri_busroutes_sf %>% 
  as.data.frame() %>% 
  mutate(ROUTE = as.character(ROUTE)) %>% 
  left_join(., as.data.frame(ri_routes_sf), by=c("ROUTE" = "route_id")) %>% 
  group_by(type) %>% 
  summarize(Number = n(), Min = min(mean_headways), Mean = mean(mean_headways), Max = max(mean_headways)) %>% 
  arrange(type) %>% 
  kableExtra::kable(longtable = T, booktabs = T, digits = 1, 
                    col.names = c("Service Type","Number of Routes",
                                  "Min",
                                  "Mean",
                                  "Max"), 
                    caption = "Headway by Service Type", 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"))
```

The average population-weighted headway for all transit in the state is `r round(ri_busHeadwayPops_df %>% filter(Group == "Total Pop") %>% select(AvgWHeadway) %>% pull(),1)` minutes. However headways, or service frequencies, vary significantly across population groups. Figure \@ref(fig:loliHeadway) below compares the population-weighted average headways for different groups across the state. Most groups have average headways shorter than the general population. However, persons over age 64 have headways that exceed that of the general population (`r round(ri_busHeadwayPops_df %>% filter(Group == "Over 64") %>% select(AvgWHeadway) %>% pull(),1)` minutes). The shortest population-weighted average headway, experienced by People of Color, as defined by Rhode Island environmental justice policy, is `r round(ri_busHeadwayPops_df %>% arrange(AvgWHeadway) %>% slice(.,1) %>% select(AvgWHeadway) %>% pull(),0)` minutes.

```{r loliHeadway, echo=FALSE, message=FALSE, warning=FALSE, fig.align = "center", fig.cap="Population-weighted average headways for groups within 400 meters of a bus stop across Rhode Island."}
ri_busHeadwayPops_df %>% 
  mutate(Group = recode(Group, "Minority" = "People of Color",
                        "RI Minority" = "RI POC",
                        "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(ri_busHeadwayPops_df[8,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\nRhode Island 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(ri_busHeadwayPops_df[8,2]), 
             linetype = "dashed") +
  geom_text(aes(x = "Under 18", y = 108, label = "Above state\naverage"),
            color = "gray48") +
  geom_text(aes(x = "Under 18", y = 80, label = "Below state\naverage"),
            color = "gray48") +
  expand_limits(y = c(60,110))

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

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

headway80thtracts <- ri_busBuff400mTracts_sf %>% 
  left_join(., ri_transit_stopHeadwayTracts_df, by = "GEOID") %>%
  filter(AvgStopHeadway >= quantile(ri_routes_sf$mean_headways,0.8))

headway80_1 <- headway80thblkgrps %>%
  st_centroid(.) %>% 
  st_intersection(ri_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) %>% 
  mutate(`Pct of Over 64` = paste0(`Pct of Over 64`,"%")) %>% 
  arrange(desc(`Over 64`))

headway80_2 <- headway80thtracts %>%
  st_centroid(.) %>% 
  st_intersection(ri_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) %>% 
  mutate(`Pct of Disabled` = paste0(`Pct of Disabled`,"%")) %>% 
  arrange(desc(Disabled))

# write out to file for use in later analysis
save(headway80thblkgrps, headway80thtracts,
     ri_busHeadwayPops_df, ri_transit_stopHeadway_df,
     file = "DATA/transport/RI/headway80th.Rds")

st_write(headway80thblkgrps,"DATA/transport/RI/headway80thblkgrps_bus.shp",
         delete_layer = TRUE)
st_write(headway80thtracts, "DATA/transport/RI/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:mapHeadwayOver64) and \@ref(fig:mapHeadwayDisabled) below highlight communities of persons over 64 or disabled persons within walking distance of transit 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 (`r quantile(ri_routes_sf$mean_headways,0.8)` minutes). Tables \@ref(tab:tabHeadwayOver64) and \@ref(tab:tabHeadwayDisabled) show breakdowns by municipality. For example, Table \@ref(tab:tabHeadwayOver64) shows that `r formatC(as.numeric(headway80_1[1,3]),big.mark=",")` people over 64 across `r headway80_1[1,2]` Block Groups in `r headway80_1[1,1]` resided in areas where the average headway of accessible transit stops exceeded `r quantile(ri_routes_sf$mean_headways,0.8)` minutes. These people represented `r headway80_1[1,4]` of people over 64 in `r headway80_1[1,1]`. Similarly, Table \@ref(tab:tabHeadwayDisabled) shows that `r formatC(as.numeric(headway80_2[1,3]),big.mark=",")` people with disabilities across `r headway80_2[1,2]` Census Tracts in `r headway80_2[1,1]` resided in areas where the average headway of accessible transit stops exceeded `r quantile(ri_routes_sf$mean_headways,0.8)` minutes. These people represented `r headway80_2[1,4]` of people with disabilities in `r headway80_2[1,1]`.
Maps and tables of other priority populations with excessive headways can be found in Appendix B. 

```{r mapHeadwayOver64, 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(ri_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(headway80thblkgrps) +
  tm_fill(col = "NewOver64", alpha = 0.9, style = "quantile",
          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(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  # tm_shape(ri_busroutes_sf) +
  # tm_lines(col = "type", palette = c("purple","red","blue"),
  #          colorNA = NULL, title.col = "Route Type", lwd = 2,
  #          alpha = 0.5) +
  tm_shape(ri_busroutes_sf) + tm_lines(lwd = 1, alpha = 0.3) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_add_legend(type = "line", lwd=1, alpha = 0.3, labels = "Transit Route") +
  tm_layout(title = paste("Persons Over 64","with Headways",paste0("Over ",round(quantile(ri_routes_sf$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/RI/figures/ri_headwayOver64.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_headwayOver64.png")
```

```{r tabHeadwayOver64, 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 %>%
  st_centroid(.) %>% 
  st_intersection(ri_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) %>% 
  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 Headways",paste0("Over ",quantile(ri_routes_sf$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 mapHeadwayDisabled, 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(ri_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(headway80thtracts) +
  tm_fill(col = "NewDisabled", alpha = 0.9, style = "quantile",
          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(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  # tm_shape(ri_busroutes_sf) +
  # tm_lines(col = "type", palette = c("purple","red","blue"),
  #          colorNA = NULL, title.col = "Route Type", lwd = 2,
  #          alpha = 0.5) +
  tm_shape(ri_busroutes_sf) + tm_lines(lwd = 1, alpha = 0.3) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_add_legend(type = "line", lwd=1, alpha = 0.3, labels = "Transit Route") +
  tm_layout(title = paste("Disabled Persons","with Headways",paste0("Over ",round(quantile(ri_routes_sf$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/RI/figures/ri_headwayDisabled.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_headwayDisabled.png")
```

```{r tabHeadwayDisabled, 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 %>%
  st_centroid(.) %>% 
  st_intersection(ri_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) %>% 
  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 Headways",paste0("Over ",quantile(ri_routes_sf$mean_headways,0.8)," 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"))

```


## 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 RI. 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
ri_blkgrp_walkability_sf <- walkability %>% 
  as.data.frame() %>% 
  select(GEOID10,SFIPS,CFIPS,TRFIPS,NatWalkInd) %>% 
  left_join(ri_blkgrps_sf,., by = c("GEOID" = "GEOID10"))

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

# Create summarized df of population weighted avg walkability by tract pop
ri_tract_walkability_df <- ri_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
ri_wAvgWalkability_df <- ri_blkgrp_walkability_sf %>% 
  as.data.frame() %>% 
  mutate(RI_LOWINC = if_else(RI_INCOME == "I",totalpopE,0)) %>% 
  mutate(RI_LOWINC = replace_na(RI_LOWINC,0)) %>% 
  mutate(RI_MINORITIES = if_else(RI_MINORITY == "M",totalpopE,0)) %>% 
  mutate(RI_MINORITIES = replace_na(RI_MINORITIES,0)) %>%
  transmute(`Total Pop` = totalpopE,
            Minority = minorityE,
            `Low Income` = num2povE,
            `Under 5` = under5E,
            `Under 18` = under18E,
            `Over 64` = over64E,
            `Limited English HH` = eng_limitE,
            `No HS Dip` = lthsE,
            `RI Low Income` = RI_LOWINC,
            `RI Minority` = RI_MINORITIES,
            NatWalkInd) %>% 
  gather(key = Group, value = Pop, `Total Pop`:`RI Minority`) %>% 
  group_by(Group) %>% 
  summarize(wNatWalkIndex = weighted.mean(NatWalkInd, Pop, na.rm = TRUE)) %>% 
  rbind(., ri_tract_walkability_df)

# save data for later analysis
save(ri_blkgrp_walkability_sf,
     ri_tract_walkability_sf, ri_wAvgWalkability_df,
     file = "DATA/transport/RI/walkability.Rds")
```

Figure \@ref(fig:mapWalkability) shows walkability scores by Census Block Group across Rhode Island.  The average population-weighted walkability score for the state is `r round(ri_wAvgWalkability_df %>% filter(Group == "Total Pop") %>% select(wNatWalkIndex) %>% pull(),1)`. However, walkability varies across the state. It is highest in and around the state's major urban centers, such as Providence, Portsmouth, and Newport, and lowest in the western suburbs.

```{r mapWalkability, echo=FALSE, message=FALSE, warning=FALSE, fig.align = "center", fig.cap="Walkability Index scores for Census Block Groups across Rhode Island."}
m <- tm_shape(ri_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(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_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/RI/figures/ri_walk.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_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 walkability index scores that are 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 Rhode Island."}
# Create lollipop graph of pop-weighted average walk index
ri_wAvgWalkability_df %>% 
  mutate(Group = recode(Group, "Minority" = "People of Color",
                        "RI Minority" = "RI POC",
                        "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(ri_wAvgWalkability_df[8,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 = 0.6, vjust=-0.8, size = 3,
            color=rgb(100,100,100, maxColorValue=255)) +
  geom_hline(yintercept = pull(ri_wAvgWalkability_df[8,2]), linetype = "dashed") +
  geom_text(aes(x = "Under 18", y = 13.5, label = "Above state\naverage"),
            color = "gray48") +
  geom_text(aes(x = "Under 18", y = 11.9, label = "Below state\naverage"),
            color = "gray48") +
  expand_limits(y = c(11.5,16))

ggsave("images/RI_Walk_graph.png")
```
```{r lowWalk, include=FALSE}
# Isolate least walkable block groups and tracts and add area variable
leastWalkable <- ri_blkgrp_walkability_sf %>% 
  filter(NatWalkInd <= 5.8) %>% 
  mutate(RI_LOWINC = if_else(RI_INCOME == "I",totalpopE,0)) %>% 
  mutate(RI_LOWINC = replace_na(RI_LOWINC,0)) %>% 
  mutate(RI_MINORITIES = if_else(RI_MINORITY == "M",totalpopE,0)) %>% 
  mutate(RI_MINORITIES = replace_na(RI_MINORITIES,0)) %>%
  mutate(SqKm = bg_area_m2/10^6)

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

leastWalkable1 <- leastWalkable %>%
  st_centroid(.) %>% 
  st_intersection(ri_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`))

leastWalkable2 <- leastWalkable %>%
  st_centroid(.) %>% 
  st_intersection(ri_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`))
```
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:mapLowWalkO64) and \@ref(fig:mapLowWalkU18) below highlight Census Block Groups across Rhode Island that are least walkable and the numbers of persons over 64 or under 18 live there. Tables \@ref(tab:tabLowWalkO64) and \@ref(tab:tabLowWalkU18) show breakdowns by municipality. For example, Table \@ref(tab:tabLowWalkO64) shows that `r formatC(as.numeric(leastWalkable1[1,3]),big.mark=",")` people over 64 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 over 64 in `r leastWalkable1[1,1]`. Similarly, Table \@ref(tab:tabLowWalkU18) shows that `r formatC(as.numeric(leastWalkable2[1,3]),big.mark=",")` people under 18 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 under 18 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 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(ri_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 = "quantile",
          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(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_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/RI/figures/ri_walkO64.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_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(ri_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`)) %>% 
  filter(`Over 64` > 0) %>% 
  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"))
```

```{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)."}
m <- tm_shape(ri_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 = "quantile",
          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(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_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/RI/figures/ri_walkU18.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_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(ri_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`)) %>% 
  filter(`Under 18` > 0) %>% 
  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"))
```



## 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 == "RI") %>% 
  st_make_valid() %>% 
  st_transform(., crs = 2840)

# Downscale LAI to Block Group by assigning tract value to all overlapping block groups (i.e. block groups with same tract GEOID prefix)
lai_ri_blkgrp <- ri_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_ri_tract <- ri_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_ri_tract_df <- lai_ri_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_ri_df <- lai_ri_blkgrp %>% 
  as.data.frame() %>% 
  mutate(RI_LOWINC = if_else(RI_INCOME == "I",totalpopE,0)) %>% 
  mutate(RI_LOWINC = replace_na(RI_LOWINC,0)) %>% 
  mutate(RI_MINORITIES = if_else(RI_MINORITY == "M",totalpopE,0)) %>% 
  mutate(RI_MINORITIES = replace_na(RI_MINORITIES,0)) %>%
  transmute(`Total Pop` = totalpopE,
            Minority = minorityE,
            `Low Income` = num2povE,
            `Under 5` = under5E,
            `Under 18` = under18E,
            `Over 64` = over64E,
            `Limited English HH` = eng_limitE,
            `No HS Dip` = lthsE,
            `RI Low Income` = RI_LOWINC,
            `RI Minority` = RI_MINORITIES,
            hh7_t) %>% 
  gather(key = Group, value = Pop, `Total Pop`:`RI Minority`) %>% 
  group_by(Group) %>% 
  summarize(wTBurden = weighted.mean(hh7_t, Pop, na.rm = TRUE)) %>% 
  rbind(., lai_ri_tract_df)

# save data for later analysis
save(lai_ri_blkgrp,
     lai_ri_tract, lai_ri_df,
     file = "DATA/transport/RI/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 Rhode Island. Transportation cost burden varies across the state. It is lowest in and around the state's major urban centers, such as Providence, Pawtucket, Woonsocket, and Newport, and highest in most of the rest of the state. In `r nrow(lai_ToverH)` Census Tract, the transportation cost burden exceeds the cost burden of housing. 

```{r mapLAI, fig.align = "center", fig.cap="Transportation cost burden for Census Tracts across Rhode Island."}
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(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_add_legend(type = "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/RI/figures/ri_tburden.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_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 that are below that of the state average of `r round(lai_ri_df %>% filter(Group == "Total Pop") %>% select(wTBurden) %>% pull(),1)`%. However, persons over age 64 tend to live in communities with transportation cost burdens that come to almost `r round(lai_ri_df %>% filter(Group == "Over 64") %>% select(wTBurden) %>% pull(),0)`% of household income. 

```{r loliLAI, fig.align = "center", fig.cap="Population-weighted average transportation cost burden as a percentage of household income for groups across Rhode Island."}
# Create lollipop graph of pop-weighted average LAI
lai_ri_df %>% 
  mutate(Group = recode(Group, "Minority" = "People of Color",
                        "RI Minority" = "RI POC",
                        "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_ri_df[8,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_ri_df[8,2]), linetype = "dashed") +
  geom_text(aes(x = "Under 18", y = 24.5, label = "Above state\naverage"),
            color = "gray48") +
  geom_text(aes(x = "Under 18", y = 21.5, label = "Below state\naverage"),
            color = "gray48") +
  expand_limits(y = c(18,25))

ggsave("images/RI_Cost_graph.png")
```
```{r HiBurden, include=FALSE}
# Isolate highest transport burden block groups and tracts and add area variable
HiBurden <- lai_ri_blkgrp %>% 
  filter(percent_rank(hh7_t) > 0.8) %>% 
  mutate(RI_LOWINC = if_else(RI_INCOME == "I",totalpopE,0)) %>% 
  mutate(RI_LOWINC = replace_na(RI_LOWINC,0)) %>% 
  mutate(RI_MINORITIES = if_else(RI_MINORITY == "M",totalpopE,0)) %>% 
  mutate(RI_MINORITIES = replace_na(RI_MINORITIES,0)) %>%
  mutate(SqKm = bg_area_m2/10^6)

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

HiBurden1 <- HiBurden %>%
  st_centroid(.) %>% 
  st_intersection(ri_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`))

HiBurden2 <- HiBurden %>%
  st_centroid(.) %>% 
  st_intersection(ri_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`))
```
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:mapHiBurdenU18) below highlight Census Block Groups across Rhode Island where the transportation cost burden is in the highest quintile (i.e., 80th percentile) and the numbers of persons over 64 or under 18 that live there. Tables \@ref(tab:tabHiBurdenO64) and \@ref(tab:tabHiBurdenU18) 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:tabHiBurdenU18) shows that `r formatC(as.numeric(HiBurden2[1,3]),big.mark=",")` people under 18 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 people under 18 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(ri_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 = "quantile",
          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(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_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/RI/figures/ri_HiBurdenO64.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_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(ri_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`)) %>% 
  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 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(ri_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 = "quantile",
          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(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_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/RI/figures/ri_HiBurdenU18.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_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(ri_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`)) %>% 
  filter(`Under 18` > 0) %>% 
  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"))
```



[^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/


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

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

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

## RIPTA transit network {-}

The transit route linework and transit stop locations for the Rhode Island Public Transit Authority (RIPTA) were downloaded from the Rhode Island Geographic Information System (RIGIS) at http://www.rigis.org/. Populations with access to transit were assumed to reside within 400 meters or 1/4 mile of a transit stop. This distance is a standard metric used to assess physical accessibility to buses. 

To identify the populations with access to transit, a 400 meter buffer was generated around all transit stops. 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("GRID 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 GRID codes excluding 22 - 24; see Table \@ref(tab:NLCDtab)). These undeveloped areas of the NLCD polygon layer were then used to erase overlapping areas of Census Block Groups. The latter step served to eliminate areas of Census Block Groups that were classified as undeveloped and thus unlikely to be occupied by residences. The remaining developed areas of Census Block Groups were assumed to contain all populations assigned to that respective Block Group. Processing of NLCD and developed Census Block Groups was done in ArcGIS Desktop 10.7. These developed Census Block Group polygons were used for identifying populations with access to public transit.


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

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 Rhode Island Public Transit Authority (RIPTA) was acquired from RIPTA GTFS - OpenMobilityData - TransitFeeds at https://transitfeeds.com/p/rhode-island-public-transit-authority/363. 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 (`r formatC(nrow(ri_transit_stops_sf),big.mark=",")` stops in total) 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 of these stops. 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. 


## 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 Rhode Island, 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., Minority pop x Headway) is then divided by the sum of the weights (i.e., total Minority population), so that 3475/45 = `r round(3475/45,2)`. The result is a weighted average 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 mapTransitAccessU18, fig.align = "center", fig.cap="Persons under 18 who are one or more miles from the nearest transit stop by Census Block Group."}
m <- tm_shape(ri_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(ri_blkgrps_sf_noBus, unit = "mi") +
  tm_fill(col = "NewUnder18", alpha = 0.9, style = "equal",
          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(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(ri_busroutes_sf) +
  tm_lines(lwd = 2, alpha = 0.5) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_add_legend(type = "line", lwd=2, alpha = 0.5, labels = "Transit Route") +
  tm_layout(title = "Persons Under 18\n1+ Miles from\nNearest Transit\nStop", 
            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/RI/figures/ri_transitU18access.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_transitU18access.png")
```

```{r tabTransitAccessU18, fig.align = "center", fig.cap="Persons under 18 in Census Block Groups one or more miles from the nearest bus stop."}
ri_blkgrps_sf_noBus %>%
  st_centroid(.) %>% 
  st_intersection(ri_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`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Persons Under 18 One or More Miles from Nearest Transit Stop", 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 mapTransitAccessU5, fig.align = "center", fig.cap="Children under 5 who are one or more miles from the nearest transit stop by Census Block Group."}
m <- tm_shape(ri_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(ri_blkgrps_sf_noBus, unit = "mi") +
  tm_fill(col = "NewUnder5", alpha = 0.9, style = "equal",
          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(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(ri_busroutes_sf) +
  tm_lines(lwd = 2, alpha = 0.5) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_add_legend(type = "line", lwd=2, alpha = 0.5, labels = "Transit Route") +
  tm_layout(title = "Children Under 5\n1+ Miles from\nNearest Transit\nStop", 
            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/RI/figures/ri_transitU5access.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_transitU5access.png")
```

```{r tabTransitAccessU5, fig.align = "center", fig.cap="Children under 5 in Census Block Groups one or more miles from the nearest bus stop."}
ri_blkgrps_sf_noBus %>%
  st_centroid(.) %>% 
  st_intersection(ri_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`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Children Under 5 One or More Miles from Nearest Transit Stop", 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 mapTransitAccessNOHS, fig.align = "center", fig.cap="Adults without a high school diploma who are one or more miles from the nearest transit stop by Census Block Group."}
m <- tm_shape(ri_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(ri_blkgrps_sf_noBus, unit = "mi") +
  tm_fill(col = "NewLths", alpha = 0.9, style = "equal",
          colorNA = "white", colorNULL = "gray", title = "Number of\nAdults without a\nHigh School Diploma",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd") +
  tm_shape(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  # tm_shape(ri_busroutes_sf) +
  # tm_lines(col = "type", palette = c("purple","red","blue"),
  #          colorNA = NULL, title.col = "Route Type", lwd = 2,
  #          alpha = 0.5) +
  tm_shape(ri_busroutes_sf) +
  tm_lines(lwd = 2, alpha = 0.5) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_add_legend(type = "line", lwd=2, alpha = 0.5, labels = "Transit Route") +
  tm_layout(title = "Adults without a\nHigh School\nDiploma\n1+ Miles from\nNearest Transit\nStop", 
            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/RI/figures/ri_transitNOHSaccess.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_transitNOHSaccess.png")
```

```{r tabTransitAccessNOHS, fig.align = "center", fig.cap="Adults without a high school diploma in Census Block Groups one or more miles from the nearest bus stop."}
ri_blkgrps_sf_noBus %>%
  st_centroid(.) %>% 
  st_intersection(ri_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) %>% 
  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 One or More Miles from Nearest Transit Stop", 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 mapTransitAccessINC, fig.align = "center", fig.cap="Low income persons who are one or more miles from the nearest transit stop by Census Block Group."}
m <- tm_shape(ri_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(ri_blkgrps_sf_noBus, unit = "mi") +
  tm_fill(col = "NewPov", alpha = 0.9, style = "equal",
          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(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  # tm_shape(ri_busroutes_sf) +
  # tm_lines(col = "type", palette = c("purple","red","blue"),
  #          colorNA = NULL, title.col = "Route Type", lwd = 2,
  #          alpha = 0.5) +
  tm_shape(ri_busroutes_sf) +
  tm_lines(lwd = 2, alpha = 0.5) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_add_legend(type = "line", lwd=2, alpha = 0.5, labels = "Transit Route") +
  tm_layout(title = "Low Income\nPersons\n1+ Miles from\nNearest Transit\nStop", 
            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/RI/figures/ri_transitINCaccess.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_transitINCaccess.png")
```

```{r tabTransitAccessINC, fig.align = "center", fig.cap="Low income persons in Census Block Groups one or more miles from the nearest bus stop."}
ri_blkgrps_sf_noBus %>%
  st_centroid(.) %>% 
  st_intersection(ri_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) %>% 
  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 One or More Miles from Nearest Transit Stop", 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 mapTransitAccessNoCar, fig.align = "center", fig.cap="Households without a car that are one or more miles from the nearest transit stop by Census Tract."}
m <- tm_shape(ri_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(ri_tracts_sf_noBus, unit = "mi") +
  tm_fill(col = "NewHHNoCar", alpha = 0.9, style = "equal",
          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(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  # tm_shape(ri_busroutes_sf) +
  # tm_lines(col = "type", palette = c("purple","red","blue"),
  #          colorNA = NULL, title.col = "Route Type", lwd = 2,
  #          alpha = 0.5) +
  tm_shape(ri_busroutes_sf) +
  tm_lines(lwd = 2, alpha = 0.5) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_add_legend(type = "line", lwd=2, alpha = 0.5, labels = "Transit Route") +
  tm_layout(title = "Households\nwithout a Car\n1+ Miles from\nNearest Transit\nStop", 
            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/RI/figures/ri_transitNoCaraccess.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_transitNoCaraccess.png")
```

```{r tabTransitAccessNoCar, fig.align = "center", fig.cap="Households without a car in Census Tracts one or more miles from the nearest bus stop."}
ri_tracts_sf_noBus %>%
  st_centroid(.) %>% 
  st_intersection(ri_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Tracts` = n(),
            `No Car HH` = sum(NewHHNoCar,na.rm = TRUE),
            `Pct of No Car HH` = round(sum(NewHHNoCar)/max(TotalNoCar)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  mutate(`Pct of No Car HH` = paste0(`Pct of No Car HH`,"%")) %>% 
  arrange(desc(`No Car HH`)) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Households without a Car One or More Miles from Nearest Transit Stop", 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"))
```

```{r mapTransitAccessMIN, fig.align = "center", fig.cap="People of Color who are one or more miles from the nearest transit stop by Census Block Group."}
m <- tm_shape(ri_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(ri_blkgrps_sf_noBus, unit = "mi") +
  tm_fill(col = "NewMinority", alpha = 0.9, style = "equal",
          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(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  # tm_shape(ri_busroutes_sf) +
  # tm_lines(col = "type", palette = c("purple","red","blue"),
  #          colorNA = NULL, title.col = "Route Type", lwd = 2,
  #          alpha = 0.5) +
  tm_shape(ri_busroutes_sf) +
  tm_lines(lwd = 2, alpha = 0.5) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_add_legend(type = "line", lwd=2, alpha = 0.5, labels = "Transit Route") +
  tm_layout(title = "People of Color\n1+ Miles from\nNearest Transit\nStop", 
            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/RI/figures/ri_transitMINaccess.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_transitMINaccess.png")
```

```{r tabTransitAccessMIN, fig.align = "center", fig.cap="People of Color in Census Block Groups one or more miles from the nearest bus stop."}
ri_blkgrps_sf_noBus %>%
  st_centroid(.) %>% 
  st_intersection(ri_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            Minority = sum(NewOver64,na.rm = TRUE),
            `Pct of Minority` = round(sum(NewMinority)/max(TotalMin)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  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 One or More Miles from Nearest Transit Stop", 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 mapTransitAccessLEH, fig.align = "center", fig.cap="Limited English speaking households that are one or more miles from the nearest transit stop by Census Block Group."}
m <- tm_shape(ri_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(ri_blkgrps_sf_noBus, unit = "mi") +
  tm_fill(col = "NewEng_limit", alpha = 0.9, style = "equal",
          colorNA = "white", colorNULL = "gray", title = "Number of\nLimited English Households",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd") +
  tm_shape(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  # tm_shape(ri_busroutes_sf) +
  # tm_lines(col = "type", palette = c("purple","red","blue"),
  #          colorNA = NULL, title.col = "Route Type", lwd = 2,
  #          alpha = 0.5) +
  tm_shape(ri_busroutes_sf) +
  tm_lines(lwd = 2, alpha = 0.5) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_add_legend(type = "line", lwd=2, alpha = 0.5, labels = "Transit Route") +
  tm_layout(title = "Limited English\nSpeaking\nHouseholds\n1+ Miles from\nNearest Transit\nStop", 
            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/RI/figures/ri_transitLEHaccess.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_transitLEHaccess.png")
```

```{r tabTransitAccessLEH, fig.align = "center", fig.cap="Limited English speaking households in Census Block Groups one or more miles from the nearest bus stop."}
ri_blkgrps_sf_noBus %>%
  st_centroid(.) %>% 
  st_intersection(ri_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Limited English HH` = sum(NewEng_limit,na.rm = TRUE),
            `Pct of Limited English HH` = round(sum(NewEng_limit)/max(TotalLEH)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  mutate(`Pct of Limited English HH` = paste0(`Pct of Limited English HH`,"%")) %>% 
  arrange(desc(`Limited English HH`)) %>% 
  filter(`Limited English HH` > 0) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Limited English Speaking Households One or More Miles from Nearest Transit Stop", 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")) 
```

Note that no Low Income or Minority communities, as defined by the state's environmental justice policy, were found to fall outside of areas with access to transit stops. 


## Excessive Headways and Priority Populations {-}

```{r mapHeadwayUnder18, 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(ri_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(headway80thblkgrps) +
  tm_fill(col = "NewUnder18", alpha = 0.9, style = "equal",
          colorNA = "white", colorNULL = "gray", title = "Number of Persons\nUnder 18",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd") +
  tm_shape(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  # tm_shape(ri_busroutes_sf) +
  # tm_lines(col = "type", palette = c("purple","red","blue"),
  #          colorNA = NULL, title.col = "Route Type", lwd = 2,
  #          alpha = 0.5) +
  tm_shape(ri_busroutes_sf) + tm_lines(lwd = 1, alpha = 0.3) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_add_legend(type = "line", lwd=1, alpha = 0.3, labels = "Transit Route") +
  tm_layout(title = paste("Persons Under 18","with Headways",paste0("Over ",round(quantile(ri_routes_sf$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/RI/figures/ri_headwayU18.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_headwayU18.png")
```

```{r tabHeadwayUnder18, 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 %>%
  st_centroid(.) %>% 
  st_intersection(ri_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`)) %>% 
  filter(`Under 18` > 0) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = paste("Persons Under 18","with Headways",paste0("Over ",quantile(ri_routes_sf$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"))

```

```{r mapHeadwayNoCar, 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(ri_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(headway80thtracts) +
  tm_fill(col = "NewHHNoCar", alpha = 0.9, style = "equal",
          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(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  # tm_shape(ri_busroutes_sf) +
  # tm_lines(col = "type", palette = c("purple","red","blue"),
  #          colorNA = NULL, title.col = "Route Type", lwd = 2,
  #          alpha = 0.5) +
  tm_shape(ri_busroutes_sf) + tm_lines(lwd = 1, alpha = 0.3) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_add_legend(type = "line", lwd=1, alpha = 0.3, labels = "Transit Route") +
  tm_layout(title = paste("Households","without a Car","with Headways",paste0("Over ",round(quantile(ri_routes_sf$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/RI/figures/ri_headwayNoCar.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_headwayNoCar.png")
```

```{r tabHeadwayNoCar, 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 %>%
  st_centroid(.) %>% 
  st_intersection(ri_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Tracts` = n(),
            `No Car HH` = sum(NewHHNoCar,na.rm = TRUE),
            `Pct of No Car HH` = round(sum(NewHHNoCar)/max(TotalNoCar)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  mutate(`Pct of No Car HH` = paste0(`Pct of No Car HH`,"%")) %>% 
  arrange(desc(`No Car HH`)) %>% 
  filter(`No Car HH` > 0) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = paste("Households without a Car","with Headways",paste0("Over ",quantile(ri_routes_sf$mean_headways,0.8)," Minutes"),sep = "\n"), 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"))

```

```{r mapHeadwayUnder5, 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(ri_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(headway80thblkgrps) +
  tm_fill(col = "NewUnder5", alpha = 0.9, style = "equal",
          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(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  # tm_shape(ri_busroutes_sf) +
  # tm_lines(col = "type", palette = c("purple","red","blue"),
  #          colorNA = NULL, title.col = "Route Type", lwd = 2,
  #          alpha = 0.5) +
  tm_shape(ri_busroutes_sf) + tm_lines(lwd = 1, alpha = 0.3) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_add_legend(type = "line", lwd=1, alpha = 0.3, labels = "Transit Route") +
  tm_layout(title = paste("Children Under 5","with Headways",paste0("Over ",round(quantile(ri_routes_sf$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/RI/figures/ri_headwayUnder5.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_headwayUnder5.png")
```

```{r tabHeadwayUnder5, 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 %>%
  st_centroid(.) %>% 
  st_intersection(ri_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`)) %>% 
  filter(`Under 5` > 0) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = paste("Children Under 5","with Headways",paste0("Over ",quantile(ri_routes_sf$mean_headways,0.8)," 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 mapHeadwayLowInc, 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(ri_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(headway80thblkgrps) +
  tm_fill(col = "NewPov", alpha = 0.9, style = "equal",
          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(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  # tm_shape(ri_busroutes_sf) +
  # tm_lines(col = "type", palette = c("purple","red","blue"),
  #          colorNA = NULL, title.col = "Route Type", lwd = 2,
  #          alpha = 0.5) +
  tm_shape(ri_busroutes_sf) + tm_lines(lwd = 1, alpha = 0.3) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_add_legend(type = "line", lwd=1, alpha = 0.3, labels = "Transit Route") +
  tm_layout(title = paste("Low Income", "Persons with","Headways Over",paste0(round(quantile(ri_routes_sf$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/RI/figures/ri_headwayLowInc.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_headwayLowInc.png")
```

```{r tabHeadwayLowInc, 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 %>%
  st_centroid(.) %>% 
  st_intersection(ri_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) %>% 
  mutate(`Pct of Low Income` = paste0(`Pct of Low Income`,"%")) %>% 
  arrange(desc(`Low Income`)) %>% 
  filter(`Low Income` > 0) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = paste("Low Income Persons","with Headways",paste0("Over ",quantile(ri_routes_sf$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 mapHeadwayNoHSDip, 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(ri_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(headway80thblkgrps) +
  tm_fill(col = "NewLths", alpha = 0.9, style = "equal",
          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(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  # tm_shape(ri_busroutes_sf) +
  # tm_lines(col = "type", palette = c("purple","red","blue"),
  #          colorNA = NULL, title.col = "Route Type", lwd = 2,
  #          alpha = 0.5) +
  tm_shape(ri_busroutes_sf) + tm_lines(lwd = 1, alpha = 0.3) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_add_legend(type = "line", lwd=1, alpha = 0.3, labels = "Transit Route") +
  tm_layout(title = paste("Adults without a","High School","Diploma with", "Headways Over",paste0(round(quantile(ri_routes_sf$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/RI/figures/ri_headwayNoHSDip.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_headwayNoHSDip.png")
```

```{r tabHeadwayNoHSDip, 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 %>%
  st_centroid(.) %>% 
  st_intersection(ri_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) %>% 
  mutate(`Pct of No HS Dip` = paste0(`Pct of No HS Dip`,"%")) %>% 
  arrange(desc(`No HS Dip`)) %>% 
  filter(`No HS Dip` > 0) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = paste("Adults without a","High School Diploma","with Headways",paste0("Over ",quantile(ri_routes_sf$mean_headways,0.8)," 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 mapHeadwayMinority, 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(ri_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(headway80thblkgrps) +
  tm_fill(col = "NewMinority", alpha = 0.9, style = "equal",
          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(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  # tm_shape(ri_busroutes_sf) +
  # tm_lines(col = "type", palette = c("purple","red","blue"),
  #          colorNA = NULL, title.col = "Route Type", lwd = 2,
  #          alpha = 0.5) +
  tm_shape(ri_busroutes_sf) + tm_lines(lwd = 1, alpha = 0.3) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_add_legend(type = "line", lwd=1, alpha = 0.3, labels = "Transit Route") +
  tm_layout(title = paste("People of Color","with Headways",paste0("Over ",round(quantile(ri_routes_sf$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/RI/figures/ri_headwayMinority.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_headwayMinority.png")
```

```{r tabHeadwayMinority, 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 %>%
  st_centroid(.) %>% 
  st_intersection(ri_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) %>% 
  mutate(`Pct of Minority` = paste0(`Pct of Minority`,"%")) %>% 
  arrange(desc(`Minority`)) %>% 
  filter(`Minority` > 0) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = paste("People of Color","with Headways",paste0("Over ",quantile(ri_routes_sf$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 mapHeadwayLEH, 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(ri_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(headway80thblkgrps) +
  tm_fill(col = "NewEng_limit", alpha = 0.9, style = "equal",
          colorNA = "white", colorNULL = "gray", title = "Number of\nLimited English Households",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd") +
  tm_shape(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  # tm_shape(ri_busroutes_sf) +
  # tm_lines(col = "type", palette = c("purple","red","blue"),
  #          colorNA = NULL, title.col = "Route Type", lwd = 2,
  #          alpha = 0.5) +
  tm_shape(ri_busroutes_sf) + tm_lines(lwd = 1, alpha = 0.3) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_add_legend(type = "line", lwd=1, alpha = 0.3, labels = "Transit Route") +
  tm_layout(title = paste("Limited English","Speaking", "Households","with Headways",paste0("Over ",round(quantile(ri_routes_sf$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/RI/figures/ri_headwayLEH.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_headwayLEH.png")
```

```{r tabHeadwayLEH, 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 %>%
  st_centroid(.) %>% 
  st_intersection(ri_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Limited English HH` = sum(NewEng_limit,na.rm = TRUE),
            `Pct of Limited English HH` = round(sum(NewEng_limit)/max(TotalLEH)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  mutate(`Pct of Limited English HH` = paste0(`Pct of Limited English HH`,"%")) %>% 
  arrange(desc(`Limited English HH`)) %>% 
  filter(`Limited English HH` > 0) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = paste("Limited English","Speaking Households","with Headways",paste0("Over ",quantile(ri_routes_sf$mean_headways,0.8)," Minutes"),sep = "\n"), 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 mapHeadwayRILowInc, fig.align = "center", fig.cap="Low income persons, as defined by the state's environmental justice policy, 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(ri_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(headway80thblkgrps) +
  tm_fill(col = "NewRI_LOWINC", alpha = 0.9, style = "equal",
          colorNA = "white", colorNULL = "gray", title = "Number of\nRI Low Income",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f", big.mark = ",")), palette = "OrRd") +
  tm_shape(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  # tm_shape(ri_busroutes_sf) +
  # tm_lines(col = "type", palette = c("purple","red","blue"),
  #          colorNA = NULL, title.col = "Route Type", lwd = 2,
  #          alpha = 0.5) +
  tm_shape(ri_busroutes_sf) + tm_lines(lwd = 1, alpha = 0.3) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_add_legend(type = "line", lwd=1, alpha = 0.3, labels = "Transit Route") +
  tm_layout(title = paste("RI Low Income", "Persons with","Headways Over",paste0(round(quantile(ri_routes_sf$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/RI/figures/ri_headwayRILowInc.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_headwayRILowInc.png")
```

```{r tabHeadwayRILowInc, fig.align = "center", fig.cap="Low income persons, as defined by the state's environmental justice policy, 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 %>%
  st_centroid(.) %>% 
  st_intersection(ri_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `RI Low Income` = sum(NewRI_LOWINC,na.rm = TRUE),
            `Pct of RI Low Income` = round(sum(NewRI_LOWINC)/max(TotalRI_LOWINC)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  mutate(`Pct of RI Low Income` = paste0(`Pct of RI Low Income`,"%")) %>% 
  arrange(desc(`RI Low Income`)) %>% 
  filter(`RI Low Income` > 0) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = paste("RI Low Income Persons","with Headways",paste0("Over ",quantile(ri_routes_sf$mean_headways,0.8)," Minutes"),sep = "\n"), align = "r", col.names = c(names(.)[1:2],"Number of RI Low Income in Block Groups","Pct of RI Low Income in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>% 
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))

```

```{r mapHeadwayRIMinority, fig.align = "center", fig.cap="People of Color, as defined by the state's environmental justice policy, 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(ri_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(headway80thblkgrps) +
  tm_fill(col = "NewRI_MINORITIES", alpha = 0.9, style = "equal",
          showNA = FALSE, title = "Number of\nRI POC",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f", big.mark = ",")), palette = "OrRd") +
  tm_shape(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  # tm_shape(ri_busroutes_sf) +
  # tm_lines(col = "type", palette = c("purple","red","blue"),
  #          colorNA = NULL, title.col = "Route Type", lwd = 2,
  #          alpha = 0.5) +
  tm_shape(ri_busroutes_sf) + tm_lines(lwd = 1, alpha = 0.3) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_add_legend(type = "line", lwd=1, alpha = 0.3, labels = "Transit Route") +
  tm_layout(title = paste("RI People", "of Color with","Headways Over",paste0(round(quantile(ri_routes_sf$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/RI/figures/ri_headwayRIMinority.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_headwayRIMinority.png")
```

```{r tabHeadwayRIMinority, fig.align = "center", fig.cap="People of Color, as defined by the state's environmental justice policy, 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 %>%
  st_centroid(.) %>% 
  st_intersection(ri_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `RI Minority` = sum(NewRI_MINORITIES,na.rm = TRUE),
            `Pct of RI Minority` = round(sum(NewRI_MINORITIES)/max(TotalRI_MINORITY)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  mutate(`Pct of RI Minority` = paste0(`Pct of RI Minority`,"%")) %>% 
  arrange(desc(`RI Minority`)) %>% 
  filter(`RI Minority` > 0) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = paste("RI People of Color","with Headways",paste0("Over ",quantile(ri_routes_sf$mean_headways,0.8)," Minutes"),sep = "\n"), align = "r", col.names = c(names(.)[1:2],"Number of RI POC in Block Groups","Pct of RI POC 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 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(ri_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 = "equal",
          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(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "Disabled Persons\nin Least Walkable\nTracts", 
            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/RI/figures/ri_walkDisabled.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_walkDisabled.png")
```

```{r tabLowWalkDisabled, fig.align = "center", fig.cap="Disabled persons by municipality in least walkable Census Tracts (i.e., Walkability Index Score less than or equal to 5.8)."}
leastWalkableTract %>%
  st_centroid(.) %>% 
  st_intersection(ri_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) %>% 
  mutate(`Pct of Disabled` = paste0(`Pct of Disabled`,"%")) %>% 
  arrange(desc(`Disabled`)) %>% 
  filter(Disabled > 0) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Disabled Persons in Least Walkable 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 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(ri_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 = "equal",
          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(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_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/RI/figures/ri_walkU5.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_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(ri_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) %>% 
  mutate(`Pct of Under 5` = paste0(`Pct of Under 5`,"%")) %>% 
  arrange(desc(`Under 5`)) %>% 
  filter(`Under 5` > 0) %>% 
  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 mapLowWalkNoHSDip, 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(ri_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 = "equal",
          colorNA = "white", colorNULL = "gray", title = "No HS Dip\nPer Square Kilometer",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd", convert2density = TRUE, area = "SqKm") +
  tm_shape(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "Adults without a\nHigh School\nDiploma in\nLeast 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/RI/figures/ri_walkNoHSDip.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_walkNoHSDip.png")
```

```{r tabLowWalkNoHSDip, 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(ri_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) %>% 
  mutate(`Pct of No HS Dip` = paste0(`Pct of No HS Dip`,"%")) %>% 
  arrange(desc(`No HS Dip`)) %>% 
  filter(`No HS Dip` > 0) %>% 
  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 mapLowWalkLowInc, 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(ri_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 = "equal",
          colorNA = "white", colorNULL = "gray", title = "Low Income Persons\nPer Square Kilometer",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd", convert2density = TRUE, area = "SqKm") +
  tm_shape(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "Low Income\nPersons in\nLeast 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/RI/figures/ri_walkLowInc.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_walkLowInc.png")
```

```{r tabLowWalkLowInc, 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(ri_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) %>% 
  mutate(`Pct of Low Income` = paste0(`Pct of Low Income`,"%")) %>% 
  arrange(desc(`Low Income`)) %>% 
  filter(`Low Income` > 0) %>% 
  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 mapLowWalkNoCar, fig.align = "center", fig.cap="Households without a car per square kilometer in least walkable Census Tracts (i.e., Walkability Index Score less than or equal to 5.8)."}
m <- tm_shape(ri_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 = "equal",
          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(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "No Car Households\nin Least Walkable\nTracts", 
            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/RI/figures/ri_walkNoCar.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_walkNoCar.png")
```

```{r tabLowWalkNoCar, fig.align = "center", fig.cap="Households without a car by municipality in least walkable Census Tracts (i.e., Walkability Index Score less than or equal to 5.8)."}
leastWalkableTract %>%
  st_centroid(.) %>% 
  st_intersection(ri_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Tracts` = n(),
            Disabled = sum(HHnoCarE,na.rm = TRUE),
            `Pct of Disabled` = round(sum(HHnoCarE)/max(TotalNoCar)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  mutate(`Pct of Disabled` = paste0(`Pct of Disabled`,"%")) %>% 
  arrange(desc(`Disabled`)) %>% 
  filter(Disabled > 0) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Households without a Car in Least Walkable Census Tracts", 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"))
```

```{r mapLowWalkMinority, 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(ri_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 = "equal",
          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(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_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/RI/figures/ri_walkMinority.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_walkMinority.png")
```

```{r tabLowWalkMinority, 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(ri_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) %>% 
  mutate(`Pct of Minority` = paste0(`Pct of Minority`,"%")) %>% 
  arrange(desc(Minority)) %>% 
  filter(Minority > 0) %>% 
  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(ri_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 = "equal",
          colorNA = "white", colorNULL = "gray", title = "Limited English Households\nPer Square Kilometer",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd", convert2density = TRUE, area = "SqKm") +
  tm_shape(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "Limited English\nSpeaking\nHouseholds in\nLeast 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/RI/figures/ri_walkLEH.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_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(ri_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Limited English HH` = sum(eng_limitE,na.rm = TRUE),
            `Pct of Limited English HH` = round(sum(eng_limitE)/max(TotalLEH)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  mutate(`Pct of Limited English HH` = paste0(`Pct of Limited English HH`,"%")) %>% 
  arrange(desc(`Limited English HH`)) %>% 
  filter(`Limited English HH` > 0) %>% 
  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"))
```

Note that no low income and minority communities, as defined by the state's environmental justice policy, were classified as least walkable communities. 


## Highest Transportation Cost Burden and Priority Populations {-}

```{r mapHiBurdenDisabled, fig.align = "center", fig.cap="Disabled Persons per square kilometer in Census Tracts with the highest transportation cost burden (i.e., transportation cost burden in the 80th percentile for the state)."}
m <- tm_shape(ri_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 = "equal",
          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(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "Disabled Persons\nin Highest\nTransportation\nCost Burden\nCensus Tracts", 
            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/RI/figures/ri_HiBurdenDisabled.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_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(ri_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) %>% 
  mutate(`Pct of Disabled` = paste0(`Pct of Disabled`,"%")) %>% 
  arrange(desc(Disabled)) %>% 
  filter(Disabled > 0) %>% 
  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(ri_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 = "equal",
          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(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_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/RI/figures/ri_HiBurdenU5.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_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(ri_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) %>% 
  mutate(`Pct of Under 5` = paste0(`Pct of Under 5`,"%")) %>% 
  arrange(desc(`Under 5`)) %>% 
  filter(`Under 5` > 0) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Persons 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 mapHiBurdenLTHS, fig.align = "center", fig.cap="Adults wihtout 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(ri_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 = "equal",
          colorNA = "white", colorNULL = "gray", title = "Adults without a\nHigh School Diploma\nPer Square Kilometer",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd", convert2density = TRUE, area = "SqKm") +
  tm_shape(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_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/RI/figures/ri_HiBurdenLTHS.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_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(ri_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) %>% 
  mutate(`Pct of No HS Dip` = paste0(`Pct of No HS Dip`,"%")) %>% 
  arrange(desc(`No HS Dip`)) %>% 
  filter(`No HS Dip` > 0) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Adults without a High Schol 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"))
```

```{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(ri_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 = "equal",
          colorNA = "white", colorNULL = "gray", title = "Low Income Persons\nPer Square Kilometer",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd", convert2density = TRUE, area = "SqKm") +
  tm_shape(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_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/RI/figures/ri_HiBurdenINC.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_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(ri_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) %>% 
  mutate(`Pct of Low Income` = paste0(`Pct of Low Income`,"%")) %>% 
  arrange(desc(`Low Income`)) %>% 
  filter(`Low Income` > 0) %>% 
  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 Tracts with the highest transportation cost burden (i.e., transportation cost burden in the 80th percentile for the state)."}
m <- tm_shape(ri_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 = "equal",
          colorNA = "white", colorNULL = "gray", title = "Households\nwithout a Car\nPer Square Kilometer",
          legend.format = list(fun = function(x) formatC(x, digits = 0, format = "f")), palette = "OrRd", convert2density = TRUE, area = "SqKm") +
  tm_shape(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "Households\nwithout a Car\nin Highest\nTransportation\nCost Burden\nCensus Tracts", 
            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/RI/figures/ri_HiBurdenNoCarHH.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_HiBurdenNoCarHH.png")
```

```{r tabHiBurdenNoCarHH, fig.align = "center", fig.cap="Households without a car 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(ri_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Tracts` = 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) %>% 
  mutate(`Pct of No Car Households` = paste0(`Pct of No Car Households`,"%")) %>% 
  arrange(desc(`No Car Households`)) %>% 
  filter(`No Car Households` > 0) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Households without a Car in Highest Transportation Cost Burden Census Tracts", 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"))
```

```{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(ri_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 = "equal",
          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(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "Minority\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/RI/figures/ri_HiBurdenMIN.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_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(ri_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) %>% 
  mutate(`Pct of Minority` = paste0(`Pct of Minority`,"%")) %>% 
  arrange(desc(Minority)) %>% 
  filter(`Minority` > 0) %>% 
  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(ri_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 = "equal",
          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(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "Limited English\nSpeaking\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/RI/figures/ri_HiBurdenLEH.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_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(ri_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `Limited English HH` = sum(eng_limitE,na.rm = TRUE),
            `Pct of Limited English HH` = round(sum(eng_limitE)/max(TotalNoHS)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  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`)) %>% 
  filter(`Limited English HH` > 0) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "Limited English Speaking Households in Highest Transportation Cost Burden Block Groups", align = "r", col.names = c(names(.)[1:2],"Number of Limited English HH in 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 mapHiBurdenRIINC, fig.align = "center", fig.cap="Low income persons, as defined by state environmental justice policy, 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(ri_state_sf_cb, unit = "mi") + 
  tm_fill(col = "gray",alpha = 0.2) +
  tm_shape(HiBurden, unit = "mi") +
  tm_fill(col = "RI_LOWINC", alpha = 0.9, style = "equal",
          colorNA = "white", colorNULL = "gray", title = "RI 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(ri_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(ri_highways) + tm_lines(col = "seashell4", lwd = 1) +
  tm_shape(I95roadSegment) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I95roadSegment2) +
  tm_symbols(shape = I95, border.lwd = NA, size = 0.1) +
  tm_shape(I195roadSegment) +
  tm_symbols(shape = I195, border.lwd = NA, size = 0.1) +
  tm_shape(I295roadSegment) +
  tm_symbols(shape = I295, border.lwd = NA, size = 0.1) +
  tm_shape(ri_towns_sf_pts) + tm_dots() +
  tm_text("NAME", size = 0.5, col = "black", 
          xmod = 0.7, ymod = 0.2, shadow = TRUE) +
  tm_scale_bar(breaks = c(0, 5, 10), text.size = 0.5, 
               position = c(0.6,0.005)) +
  tm_layout(title = "RI 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/RI/figures/ri_HiBurdenRIINC.png",
          height = 8, width = 8, units = "in", dpi = 600)

knitr::include_graphics("DATA/transport/RI/figures/ri_HiBurdenRIINC.png")
```

```{r tabHiBurdenRIINC, fig.align = "center", fig.cap="Low income persons, as defined by state environmental justice policy, 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(ri_towns_sf,.) %>% 
  as.data.frame() %>% 
  left_join(., townpopsdf,by="NAME") %>%
  group_by(NAME) %>%
  summarize(`Block Groups` = n(),
            `RI Low Income` = sum(RI_LOWINC,na.rm = TRUE),
            `Pct of RI Low Income` = round(sum(RI_LOWINC)/max(TotalRI_LOWINC)*100,1)) %>% 
  rename("City/Town" = NAME) %>% 
  mutate(`Pct of RI Low Income` = if_else(`Pct of RI Low Income` > 100, 100,`Pct of RI Low Income`)) %>% 
  mutate(`Pct of RI Low Income` = paste0(`Pct of RI Low Income`,"%")) %>% 
  arrange(desc(`RI Low Income`)) %>% 
  filter(`RI Low Income` > 0) %>% 
  kableExtra::kable(longtable = T, booktabs = T,
                    format.args = list(big.mark = ','), 
                    digits = 1,
                    caption = "RI Low Income in Highest Transportation Cost Burden Block Groups", align = "r", col.names = c(names(.)[1:2],"Number of RI Low Income in Block Groups","Pct of RI Low Income in City/Town")) %>% 
  kableExtra::column_spec(3:4, width = "4cm") %>%
  kableExtra::kable_styling(latex_options = c("repeat_header","striped"))
```

Note that no minority communities, as defined by the state's environmental justice policy, were classified as highest transportation cost burden communities.