From d1b1fd0253cec41c965730a3b312eb083d844bb5 Mon Sep 17 00:00:00 2001
From: syounkin <syounkin@wisc.edu>
Date: Wed, 13 Nov 2024 14:00:07 -0600
Subject: [PATCH 01/11] Created route_to_school.Rmd

This script is a spin-off of cycle_route_analysis_brouter.Rmd.
---
 Makefile            |   5 +-
 R/data/.gitignore   |   4 +
 R/functions.R       |  47 +++++++++++
 route_to_school.Rmd | 198 ++++++++++++++++++++++++++++++++++++++++++++
 4 files changed, 253 insertions(+), 1 deletion(-)
 create mode 100644 R/data/.gitignore
 create mode 100644 route_to_school.Rmd

diff --git a/Makefile b/Makefile
index ee206b4..009db84 100644
--- a/Makefile
+++ b/Makefile
@@ -2,7 +2,7 @@ all: data containers cycle
 
 data: osrm-data brouter-data
 containers: osrm-container brouter-container
-cycle: cycle_brouter
+cycle: route_to_school
 
 walk: route_analysis.Rmd
 	R -e 'library("rmarkdown"); old_path <- Sys.getenv("PATH"); Sys.setenv(PATH = paste(old_path, "/usr/local/bin", sep = ":")); rmarkdown::render(knit_root_dir = "./", output_dir = "./html", input = "./route_analysis.Rmd", output_file = "./html/route_analysis.html")'
@@ -13,6 +13,9 @@ cycle_osrm: cycling_route_analysis.Rmd
 cycle_brouter: cycling_route_analysis_brouter.Rmd
 	R -e 'library("rmarkdown"); old_path <- Sys.getenv("PATH"); Sys.setenv(PATH = paste(old_path, "/usr/local/bin", sep = ":")); rmarkdown::render(knit_root_dir = "./", output_dir = "./html", input = "./cycling_route_analysis_brouter.Rmd", output_file = "./html/cycling_route_analysis.html")'
 
+route_to_school: route_to_school.Rmd
+	R -e 'library("rmarkdown"); old_path <- Sys.getenv("PATH"); Sys.setenv(PATH = paste(old_path, "/usr/local/bin", sep = ":")); rmarkdown::render(knit_root_dir = "./", output_dir = "./html", input = "./route_to_school.Rmd", output_file = "./html/route_to_school.html")'
+
 osrm-container: ./docker/osrm/docker-compose.yml
 	cd ./docker/osrm/; docker compose up -d
 
diff --git a/R/data/.gitignore b/R/data/.gitignore
new file mode 100644
index 0000000..5e7d273
--- /dev/null
+++ b/R/data/.gitignore
@@ -0,0 +1,4 @@
+# Ignore everything in this directory
+*
+# Except this file
+!.gitignore
diff --git a/R/functions.R b/R/functions.R
index c24604a..0f3f2de 100644
--- a/R/functions.R
+++ b/R/functions.R
@@ -32,3 +32,50 @@ getLTSForRoute <- function(i) {
 
   return(result)
 }
