Route Analysis

date()

## [1] "Thu Oct 31 16:37:36 2024"

rm(list=ls())
library(tidyverse)
library(ggmap)
library(sf)
library(osrm)
library(smoothr)
library(magick)
library(ggnewscale)
library(rsvg)
fig.height <- 6
set.seed(1)

Main R script

## school focus
school_focus <- data.frame(name = c("East High School"), NCES_CODE = c("550852000925"))

## walk boundary
walk_boundary_m <- 1.5 * 1609

## load school locations

WI_schools <- st_transform(st_read(dsn = "data/Schools/Wisconsin_Public_Schools_-5986231931870160084.gpkg"), crs = 4326)

## Reading layer `WI_Public_Schools' from data source 
##   `/home/ben/Bike Fed/Documents/data_analysis/route_analysis/data/Schools/Wisconsin_Public_Schools_-5986231931870160084.gpkg' 
##   using driver `GPKG'
## replacing null geometries with empty geometries
## Simple feature collection with 2274 features and 32 fields (with 120 geometries empty)
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 298957.7 ymin: 225953.4 xmax: 760547.6 ymax: 704946.8
## Projected CRS: NAD83(HARN) / Wisconsin Transverse Mercator

WI_schools <- WI_schools %>% mutate(geom = SHAPE)

## load addresses
addresses <- read_csv(file="data/addresses/Addresses_Students_EastHS_2024_GeocodeResults.csv") %>%
  filter(lat > 0) %>%
  st_as_sf(coords=c("lon","lat"), crs=4326) # remember x=lon and y=lat

## set osrm options
options(osrm.server = "http://127.0.0.1:5000/")
options(osrm.profile = "walk")

register_stadiamaps(key = substr(read_file(file = "api_keys/stadia_api_key"), 1, 36))

## subset addresses within 1.5 miles
walk_boundary_poly <- fill_holes(st_make_valid(osrmIsodistance(
  loc = WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE),
  breaks = c(walk_boundary_m),
  res = 80)
), units::set_units(1, km^2))

addresses_near <- st_intersection(addresses, walk_boundary_poly)

## load bike tls
bike_lts <- st_read("data/bike_lts/bike_lts_DANE.geojson")

## Reading layer `Bike_LTS' from data source 
##   `/home/ben/Bike Fed/Documents/data_analysis/route_analysis/data/bike_lts/bike_lts_DANE.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 43355 features and 5 fields
## Geometry type: MULTILINESTRING
## Dimension:     XY
## Bounding box:  xmin: -89.84212 ymin: 42.8404 xmax: -89.0082 ymax: 43.2942
## Geodetic CRS:  WGS 84

bike_lts[["lts"]] <- as.factor(bike_lts$LTS_F)

bike_lts_scale <- data.frame(code = c(1, 2, 3, 4, 9),
                             color = c("#1a9641",
                                       "#a6d96a",
                                       "#fdae61",
                                       "#d7191c",
                                       "#d7191c"))

## calculate routes
routes <- list(NULL)
for(i in addresses_near$number) {
 routes[[i]] <- osrmRoute(
      src = addresses_near %>% filter(number == i),
      dst = WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE))
  message(paste0("done - ", i, "of", max(addresses_near$number)))
}
routes <- bind_rows(routes)

## combine routes
# Count the routes that intersect or overlap with each segment of the bike_tls network.
# The intersections have a buffer of 20m
bike_lts_buffer <- st_buffer(st_intersection(bike_lts, walk_boundary_poly), 20)

bike_lts_buffer["student_use"] <- unlist(lapply(st_intersects(bike_lts_buffer, routes), length))

bike_lts <- st_join(bike_lts, bike_lts_buffer %>% select(OBJECTID, student_use))

## make maps
# load logo
logo <- image_read(path = "other/BFW_Logo_180_x_200_transparent_background.png")
school_symbol <- image_read_svg(path = "other/school_FILL0_wght400_GRAD0_opsz24.svg")


bbox <- st_bbox(st_buffer(walk_boundary_poly, dist = 500))
bbox <- c(left = as.double(bbox[1]),
          bottom = as.double(bbox[2]),
          right = as.double(bbox[3]),
          top = as.double(bbox[4]))

#get basemap
basemap <- get_stadiamap(bbox = bbox, zoom = 15, maptype = "stamen_toner_lite")

# generate map
ggmap(basemap) +
  labs(title = paste0("Walking routes for students at ",
    school_focus %>% pull(name)),
    subtitle = "only showing routes within the 1.5 walk boundary",
    x = NULL,
    y = NULL,
    color = NULL,
    linewidth = "Potential student walkers") +
  theme(axis.text=element_blank(),
        axis.ticks=element_blank(),
        plot.caption = element_text(color = "grey")) +
  geom_sf(data = walk_boundary_poly,
          inherit.aes = FALSE,
          aes(color = paste0(1.5, " mile walking boundary")),
          fill = NA,
          linewidth = 1) +
  scale_color_manual(values = "blue", name = NULL) +
  new_scale_color() +
  geom_sf(data = bike_lts %>% filter(!is.na(student_use), student_use > 3),
          inherit.aes = FALSE,
          aes(linewidth = student_use),
          color = "mediumvioletred",
          fill = NA) +
  scale_linewidth_continuous(range = c(0, 3)) +
  annotation_raster(school_symbol,
                    # Position adjustments here using plot_box$max/min/range
                    ymin = as.double((WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE) %>% pull(geom))[[1]])[2] - 0.001,
                    ymax = as.double((WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE) %>% pull(geom))[[1]])[2] + 0.001,
                    xmin = as.double((WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE) %>% pull(geom))[[1]])[1] - 0.0015,
                    xmax = as.double((WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE) %>% pull(geom))[[1]])[1] + 0.0015) +
  geom_sf_label(data = WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE),
                inherit.aes = FALSE,
                mapping = aes(label = school_focus %>% pull(name)),
                nudge_y = 0.0015,
                label.size = 0.04,
                size = 2)

ggsave(file = paste0("figures/",
                     school_focus %>% pull(name),
                     " Routes.pdf"),
       title = paste0(school_focus %>% pull(name), " Walking Routes"),
       device = pdf,
       height = 8.5,
       width = 11,
       units = "in",
       create.dir = TRUE)

# generate map
ggmap(basemap) +
  labs(title = paste0("Walking routes for students at ",
                      school_focus %>% pull(name)),
       subtitle = "only showing routes within the walk boundary",
       x = NULL,
       y = NULL,
       color = NULL,
       linewidth = "Potential student walkers") +
  theme(axis.text=element_blank(),
        axis.ticks=element_blank(),
        plot.caption = element_text(color = "grey")) +
  geom_sf(data = walk_boundary_poly,
          inherit.aes = FALSE,
          aes(color = paste0(1.5, " mile walking boundary")),
          fill = NA,
          linewidth = 1) +
  scale_color_manual(values = "blue", name = NULL) +
  new_scale_color() +
  geom_sf(data = bike_lts %>% filter(!is.na(student_use), student_use > 0),
         inherit.aes = FALSE,
         aes(color = lts,
             linewidth = student_use)) +
  scale_color_manual(values = bike_lts_scale$color, name = "Bike Level of Traffic Stress") +
  scale_linewidth_continuous(range = c(0, 3)) +
  annotation_raster(school_symbol,
                    # Position adjustments here using plot_box$max/min/range
                    ymin = as.double((WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE) %>% pull(geom))[[1]])[2] - 0.001,
                    ymax = as.double((WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE) %>% pull(geom))[[1]])[2] + 0.001,
                    xmin = as.double((WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE) %>% pull(geom))[[1]])[1] - 0.0015,
                    xmax = as.double((WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE) %>% pull(geom))[[1]])[1] + 0.0015) +
  geom_sf_label(data = WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE),
                inherit.aes = FALSE,
                mapping = aes(label = school_focus %>% pull(name)),
                nudge_y = 0.0015,
                label.size = 0.04,
                size = 2)

ggsave(file = paste0("figures/",
                     school_focus %>% pull(name),
                     " Routes - Traffic Stress.pdf"),
       title = paste0(school_focus %>% pull(name), " Walking Routes - Traffic Stress"),
       device = pdf,
       height = 8.5,
       width = 11,
       units = "in",
       create.dir = TRUE)