+
+
+routeChar <- function(route){
+
+text <- as.data.frame(route)$messages
+text <- gsub(x = text, pattern = "\\\"", replacement = "")
+text <- gsub(x = text, pattern = "\ ", replacement = "")
+text <- gsub(x = text, pattern = "\\[\\[", replacement = "")
+text <- gsub(x = text, pattern = "\\]\\]", replacement = "")
+foobar <- strsplit(text, split = "],[", fixed = TRUE)
+x <- lapply(foobar, function(x){strsplit(x, split = ",", fixed = TRUE)})
+xx <- unlist(x)
+m <- matrix(xx, ncol = 13, byrow = TRUE)
+df <- data.frame(m[-1,])
+names(df) <- m[1,]
+
+
+df2 <- within(df, {
+    Time <- as.numeric(Time)
+    stageTime <- diff(c(0,Time))
+    path <- grepl("highway=path", df$WayTags)
+    residential <- grepl("highway=residential", df$WayTags)
+    footway <- grepl("highway=footway", df$WayTags)
+    primary <- grepl("highway=primary", df$WayTags)
+    service <- grepl("highway=service", df$WayTags)
+    cycleway <- grepl("highway=cycleway", df$WayTags)
+    bike <- grepl("bicycle=designated", df$WayTags)
+})
+
+
+foo <- function(x){
+    ifelse(x$path, "path", ifelse(x$residential, "residential", ifelse(x$footway, "footway", ifelse(x$primary, "primary", ifelse(x$service, "service", ifelse(x$cycleway, "cycleway", "other"))))))
+}
+
+df2 <- cbind(df2, highway = foo(df2))
+df2 <- df2 %>% group_by(highway) %>% summarize(T = sum(stageTime))
+
+df2 <- df2 %>% filter(!is.na(highway))
+
+
+    if(!("cycleway" %in% df2$highway)){
+        return(0)
+    }else{
+        return(df2[df2$highway == "cycleway",]$T)
+    }
+
+}
diff --git a/route_to_school.Rmd b/route_to_school.Rmd
new file mode 100644
index 0000000..0758473
--- /dev/null
+++ b/route_to_school.Rmd
@@ -0,0 +1,198 @@
+---
+title: "East High Cycling Routes"
+output:
+  html_document:
+    toc: true
+    toc_depth: 5
+    toc_float:
+      collapsed: false
+      smooth_scroll: true
+editor_options:
+  chunk_output_type: console
+---
+
+# Input Data & Configuration
+
+## Libraries
+
+```{r libs, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
+date()
+rm(list=ls())
+library(tidyverse)
+library(ggmap)
+library(sf)
+library(osrm)
+library(smoothr)
+library(magick)
+library(ggnewscale)
+library(rsvg)
+library(httr)
+library(jsonlite)
+library(parallel)
+fig.height <- 6
+set.seed(1)
+makePlots <- TRUE
+source("./R/functions.R")
+```
+
+# External sources configurations
+
+## Open Source Routing Machine (OSRM)
+
+```{r osrm, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
+# Set url and profile of OSRM server
+options(osrm.server = "http://127.0.0.1:5001/")
+options(osrm.profile = "bike")
+```
+
+## Brouter options
+```{r brouter, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
+# Set url and profile of brouter server
+brouter_url <- "http://127.0.0.1:17777/brouter"
+brouter_profile <- "safety"
+```
+
+## Stadia Maps API Key
+
+```{r stadiamaps, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
+register_stadiamaps(key = substr(read_file(file = "api_keys/stadia_api_key"), 1, 36))
+```
+
+
+# Analysis
+
+## Create Bikeable Region Using OSRM
+
+```{r boundary, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+WI_schools <- st_transform(st_read(dsn = "data/Schools/Wisconsin_Public_Schools_-5986231931870160084.gpkg"), crs = 4326)
+WI_schools <- WI_schools %>% mutate(geom = SHAPE)
+
+school_focus <- data.frame(name = c("East High School"), NCES_CODE = c("550852000925"))
+#school_focus <- data.frame(name = c("IMAP"), NCES_CODE = c("550008203085"))
+school_location <- WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE)
+
+radius <- 3 # miles
+levels <- c(1)
+res <- 100
+threshold <- units::set_units(1, km^2)
+
+cycle_boundary_m <- radius*1609
+
+cycle_boundary_poly <- osrmIsodistance( loc = school_location, breaks = cycle_boundary_m, res = res )
+
+cycle_boundary_poly <- st_make_valid(cycle_boundary_poly)
+
+cycle_boundary_poly <- fill_holes(cycle_boundary_poly, threshold)
+
+cycle_boundary_poly <- st_transform(cycle_boundary_poly, crs = 4326)
+
+saveRDS(cycle_boundary_poly, "./R/data/cycle_boundary_poly.rds")
+```
+
+# Create Grid Over Bikeable Region
+
+```{r grid, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+cellsize <- 5e-3
+grid <- st_intersection(cycle_boundary_poly, st_make_grid(cycle_boundary_poly, cellsize = cellsize, what = "polygons", square = FALSE))
+```
+
+# Compute Routes from Cell Centroids to School with brouter
+
+```{r routes, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+grid_pts <- st_centroid(grid)
+grid_coods <- st_coordinates(grid_pts)
+school_focus_location <- WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE) %>% select(LAT, LON)
+
+routes <- list(NULL)
+for(i in 1:nrow(grid_coods) ) {
+  query <- paste0(
+    brouter_url,
+    "?lonlats=", grid_coods[i,1], ",",grid_coods[i,2], "|",
+    school_focus_location$LON, ",", school_focus_location$LAT,
+    "&profile=", brouter_profile,
+    "&alternativeidx=0&format=geojson"
+  )
+  response <- GET(query)
+  route_run <- st_read(content <- content(response, as = "text"), quiet = TRUE)
+  route_run[["student_number"]] <- i
+  routes[[i]] <- route_run
+}
+
+routes <- st_transform(bind_rows(routes), crs = 4326)
+```
+
+Notes:
+
+- What does `st_transform(bind_rows(routes), crs = 4326)` do?
+
+
+# Generate Map for Total Time
+
+## Set boundaries and get basemap
+```{r basemap, eval = makePlots, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
+bbox <- st_bbox(st_buffer(cycle_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]))
+
+basemap <- get_stadiamap(bbox = bbox, zoom = 15, maptype = "stamen_toner_lite")
+```
+
+## Route Characteristic Map
+
+```{r sandbox3, eval = makePlots, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+track.length.vec <- routes %>% pull(track.length)
+grid <- cbind(grid, track.length = as.numeric(track.length.vec)/1609)
+
+total.time.vec <- routes %>% pull(total.time)
+grid <- cbind(grid, total.time = as.numeric(total.time.vec)/60)
+
+total.energy.vec <- routes %>% pull(total.energy)
+grid <- cbind(grid, total.energy = as.numeric(total.energy.vec))
+
+gg1 <- ggmap(basemap) + geom_sf(data = grid, aes(fill = total.time), inherit.aes = FALSE)
+
+ggsave(gg1, filename = "./figures/route-characteristics.pdf", width = 11, height = 8, units = "in")
+
+gg1
+```
+## Routes Map
+
+```{r sandbox3b, eval = makePlots, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+gg2 <- ggmap(basemap) + geom_sf(data = routes, aes(color = "red"), inherit.aes = FALSE)
+ggsave(gg2, filename = "./figures/routes.pdf", width = 11, height = 8, units = "in")
+gg2
+```
+
+# Available Route Data
+
+## Investigatioin of Messages Data
+
+```{r sandbox4, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+x.vec <- c()
+for(j in 1:nrow(routes)){
+    foobar <- routeChar(routes[j, "messages"])
+    x.vec <- c(x.vec, foobar)
+}
+
+new.df <- cbind(grid, T.cycleway = x.vec)
+
+gg3 <- ggmap(basemap) + geom_sf(data = new.df, aes(fill= T.cycleway/60), inherit.aes = FALSE)
+ggsave(gg3, filename = "./figures/routes.pdf", width = 11, height = 8, units = "in")
+gg3
+```
+
+# Message Data?
+
+What information can we pull out of the messages data?
+
+```{r sandbox5, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+routes[1,"messages"]
+```
+
+
+```{r chunklast, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+date()
+sessionInfo()
+```

From 348608fc5a047037140ffd085965ac52cf0e8901 Mon Sep 17 00:00:00 2001
From: syounkin <syounkin@wisc.edu>
Date: Wed, 13 Nov 2024 14:24:15 -0600
Subject: [PATCH 02/11] Updates to route_to_school.Rmd

---
 route_to_school.Rmd | 22 ++++++++--------------
 1 file changed, 8 insertions(+), 14 deletions(-)

diff --git a/route_to_school.Rmd b/route_to_school.Rmd
index 0758473..e5819fa 100644
--- a/route_to_school.Rmd
+++ b/route_to_school.Rmd
@@ -92,7 +92,7 @@ saveRDS(cycle_boundary_poly, "./R/data/cycle_boundary_poly.rds")
 # Create Grid Over Bikeable Region
 
 ```{r grid, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
-cellsize <- 5e-3
+cellsize <- 2.5e-3
 grid <- st_intersection(cycle_boundary_poly, st_make_grid(cycle_boundary_poly, cellsize = cellsize, what = "polygons", square = FALSE))
 ```
 
@@ -146,12 +146,12 @@ track.length.vec <- routes %>% pull(track.length)
 grid <- cbind(grid, track.length = as.numeric(track.length.vec)/1609)
 
 total.time.vec <- routes %>% pull(total.time)
-grid <- cbind(grid, total.time = as.numeric(total.time.vec)/60)
+grid <- cbind(grid, total.time = as.numeric(total.time.vec))
 
 total.energy.vec <- routes %>% pull(total.energy)
 grid <- cbind(grid, total.energy = as.numeric(total.energy.vec))
 
-gg1 <- ggmap(basemap) + geom_sf(data = grid, aes(fill = total.time), inherit.aes = FALSE)
+gg1 <- ggmap(basemap) + geom_sf(data = subset(grid, track.length > 1), aes(fill = total.time), inherit.aes = FALSE)
 
 ggsave(gg1, filename = "./figures/route-characteristics.pdf", width = 11, height = 8, units = "in")
 
@@ -176,21 +176,15 @@ for(j in 1:nrow(routes)){
     x.vec <- c(x.vec, foobar)
 }
 
-new.df <- cbind(grid, T.cycleway = x.vec)
+grid <- cbind(grid, T.cycleway = x.vec)
+grid <- cbind( grid, pct.cycleway = 100*grid$T.cycleway/grid$total.time)
 