ggmap(basemap) +
  labs(title = paste0("Student homes at ",
                      school_focus %>% pull(name)),
       x = NULL,
       y = NULL,
       color = NULL,
       fill = "How many students live there") +
  theme(axis.text=element_blank(),
        axis.ticks=element_blank(),
        plot.caption = element_text(color = "grey")) +
  geom_hex(data = addresses %>% extract(geometry, into = c('Lat', 'Lon'), '\\((.*),(.*)\\)', conv = T),
           aes(x = Lat,
               y = Lon),
           alpha = 0.7) +
  scale_fill_distiller(palette = "YlOrRd", direction = "reverse") +
  geom_sf(data = walk_boundary_poly,
          inherit.aes = FALSE,
          aes(color = paste0(1.5, " mile walking boundary")),
          fill = NA,
          linewidth = 1) +
  scale_color_manual(values = "blue", name = NULL) +
  new_scale_color() +
  annotation_raster(school_symbol,
                    # Position adjustments here using plot_box$max/min/range
                    ymin = as.double((WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE) %>% pull(geom))[[1]])[2] - 0.001,
                    ymax = as.double((WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE) %>% pull(geom))[[1]])[2] + 0.001,
                    xmin = as.double((WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE) %>% pull(geom))[[1]])[1] - 0.0015,
                    xmax = as.double((WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE) %>% pull(geom))[[1]])[1] + 0.0015) +
  geom_sf_label(data = WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE),
                inherit.aes = FALSE,
                mapping = aes(label = school_focus %>% pull(name)),
                nudge_y = 0.0015,
                label.size = 0.04,
                size = 2)

ggsave(file = paste0("figures/",
                     school_focus %>% pull(name),
                     " Addresses.pdf"),
       title = paste0(school_focus %>% pull(name), " Addresses"),
       device = pdf,
       height = 8.5,
       width = 11,
       units = "in",
       create.dir = TRUE)

Appendix

date()

## [1] "Thu Oct 31 16:38:18 2024"

sessionInfo()

## R version 4.1.2 (2021-11-01)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Pop!_OS 22.04 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] rsvg_2.6.1       ggnewscale_0.5.0 magick_2.8.5     smoothr_1.0.1   
##  [5] osrm_4.2.0       sf_1.0-18        ggmap_4.0.0      lubridate_1.9.3 
##  [9] forcats_1.0.0    stringr_1.5.1    dplyr_1.1.4      purrr_1.0.2     
## [13] readr_2.1.5      tidyr_1.3.1      tibble_3.2.1     ggplot2_3.5.1   
## [17] tidyverse_2.0.0 
## 
## loaded via a namespace (and not attached):
##  [1] httr_1.4.7            sass_0.4.9            bit64_4.5.2          
##  [4] vroom_1.6.5           jsonlite_1.8.9        bslib_0.8.0          
##  [7] highr_0.11            yaml_2.3.10           lattice_0.20-45      
## [10] pillar_1.9.0          glue_1.8.0            digest_0.6.37        
## [13] RColorBrewer_1.1-3    colorspace_2.1-1      htmltools_0.5.8.1    
## [16] plyr_1.8.9            pkgconfig_2.0.3       s2_1.1.7             
## [19] scales_1.3.0          terra_1.7-83          jpeg_0.1-10          
## [22] tzdb_0.4.0            timechange_0.3.0      proxy_0.4-27         
## [25] generics_0.1.3        farver_2.1.2          googlePolylines_0.8.5
## [28] cachem_1.1.0          withr_3.0.2           hexbin_1.28.4        
## [31] cli_3.6.3             magrittr_2.0.3        crayon_1.5.3         
## [34] evaluate_1.0.1        fansi_1.0.6           class_7.3-20         
## [37] tools_4.1.2           hms_1.1.3             RgoogleMaps_1.5.1    
## [40] mapiso_0.3.0          lifecycle_1.0.4       munsell_0.5.1        
## [43] isoband_0.2.7         compiler_4.1.2        jquerylib_0.1.4      
## [46] e1071_1.7-16          rlang_1.1.4           classInt_0.4-10      
## [49] units_0.8-5           grid_4.1.2            rstudioapi_0.17.1    
## [52] bitops_1.0-9          labeling_0.4.3        rmarkdown_2.28       
## [55] wk_0.9.4              gtable_0.3.6          codetools_0.2-18     
## [58] DBI_1.2.3             curl_5.2.3            RcppSimdJson_0.1.12  
## [61] R6_2.5.1              knitr_1.48            fastmap_1.2.0        
## [64] bit_4.5.0             utf8_1.2.4            KernSmooth_2.23-20   
## [67] stringi_1.8.4         parallel_4.1.2        Rcpp_1.0.13          
## [70] vctrs_0.6.5           png_0.1-8             tidyselect_1.2.1     
## [73] xfun_0.48