-gg3 <- ggmap(basemap) + geom_sf(data = new.df, aes(fill= T.cycleway/60), inherit.aes = FALSE)
-ggsave(gg3, filename = "./figures/routes.pdf", width = 11, height = 8, units = "in")
+gg3 <- ggmap(basemap) + geom_sf(data = subset(grid, track.length > 1), aes(fill= pct.cycleway), inherit.aes = FALSE)
+ggsave(gg3, filename = "./figures/cycleway.pdf", width = 11, height = 8, units = "in")
 gg3
 ```
 
-# Message Data?
-
-What information can we pull out of the messages data?
-
-```{r sandbox5, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
-routes[1,"messages"]
-```
-
+# Archive
 
 ```{r chunklast, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
 date()

From 099d94359ce070092501214d3bbc33e819b547b6 Mon Sep 17 00:00:00 2001
From: syounkin <syounkin@wisc.edu>
Date: Wed, 13 Nov 2024 17:15:51 -0600
Subject: [PATCH 03/11] Updates to route_to_school.Rmd

---
 route_to_school.Rmd | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/route_to_school.Rmd b/route_to_school.Rmd
index e5819fa..d1e8c2b 100644
--- a/route_to_school.Rmd
+++ b/route_to_school.Rmd
@@ -93,7 +93,8 @@ saveRDS(cycle_boundary_poly, "./R/data/cycle_boundary_poly.rds")
 
 ```{r grid, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
 cellsize <- 2.5e-3
-grid <- st_intersection(cycle_boundary_poly, st_make_grid(cycle_boundary_poly, cellsize = cellsize, what = "polygons", square = FALSE))
+grid <- st_make_grid(cycle_boundary_poly, cellsize = cellsize, what = "polygons", square = FALSE)
+grid <- st_intersection(cycle_boundary_poly, grid)
 ```
 
 # Compute Routes from Cell Centroids to School with brouter
@@ -136,10 +137,11 @@ bbox <- c(left = as.double(bbox[1]),
           right = as.double(bbox[3]),
           top = as.double(bbox[4]))
 
-basemap <- get_stadiamap(bbox = bbox, zoom = 15, maptype = "stamen_toner_lite")
+zoom.level <- 12
+basemap <- get_stadiamap(bbox = bbox, zoom = zoom.level, maptype = "stamen_toner_lite")
 ```
 
-## Route Characteristic Map
+## Total Trip Time Map
 
 ```{r sandbox3, eval = makePlots, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
 track.length.vec <- routes %>% pull(track.length)
@@ -157,6 +159,7 @@ ggsave(gg1, filename = "./figures/route-characteristics.pdf", width = 11, height
 
 gg1
 ```
+
 ## Routes Map
 
 ```{r sandbox3b, eval = makePlots, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
@@ -167,7 +170,7 @@ gg2
 
 # Available Route Data
 
-## Investigatioin of Messages Data
+## Compute Percent of Trip on Cycleway
 
 ```{r sandbox4, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
 x.vec <- c()

From afc6257b6179c09cfd6c3127b6fcbae0b4e7385b Mon Sep 17 00:00:00 2001
From: syounkin <syounkin@wisc.edu>
Date: Wed, 20 Nov 2024 09:51:08 -0600
Subject: [PATCH 04/11] Minor updates to route_to_school.Rmd

---
 route_to_school.Rmd | 21 +++++++++------------
 1 file changed, 9 insertions(+), 12 deletions(-)

diff --git a/route_to_school.Rmd b/route_to_school.Rmd
index d1e8c2b..c4db7e2 100644
--- a/route_to_school.Rmd
+++ b/route_to_school.Rmd
@@ -15,7 +15,7 @@ editor_options:
 
 ## Libraries
 
-```{r libs, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
+```{r libs, eval = TRUE, echo = FALSE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
 date()
 rm(list=ls())
 library(tidyverse)
@@ -31,7 +31,6 @@ library(jsonlite)
 library(parallel)
 fig.height <- 6
 set.seed(1)
-makePlots <- TRUE
 source("./R/functions.R")
 ```
 
@@ -40,14 +39,12 @@ source("./R/functions.R")
 ## Open Source Routing Machine (OSRM)
 
 ```{r osrm, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
-# Set url and profile of OSRM server
 options(osrm.server = "http://127.0.0.1:5001/")
 options(osrm.profile = "bike")
 ```
 
 ## Brouter options
 ```{r brouter, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
-# Set url and profile of brouter server
 brouter_url <- "http://127.0.0.1:17777/brouter"
 brouter_profile <- "safety"
 ```
@@ -71,7 +68,7 @@ school_focus <- data.frame(name = c("East High School"), NCES_CODE = c("55085200
 #school_focus <- data.frame(name = c("IMAP"), NCES_CODE = c("550008203085"))
 school_location <- WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE)
 
-radius <- 3 # miles
+radius <- 4 # miles
 levels <- c(1)
 res <- 100
 threshold <- units::set_units(1, km^2)
@@ -92,7 +89,7 @@ saveRDS(cycle_boundary_poly, "./R/data/cycle_boundary_poly.rds")
 # Create Grid Over Bikeable Region
 
 ```{r grid, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
-cellsize <- 2.5e-3
+cellsize <- 5e-3
 grid <- st_make_grid(cycle_boundary_poly, cellsize = cellsize, what = "polygons", square = FALSE)
 grid <- st_intersection(cycle_boundary_poly, grid)
 ```
@@ -130,7 +127,7 @@ Notes:
 # Generate Map for Total Time
 
 ## Set boundaries and get basemap
-```{r basemap, eval = makePlots, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
+```{r basemap, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
 bbox <- st_bbox(st_buffer(cycle_boundary_poly, dist = 500))
 bbox <- c(left = as.double(bbox[1]),
           bottom = as.double(bbox[2]),
@@ -143,7 +140,7 @@ basemap <- get_stadiamap(bbox = bbox, zoom = zoom.level, maptype = "stamen_toner
 
 ## Total Trip Time Map
 
-```{r sandbox3, eval = makePlots, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+```{r sandbox3, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
 track.length.vec <- routes %>% pull(track.length)
 grid <- cbind(grid, track.length = as.numeric(track.length.vec)/1609)
 
@@ -162,13 +159,13 @@ gg1
 
 ## Routes Map
 
-```{r sandbox3b, eval = makePlots, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+```{r sandbox3b, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
 gg2 <- ggmap(basemap) + geom_sf(data = routes, aes(color = "red"), inherit.aes = FALSE)
 ggsave(gg2, filename = "./figures/routes.pdf", width = 11, height = 8, units = "in")
 gg2
 ```
 
-# Available Route Data
+# Route Characteristics
 
 ## Compute Percent of Trip on Cycleway
 
@@ -180,9 +177,9 @@ for(j in 1:nrow(routes)){
 }
 
 grid <- cbind(grid, T.cycleway = x.vec)
-grid <- cbind( grid, pct.cycleway = 100*grid$T.cycleway/grid$total.time)
+grid <- cbind( grid, not.cycleway = (grid$total.time - grid$T.cycleway)/60)
 
-gg3 <- ggmap(basemap) + geom_sf(data = subset(grid, track.length > 1), aes(fill= pct.cycleway), inherit.aes = FALSE)
+gg3 <- ggmap(basemap) + geom_sf(data = grid, aes(fill= not.cycleway), inherit.aes = FALSE) + scale_fill_gradient(low = "yellow", high = "red", limits = c(0,17), na.value = NA)
 ggsave(gg3, filename = "./figures/cycleway.pdf", width = 11, height = 8, units = "in")
 gg3
 ```

From 7d7a9804004f4c5ec07f69a85698a8a9d90d7baf Mon Sep 17 00:00:00 2001
From: syounkin <syounkin@wisc.edu>
Date: Thu, 21 Nov 2024 11:32:27 -0600
Subject: [PATCH 05/11] Created WI-schools-cycle.Rmd

This script loops through WI schools and computes the cycling route from each grid cell to school. These routes and grids are saved to rds objects to be analyzed in a separate script.
---
 Makefile             |   5 +-
 WI-schools-cycle.Rmd | 131 +++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 135 insertions(+), 1 deletion(-)
 create mode 100644 WI-schools-cycle.Rmd

diff --git a/Makefile b/Makefile
index d21affa..e69b36e 100644
--- a/Makefile
+++ b/Makefile
@@ -2,7 +2,7 @@ all: data containers cycle
 
 data: osrm-data brouter-data
 containers: osrm-container brouter-container
-cycle: route_to_school
+cycle: WI-schools-cycle
 
 walk: route_analysis.Rmd
 	R -e 'library("rmarkdown"); old_path <- Sys.getenv("PATH"); Sys.setenv(PATH = paste(old_path, "/usr/local/bin", sep = ":")); rmarkdown::render(knit_root_dir = "./", output_dir = "./html", input = "./route_analysis.Rmd", output_file = "./html/route_analysis.html")'
@@ -16,6 +16,9 @@ cycle_brouter: cycling_route_analysis_brouter.Rmd
 route_to_school: route_to_school.Rmd
 	R -e 'library("rmarkdown"); old_path <- Sys.getenv("PATH"); Sys.setenv(PATH = paste(old_path, "/usr/local/bin", sep = ":")); rmarkdown::render(knit_root_dir = "./", output_dir = "./html", input = "./route_to_school.Rmd", output_file = "./html/route_to_school.html")'
 
+WI-schools-cycle: WI-schools-cycle.Rmd
+	R -e 'library("rmarkdown"); old_path <- Sys.getenv("PATH"); Sys.setenv(PATH = paste(old_path, "/usr/local/bin", sep = ":")); rmarkdown::render(knit_root_dir = "./", output_dir = "./html", input = "./WI-schools-cycle.Rmd", output_file = "./html/WI-schools-cycle.html")'
+
 osrm-container: ./docker/osrm/docker-compose.yml
 	cd ./docker/osrm/; docker compose up -d
 
diff --git a/WI-schools-cycle.Rmd b/WI-schools-cycle.Rmd
new file mode 100644
index 0000000..de4c0b7
--- /dev/null
+++ b/WI-schools-cycle.Rmd
@@ -0,0 +1,131 @@
+---
+title: "Wisconsin Cycling to School"
+output:
+  html_document:
+    toc: true
+    toc_depth: 5
+    toc_float:
+      collapsed: false
+      smooth_scroll: true
+editor_options:
+  chunk_output_type: console
+---
+
+# Input Data & Configuration
+
+## Libraries
+
+```{r libs, eval = TRUE, echo = FALSE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
+date()
+rm(list=ls())
+library(tidyverse)
+library(ggmap)
+library(sf)
+library(osrm)
+library(smoothr)
+library(magick)
+library(ggnewscale)
+library(rsvg)
+library(httr)
+library(jsonlite)
+library(parallel)
+fig.height <- 6
+set.seed(1)
+source("./R/functions.R")
+runLoop <- TRUE
+```
+## Configuration
+```{r config, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+options(osrm.server = "http://127.0.0.1:5001/")
+options(osrm.profile = "bike")
+brouter_url <- "http://127.0.0.1:17777/brouter"
+brouter_profile <- "safety"
+register_stadiamaps(key = substr(read_file(file = "api_keys/stadia_api_key"), 1, 36))
+WI_schools <- st_transform(st_read(dsn = "data/Schools/Wisconsin_Public_Schools_-5986231931870160084.gpkg"), crs = 4326)
+WI_schools <- WI_schools %>% mutate(geom = SHAPE)
+WI_schools <- subset(WI_schools, !is.na(LAT) & !is.na(LON) & GRADE_RANGE == "09-12")
+```
+## Loop through WI Schools
+```{r mainloop, eval = runLoop, echo = FALSE, results = "show", warning = FALSE, error = TRUE, message = TRUE}
+radius <- 3 # miles
+levels <- c(1)
+res <- 100
+threshold <- units::set_units(1, km^2)
+
+gridList <- list()
+routesList <- list()
+
+indexVec <- 1:nrow(WI_schools)
+
+for(j in indexVec){
+
+    school_location <- WI_schools[j,]
+
+    cycle_boundary_m <- radius*1609
+    cycle_boundary_poly <- osrmIsodistance( loc = school_location, breaks = cycle_boundary_m, res = res )
+    cycle_boundary_poly <- st_make_valid(cycle_boundary_poly)
+    cycle_boundary_poly <- fill_holes(cycle_boundary_poly, threshold)
+    cycle_boundary_poly <- st_transform(cycle_boundary_poly, crs = 4326)
+
+    cellsize <- 1e-2
+    grid <- st_make_grid(cycle_boundary_poly, cellsize = cellsize, what = "polygons", square = FALSE)
+    grid <- st_intersection(cycle_boundary_poly, grid)
+
+    grid <- st_make_valid(grid)
+
+    grid_pts <- st_centroid(grid)
+
+    grid_coods <- st_coordinates(grid_pts)
+    school_focus_location <- school_location %>% select(LAT, LON)
+
+    routes <- list(NULL)
+    for(i in 1:nrow(grid_coods) ) {
+        query <- paste0(
+            brouter_url,
+            "?lonlats=", grid_coods[i,1], ",",grid_coods[i,2], "|",
+            school_focus_location$LON, ",", school_focus_location$LAT,
+            "&profile=", brouter_profile,
+            "&alternativeidx=0&format=geojson"
+        )
+        response <- GET(query)
+
+        if( response$status_code == "200" ){
+            route_run <- st_read(content <- content(response, as = "text"), quiet = TRUE)
+            route_run[["student_number"]] <- i
+            routes[[i]] <- route_run
+        }else{
+            routes[[i]] <- NA
+        }
+    }
+
+    bad.cell <- which(is.na(routes))
+
+    if(length(bad.cell) > 0){
+        routes <- routes[-bad.cell]
+        grid <- grid[-bad.cell]
+    }
+
+
+    if(length(routes) > 0){
+        routes <- st_transform(bind_rows(routes), crs = 4326)
+    }else{
+        routes <- NA
+    }
+
+    gridList[[j]] <- grid
+    routesList[[j]] <- routes
+
+}
+
+saveRDS(gridList, "./R/data/gridList.rds")
+saveRDS(routesList, "./R/data/routesList.rds")
+```
+
+
+
+# Archive
+
+```{r chunklast, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+date()
+sessionInfo()
+```

From 0c9b3e070fb3545cda95a4b1a998be1f0f28e342 Mon Sep 17 00:00:00 2001
From: syounkin <syounkin@wisc.edu>
Date: Thu, 21 Nov 2024 14:48:22 -0600
Subject: [PATCH 06/11] Update to WI-schools-cycle.Rmd

Trying to get through all grade 9-12 schools in WI. Some schools have no routes, no coordinates, etc. Working my way through all the special cases.
---
 R/functions.R        | 70 +++++++++++++++++++++++++-------------------
 WI-schools-cycle.Rmd | 59 +++++++++++++++++++++++++++++++++----
 2 files changed, 93 insertions(+), 36 deletions(-)

diff --git a/R/functions.R b/R/functions.R
index 1e70ff8..b6dc7d6 100644
--- a/R/functions.R
+++ b/R/functions.R
@@ -36,40 +36,50 @@ getLTSForRoute <- function(i, route_table) {
 
 routeChar <- function(route){
 
-text <- as.data.frame(route)$messages
-text <- gsub(x = text, pattern = "\\\"", replacement = "")
-text <- gsub(x = text, pattern = "\ ", replacement = "")
-text <- gsub(x = text, pattern = "\\[\\[", replacement = "")
-text <- gsub(x = text, pattern = "\\]\\]", replacement = "")
-foobar <- strsplit(text, split = "],[", fixed = TRUE)
-x <- lapply(foobar, function(x){strsplit(x, split = ",", fixed = TRUE)})
-xx <- unlist(x)
-m <- matrix(xx, ncol = 13, byrow = TRUE)
-df <- data.frame(m[-1,])
-names(df) <- m[1,]
+    if(is.na(route$messages)){
+        return(NA)
+    }
+
+    text <- route$messages
+    text <- gsub(x = text, pattern = "\\\"", replacement = "")
+    text <- gsub(x = text, pattern = "\ ", replacement = "")
+    text <- gsub(x = text, pattern = "\\[\\[", replacement = "")
+    text <- gsub(x = text, pattern = "\\]\\]", replacement = "")
+    foobar <- strsplit(text, split = "],[", fixed = TRUE)
+    x <- lapply(foobar, function(x){strsplit(x, split = ",", fixed = TRUE)})
+    xx <- unlist(x)
+    m <- matrix(xx, ncol = 13, byrow = TRUE)
+    names.vec <- m[1,]
+
+    if(nrow(m) == 2){
+        df <- data.frame(t(m[-1,]))
+    }else{
+        df <- data.frame(m[-1,])
+    }
+
+    names(df) <- names.vec
+
+    df2 <- within(df, {
+        Time <- as.numeric(Time)
+        stageTime <- diff(c(0,Time))
+        path <- grepl("highway=path", df$WayTags)
+        residential <- grepl("highway=residential", df$WayTags)
+        footway <- grepl("highway=footway", df$WayTags)
+        primary <- grepl("highway=primary", df$WayTags)
+        service <- grepl("highway=service", df$WayTags)
+        cycleway <- grepl("highway=cycleway", df$WayTags)
+        bike <- grepl("bicycle=designated", df$WayTags)
+    })
 
 
-df2 <- within(df, {
-    Time <- as.numeric(Time)
-    stageTime <- diff(c(0,Time))
-    path <- grepl("highway=path", df$WayTags)
-    residential <- grepl("highway=residential", df$WayTags)
-    footway <- grepl("highway=footway", df$WayTags)
-    primary <- grepl("highway=primary", df$WayTags)
-    service <- grepl("highway=service", df$WayTags)
-    cycleway <- grepl("highway=cycleway", df$WayTags)
-    bike <- grepl("bicycle=designated", df$WayTags)
-})
+    foo <- function(x){
+        ifelse(x$path, "path", ifelse(x$residential, "residential", ifelse(x$footway, "footway", ifelse(x$primary, "primary", ifelse(x$service, "service", ifelse(x$cycleway, "cycleway", "other"))))))
+    }
 
+    df2 <- cbind(df2, highway = foo(df2))
+    df2 <- df2 %>% group_by(highway) %>% summarize(T = sum(stageTime))
 
-foo <- function(x){
-    ifelse(x$path, "path", ifelse(x$residential, "residential", ifelse(x$footway, "footway", ifelse(x$primary, "primary", ifelse(x$service, "service", ifelse(x$cycleway, "cycleway", "other"))))))
-}
-
-df2 <- cbind(df2, highway = foo(df2))
-df2 <- df2 %>% group_by(highway) %>% summarize(T = sum(stageTime))
-
-df2 <- df2 %>% filter(!is.na(highway))
+    df2 <- df2 %>% filter(!is.na(highway))
 
 
     if(!("cycleway" %in% df2$highway)){
diff --git a/WI-schools-cycle.Rmd b/WI-schools-cycle.Rmd
index de4c0b7..5e75b19 100644
--- a/WI-schools-cycle.Rmd
+++ b/WI-schools-cycle.Rmd
@@ -34,7 +34,9 @@ set.seed(1)
 source("./R/functions.R")
 runLoop <- TRUE
 ```
+
 ## Configuration
+
 ```{r config, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
 options(osrm.server = "http://127.0.0.1:5001/")
 options(osrm.profile = "bike")
@@ -43,9 +45,18 @@ brouter_profile <- "safety"
 register_stadiamaps(key = substr(read_file(file = "api_keys/stadia_api_key"), 1, 36))
 WI_schools <- st_transform(st_read(dsn = "data/Schools/Wisconsin_Public_Schools_-5986231931870160084.gpkg"), crs = 4326)
 WI_schools <- WI_schools %>% mutate(geom = SHAPE)
+```
+
+## Subset Schools
+
+```{r subsetSchools, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
 WI_schools <- subset(WI_schools, !is.na(LAT) & !is.na(LON) & GRADE_RANGE == "09-12")
 ```
+
+Non-virtual, grades 9-12.
+
 ## Loop through WI Schools
+
 ```{r mainloop, eval = runLoop, echo = FALSE, results = "show", warning = FALSE, error = TRUE, message = TRUE}
 radius <- 3 # miles
 levels <- c(1)
@@ -56,6 +67,9 @@ gridList <- list()
 routesList <- list()
 
 indexVec <- 1:nrow(WI_schools)
+#indexVec <- 1:50
+jj <- 1;
+bad.school.vec <- c()
 
 for(j in indexVec){
 
@@ -91,7 +105,6 @@ for(j in indexVec){
 
         if( response$status_code == "200" ){
             route_run <- st_read(content <- content(response, as = "text"), quiet = TRUE)
-            route_run[["student_number"]] <- i
             routes[[i]] <- route_run
         }else{
             routes[[i]] <- NA
@@ -102,26 +115,60 @@ for(j in indexVec){
 
     if(length(bad.cell) > 0){
         routes <- routes[-bad.cell]
-        grid <- grid[-bad.cell]
+        grid <- grid[-bad.cell,]
     }
 
-
     if(length(routes) > 0){
         routes <- st_transform(bind_rows(routes), crs = 4326)
+        gridList[[jj]] <- grid
+        routesList[[jj]] <- routes
+        jj <- jj + 1
     }else{
         routes <- NA
+        bad.school.vec <- c(j, bad.school.vec)
+        cat( WI_schools$SCHOOL[j], "has zero routes to school and has been removed from analysis.\n")
     }
 
-    gridList[[j]] <- grid
-    routesList[[j]] <- routes
-
 }
 
+WI_schools <- WI_schools[-bad.school.vec,]
+
+saveRDS(WI_schools, "./R/data/WI_schools.rds")
 saveRDS(gridList, "./R/data/gridList.rds")
 saveRDS(routesList, "./R/data/routesList.rds")
 ```
 
+```{r showLists, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+routesList <- readRDS(file = "./R/data/routesList.rds")
+gridList <- readRDS(file = "./R/data/gridList.rds")
 
+not.cycleway.vec <- c()
+
+for(j in 1:length(gridList)){
+
+    grid <- gridList[[j]]
+    routes <- routesList[[j]]
+
+    total.time.vec <- routes %>% pull(total.time)
+    grid <- cbind(grid, total.time = as.numeric(total.time.vec))
+
+    x.vec <- c()
+    for( i in 1:nrow(grid) ){
+        route <- routes[i,"messages"] # Grid cell i to school j
+        x <- routeChar(route)
+        x.vec <- c(x.vec, x)
+    }
+
+    grid <- cbind(grid, T.cycleway = x.vec)
+    grid <- cbind( grid, not.cycleway = (grid$total.time - grid$T.cycleway)/60)
+    gridList[[j]] <- grid
+
+    not.cycleway.vec <- c(not.cycleway.vec, median(grid$not.cycleway))
+
+}
+
+ggplot(data.frame(not.cycleway = not.cycleway.vec), aes(not.cycleway)) + geom_histogram(fill = "grey", color = "black") + theme_bw()
+```
 
 # Archive
 

From 5c12c14f8dd16dbf8f1e80d94be65e600a6911ae Mon Sep 17 00:00:00 2001
From: syounkin <syounkin@wisc.edu>
Date: Thu, 21 Nov 2024 15:13:37 -0600
Subject: [PATCH 07/11] Bug fix in WI-schools-cycle.Rmd

---
 WI-schools-cycle.Rmd | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/WI-schools-cycle.Rmd b/WI-schools-cycle.Rmd
index 5e75b19..bd8f410 100644
--- a/WI-schools-cycle.Rmd
+++ b/WI-schools-cycle.Rmd
@@ -67,7 +67,7 @@ gridList <- list()
 routesList <- list()
 
 indexVec <- 1:nrow(WI_schools)
-#indexVec <- 1:50
+
 jj <- 1;
 bad.school.vec <- c()
 
@@ -131,7 +131,9 @@ for(j in indexVec){
 
 }
 
-WI_schools <- WI_schools[-bad.school.vec,]
+if(length(bad.school.vec) > 0){
+    WI_schools <- WI_schools[-bad.school.vec,]
+}
 
 saveRDS(WI_schools, "./R/data/WI_schools.rds")
 saveRDS(gridList, "./R/data/gridList.rds")

From a1be8f5f25074b6bcc1b7924bf97539a24c94694 Mon Sep 17 00:00:00 2001
From: syounkin <syounkin@wisc.edu>
Date: Thu, 21 Nov 2024 17:10:46 -0600
Subject: [PATCH 08/11] Update to WI-schools-cycle.Rmd

---
 WI-schools-cycle.Rmd | 57 ++++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 55 insertions(+), 2 deletions(-)

diff --git a/WI-schools-cycle.Rmd b/WI-schools-cycle.Rmd
index bd8f410..3f8b5aa 100644
--- a/WI-schools-cycle.Rmd
+++ b/WI-schools-cycle.Rmd
@@ -32,7 +32,7 @@ library(parallel)
 fig.height <- 6
 set.seed(1)
 source("./R/functions.R")
-runLoop <- TRUE
+runLoop <- FALSE
 ```
 
 ## Configuration
@@ -140,9 +140,12 @@ saveRDS(gridList, "./R/data/gridList.rds")
 saveRDS(routesList, "./R/data/routesList.rds")
 ```
 
-```{r showLists, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+## Read List Data
+
+```{r readLists, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
 routesList <- readRDS(file = "./R/data/routesList.rds")
 gridList <- readRDS(file = "./R/data/gridList.rds")
+WI_schools <- readRDS(file = "./R/data/WI_schools.rds")
 
 not.cycleway.vec <- c()
 
@@ -168,8 +171,58 @@ for(j in 1:length(gridList)){
     not.cycleway.vec <- c(not.cycleway.vec, median(grid$not.cycleway))
 
 }
+```
 
+## Plot List Data
+
+### Best & Worst Schools
+
+```{r plots, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
 ggplot(data.frame(not.cycleway = not.cycleway.vec), aes(not.cycleway)) + geom_histogram(fill = "grey", color = "black") + theme_bw()
+
+register_stadiamaps(key = substr(read_file(file = "api_keys/stadia_api_key"), 1, 36))
+
+zoom.level <- 15
+k <- 306 #247
+
+bbox <- st_bbox(st_buffer(gridList[[k]], dist = 500))
+bbox <- c(left = as.double(bbox[1]),
+          bottom = as.double(bbox[2]),
+          right = as.double(bbox[3]),
+          top = as.double(bbox[4]))
+basemap <- get_stadiamap(bbox = bbox, zoom = zoom.level, maptype = "stamen_toner_lite")
+
+ggmap(basemap) + geom_sf(data = gridList[[k]], aes(fill= not.cycleway), inherit.aes = FALSE) + scale_fill_gradient(low = "yellow", high = "red", limits = c(0,17), na.value = NA)
+
+k <- 247 #306
+
+bbox <- st_bbox(st_buffer(gridList[[k]], dist = 500))
+bbox <- c(left = as.double(bbox[1]),
+          bottom = as.double(bbox[2]),
+          right = as.double(bbox[3]),
+          top = as.double(bbox[4]))
+basemap <- get_stadiamap(bbox = bbox, zoom = zoom.level, maptype = "stamen_toner_lite")
+
+ggmap(basemap) + geom_sf(data = gridList[[k]], aes(fill= not.cycleway), inherit.aes = FALSE) + scale_fill_gradient(low = "yellow", high = "red", limits = c(0,17), na.value = NA)
+```
+
+### Statewide Map
+
+```{r plots2, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+
+D <- cbind(WI_schools, not.cycleway.vec)
+
+zoom.level <- 8
+
+bbox <- st_bbox(st_buffer(D, dist = 10e3))
+bbox <- c(left = as.double(bbox[1]),
+          bottom = as.double(bbox[2]),
+          right = as.double(bbox[3]),
+          top = as.double(bbox[4]))
+basemap <- get_stadiamap(bbox = bbox, zoom = zoom.level, maptype = "stamen_toner_lite")
+
+ggmap(basemap) + geom_sf(data = D, aes(size = 2, color = not.cycleway.vec), inherit.aes = FALSE) + scale_color_gradient(low = "yellow", high = "red", na.value = NA) # , limits = c(0,17)
+
 ```
 
 # Archive

From a0b7d2fdf90b89fc5ff0ff863e572292ad16f578 Mon Sep 17 00:00:00 2001
From: syounkin <syounkin@wisc.edu>
Date: Fri, 22 Nov 2024 10:24:46 -0600
Subject: [PATCH 09/11] Added Table to WI-schools-cycle.Rmd

Used package reactable to create a table of non.cycleway values. North High in Eau Claire had the shortest median duration of 4.9 min. and River Ridge High the longest at 17.9 min.
---
 WI-schools-cycle.Rmd | 77 +++++++++++++++++++++++++++++---------------
 1 file changed, 51 insertions(+), 26 deletions(-)

diff --git a/WI-schools-cycle.Rmd b/WI-schools-cycle.Rmd
index 3f8b5aa..b2c76aa 100644
--- a/WI-schools-cycle.Rmd
+++ b/WI-schools-cycle.Rmd
@@ -22,13 +22,14 @@ library(tidyverse)
 library(ggmap)
 library(sf)
 library(osrm)
-library(smoothr)
-library(magick)
-library(ggnewscale)
-library(rsvg)
-library(httr)
-library(jsonlite)
-library(parallel)
+## library(smoothr)
+## library(magick)
+## library(ggnewscale)
+## library(rsvg)
+## library(httr)
+## library(jsonlite)
+## library(parallel)
+library(reactable)
 fig.height <- 6
 set.seed(1)
 source("./R/functions.R")
@@ -37,13 +38,13 @@ runLoop <- FALSE
 
 ## Configuration
 
-```{r config, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+```{r config, eval = TRUE, echo = FALSE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
 options(osrm.server = "http://127.0.0.1:5001/")
 options(osrm.profile = "bike")
 brouter_url <- "http://127.0.0.1:17777/brouter"
 brouter_profile <- "safety"
 register_stadiamaps(key = substr(read_file(file = "api_keys/stadia_api_key"), 1, 36))
-WI_schools <- st_transform(st_read(dsn = "data/Schools/Wisconsin_Public_Schools_-5986231931870160084.gpkg"), crs = 4326)
+WI_schools <- st_transform(st_read(dsn = "data/Schools/Wisconsin_Public_Schools_-5986231931870160084.gpkg", quiet = TRUE), crs = 4326)
 WI_schools <- WI_schools %>% mutate(geom = SHAPE)
 ```
 
@@ -142,7 +143,7 @@ saveRDS(routesList, "./R/data/routesList.rds")
 
 ## Read List Data
 
-```{r readLists, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+```{r readLists, eval = TRUE, echo = FALSE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
 routesList <- readRDS(file = "./R/data/routesList.rds")
 gridList <- readRDS(file = "./R/data/gridList.rds")
 WI_schools <- readRDS(file = "./R/data/WI_schools.rds")
@@ -175,26 +176,22 @@ for(j in 1:length(gridList)){
 
 ## Plot List Data
 
-### Best & Worst Schools
+### Shortest and Longest Median Non-Cycleway Duration
 
-```{r plots, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+First we take a look at the schools with the shortest and longest
+median time on cycleway.
+
+#### Longest
+
+
+
+```{r worst, eval = TRUE, echo = FALSE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
 ggplot(data.frame(not.cycleway = not.cycleway.vec), aes(not.cycleway)) + geom_histogram(fill = "grey", color = "black") + theme_bw()
 
 register_stadiamaps(key = substr(read_file(file = "api_keys/stadia_api_key"), 1, 36))
-
+k <- 306
 zoom.level <- 15
-k <- 306 #247
 
-bbox <- st_bbox(st_buffer(gridList[[k]], dist = 500))
-bbox <- c(left = as.double(bbox[1]),
-          bottom = as.double(bbox[2]),
-          right = as.double(bbox[3]),
-          top = as.double(bbox[4]))
-basemap <- get_stadiamap(bbox = bbox, zoom = zoom.level, maptype = "stamen_toner_lite")
-
-ggmap(basemap) + geom_sf(data = gridList[[k]], aes(fill= not.cycleway), inherit.aes = FALSE) + scale_fill_gradient(low = "yellow", high = "red", limits = c(0,17), na.value = NA)
-
-k <- 247 #306
 
 bbox <- st_bbox(st_buffer(gridList[[k]], dist = 500))
 bbox <- c(left = as.double(bbox[1]),
@@ -206,9 +203,28 @@ basemap <- get_stadiamap(bbox = bbox, zoom = zoom.level, maptype = "stamen_toner
 ggmap(basemap) + geom_sf(data = gridList[[k]], aes(fill= not.cycleway), inherit.aes = FALSE) + scale_fill_gradient(low = "yellow", high = "red", limits = c(0,17), na.value = NA)
 ```
 
+The longest is `r WI_schools[k,] |> pull(SCHOOL)`.
+
+#### Shortest
+
+```{r best, eval = TRUE, echo = FALSE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+k <- 247
+
+bbox <- st_bbox(st_buffer(gridList[[k]], dist = 500))
+bbox <- c(left = as.double(bbox[1]),
+          bottom = as.double(bbox[2]),
+          right = as.double(bbox[3]),
+          top = as.double(bbox[4]))
+basemap <- get_stadiamap(bbox = bbox, zoom = zoom.level, maptype = "stamen_toner_lite")
+
+ggmap(basemap) + geom_sf(data = gridList[[k]], aes(fill= not.cycleway), inherit.aes = FALSE) + scale_fill_gradient(low = "yellow", high = "red", limits = c(0,17), na.value = NA)
+```
+
+The shortest is `r WI_schools[k,] |> pull(SCHOOL)`.
+
 ### Statewide Map
 
-```{r plots2, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+```{r plots2, eval = TRUE, echo = FALSE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
 
 D <- cbind(WI_schools, not.cycleway.vec)
 
@@ -225,9 +241,18 @@ ggmap(basemap) + geom_sf(data = D, aes(size = 2, color = not.cycleway.vec), inhe
 
 ```
 
+### Statewide Table
+
+```{r table, eval = TRUE, echo = FALSE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+D_table <- as.data.frame(cbind(WI_schools, not.cycleway.vec))
+D <- select(D_table, SCHOOL, COUNTY, NONCYCLEWAY = not.cycleway.vec)
+D <- D |> mutate(NONCYCLEWAY = round(NONCYCLEWAY,1))
+reactable(D)
+```
+
 # Archive
 
-```{r chunklast, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+```{r chunklast, eval = TRUE, echo = FALSE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
 date()
 sessionInfo()
 ```

From 15656c9ca0c7d038ae765973db379dec8013f8cc Mon Sep 17 00:00:00 2001
From: syounkin <syounkin@wisc.edu>
Date: Fri, 22 Nov 2024 11:01:29 -0600
Subject: [PATCH 10/11] Moved route_to_school.Rmd to archive

---
 route_to_school.Rmd => archive/route_to_school.Rmd | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 rename route_to_school.Rmd => archive/route_to_school.Rmd (100%)

diff --git a/route_to_school.Rmd b/archive/route_to_school.Rmd
similarity index 100%
rename from route_to_school.Rmd
rename to archive/route_to_school.Rmd

From d0eecdc2c5742e430aeb4f75352d63cb86decb88 Mon Sep 17 00:00:00 2001
From: syounkin <syounkin@wisc.edu>
Date: Fri, 22 Nov 2024 11:08:43 -0600
Subject: [PATCH 11/11] Added descriptive text to WI-schools-cycle.Rmd

---
 Makefile             |  3 +-
 WI-schools-cycle.Rmd | 74 +++++++++++++++++++++++++++++++-------------
 2 files changed, 55 insertions(+), 22 deletions(-)

diff --git a/Makefile b/Makefile
index e69b36e..eb42337 100644
--- a/Makefile
+++ b/Makefile
@@ -2,7 +2,8 @@ all: data containers cycle
 
 data: osrm-data brouter-data
 containers: osrm-container brouter-container
-cycle: WI-schools-cycle
+WI-cycle: WI-schools-cycle
+cycle: cycle_brouter
 
 walk: route_analysis.Rmd
 	R -e 'library("rmarkdown"); old_path <- Sys.getenv("PATH"); Sys.setenv(PATH = paste(old_path, "/usr/local/bin", sep = ":")); rmarkdown::render(knit_root_dir = "./", output_dir = "./html", input = "./route_analysis.Rmd", output_file = "./html/route_analysis.html")'
diff --git a/WI-schools-cycle.Rmd b/WI-schools-cycle.Rmd
index b2c76aa..b7971a5 100644
--- a/WI-schools-cycle.Rmd
+++ b/WI-schools-cycle.Rmd
@@ -11,24 +11,20 @@ editor_options:
   chunk_output_type: console
 ---
 
+```{r preCode, eval = TRUE, echo = FALSE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
+date()
+```
+
 # Input Data & Configuration
 
 ## Libraries
 
 ```{r libs, eval = TRUE, echo = FALSE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
-date()
 rm(list=ls())
 library(tidyverse)
 library(ggmap)
 library(sf)
 library(osrm)
-## library(smoothr)
-## library(magick)
-## library(ggnewscale)
-## library(rsvg)
-## library(httr)
-## library(jsonlite)
-## library(parallel)
 library(reactable)
 fig.height <- 6
 set.seed(1)
@@ -38,6 +34,8 @@ runLoop <- FALSE
 
 ## Configuration
 
+Set configuration parameters for OSRM, brouter, and stadiamaps.
+
 ```{r config, eval = TRUE, echo = FALSE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
 options(osrm.server = "http://127.0.0.1:5001/")
 options(osrm.profile = "bike")
@@ -48,16 +46,32 @@ WI_schools <- st_transform(st_read(dsn = "data/Schools/Wisconsin_Public_Schools_
 WI_schools <- WI_schools %>% mutate(geom = SHAPE)
 ```
 
+# Analysis
+
+We focus on the statistic *non-cycleway duration* in this analysis. It
+is computed as the duration (in minutes) of the bike trip to school
+(brouter, safety) for each grid cell in the school's bikeable
+area. A bikeable area is defined as the region within 3 miles of
+school by bike (OSRM).
+
 ## Subset Schools
 
 ```{r subsetSchools, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
 WI_schools <- subset(WI_schools, !is.na(LAT) & !is.na(LON) & GRADE_RANGE == "09-12")
 ```
 
-Non-virtual, grades 9-12.
+We keep only schools with coordinates (non-virtual) and, for
+simplicity and efficiency of the initial analysis, we keep only
+schools with grades 9-12.
 
 ## Loop through WI Schools
 
+For each school we compute the grid and the routes sf objects and save
+them as lists as R data files, _gridList.rds_ and
+_routesList.rds_. These will then be analyzed downstream and this loop
+need not be run again. It took around 40 minutes to run. The code here
+is suppressed because it is long and ugly.
+
 ```{r mainloop, eval = runLoop, echo = FALSE, results = "show", warning = FALSE, error = TRUE, message = TRUE}
 radius <- 3 # miles
 levels <- c(1)
@@ -176,18 +190,31 @@ for(j in 1:length(gridList)){
 
 ## Plot List Data
 
-### Shortest and Longest Median Non-Cycleway Duration
+### Median Non-Cycleway Duration
 
-First we take a look at the schools with the shortest and longest
-median time on cycleway.
+#### Histogram
+
+First we investigate the distribution of median non-cycleway duration
+across school. Recall that we are considering now only schools grade
+9-12.
+
+```{r hist, eval = TRUE, echo = FALSE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
+ggplot(data.frame(not.cycleway = not.cycleway.vec), aes(not.cycleway)) + geom_histogram(fill = "orange", color = "black") + theme_bw()
+```
+
+Next, we take a look at the schools with the shortest and longest
+median time on cycleway. Note that the analysis is peformed across a
+gridded area and not with respect to where students live. The median
+non-cycleway duration is computed across grid cells, not students.
+
+Note too that this statistics was computed by parsing the *messages*
+field of the route returned by brouter. I am not sure if there is a
+better way to do this. Within the messages field there is information
+on highway type, surface, etc for each segment of the route.
 
 #### Longest
 
-
-
-```{r worst, eval = TRUE, echo = FALSE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
-ggplot(data.frame(not.cycleway = not.cycleway.vec), aes(not.cycleway)) + geom_histogram(fill = "grey", color = "black") + theme_bw()
-
+```{r worst, eval = TRUE, echo = FALSE, results = "show", warning = TRUE, error = TRUE, message = FALSE}
 register_stadiamaps(key = substr(read_file(file = "api_keys/stadia_api_key"), 1, 36))
 k <- 306
 zoom.level <- 15
@@ -207,7 +234,7 @@ The longest is `r WI_schools[k,] |> pull(SCHOOL)`.
 
 #### Shortest
 
-```{r best, eval = TRUE, echo = FALSE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+```{r best, eval = TRUE, echo = FALSE, results = "show", warning = TRUE, error = TRUE, message = FALSE}
 k <- 247
 
 bbox <- st_bbox(st_buffer(gridList[[k]], dist = 500))
@@ -224,7 +251,7 @@ The shortest is `r WI_schools[k,] |> pull(SCHOOL)`.
 
 ### Statewide Map
 
-```{r plots2, eval = TRUE, echo = FALSE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+```{r plots2, eval = TRUE, echo = FALSE, results = "show", warning = TRUE, error = TRUE, message = FALSE}
 
 D <- cbind(WI_schools, not.cycleway.vec)
 
@@ -243,16 +270,21 @@ ggmap(basemap) + geom_sf(data = D, aes(size = 2, color = not.cycleway.vec), inhe
 
 ### Statewide Table
 
+The values shown above can be seen below in this clickable table.
+
 ```{r table, eval = TRUE, echo = FALSE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
 D_table <- as.data.frame(cbind(WI_schools, not.cycleway.vec))
 D <- select(D_table, SCHOOL, COUNTY, NONCYCLEWAY = not.cycleway.vec)
 D <- D |> mutate(NONCYCLEWAY = round(NONCYCLEWAY,1))
 reactable(D)
 ```
+```{r date, eval = TRUE, echo = FALSE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+date()
+```
 
 # Archive
 
-```{r chunklast, eval = TRUE, echo = FALSE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
-date()
+## Session Info
+```{r sessionInfo, eval = TRUE, echo = FALSE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
 sessionInfo()
 ```