Compare commits

..

14 Commits

Author SHA1 Message Date
bvarick
f594e4cbd6
Merge pull request #19 from syounkin/sgy
Updates to WI-schools-cycle.Rmd
2024-11-27 08:31:44 -06:00
syounkin
d0eecdc2c5 Added descriptive text to WI-schools-cycle.Rmd 2024-11-22 11:08:43 -06:00
syounkin
15656c9ca0 Moved route_to_school.Rmd to archive 2024-11-22 11:01:29 -06:00
syounkin
a0b7d2fdf9 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.
2024-11-22 10:24:46 -06:00
bvarick
bfa8e3880f
Merge pull request #18 from syounkin/sgy
Created New Script route_to_school.Rmd
2024-11-21 17:12:56 -06:00
syounkin
a1be8f5f25 Update to WI-schools-cycle.Rmd 2024-11-21 17:10:46 -06:00
syounkin
5c12c14f8d Bug fix in WI-schools-cycle.Rmd 2024-11-21 15:13:37 -06:00
syounkin
0c9b3e070f 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.
2024-11-21 14:48:22 -06:00
syounkin
7d7a980400 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.
2024-11-21 11:32:27 -06:00
syounkin
912e12de2d Merge branch 'main' into sgy 2024-11-20 09:52:58 -06:00
syounkin
afc6257b61 Minor updates to route_to_school.Rmd 2024-11-20 09:51:08 -06:00
syounkin
099d94359c Updates to route_to_school.Rmd 2024-11-13 17:15:51 -06:00
syounkin
348608fc5a Updates to route_to_school.Rmd 2024-11-13 14:24:15 -06:00
syounkin
d1b1fd0253 Created route_to_school.Rmd
This script is a spin-off of cycle_route_analysis_brouter.Rmd.
2024-11-13 14:00:07 -06:00
5 changed files with 550 additions and 0 deletions

View File

@ -2,6 +2,7 @@ all: data containers cycle
data: osrm-data brouter-data data: osrm-data brouter-data
containers: osrm-container brouter-container containers: osrm-container brouter-container
WI-cycle: WI-schools-cycle
cycle: cycle_brouter cycle: cycle_brouter
walk: route_analysis.Rmd walk: route_analysis.Rmd
@ -13,6 +14,12 @@ cycle_osrm: cycling_route_analysis.Rmd
cycle_brouter: cycling_route_analysis_brouter.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")' 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")'
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 osrm-container: ./docker/osrm/docker-compose.yml
cd ./docker/osrm/; docker compose up -d cd ./docker/osrm/; docker compose up -d

4
R/data/.gitignore vendored Normal file
View File

@ -0,0 +1,4 @@
# Ignore everything in this directory
*
# Except this file
!.gitignore

View File

@ -32,3 +32,60 @@ getLTSForRoute <- function(i, route_table) {
return(result) return(result)
} }
routeChar <- function(route){
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)
})
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)
}
}

290
WI-schools-cycle.Rmd Normal file
View File

@ -0,0 +1,290 @@
---
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
---
```{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}
rm(list=ls())
library(tidyverse)
library(ggmap)
library(sf)
library(osrm)
library(reactable)
fig.height <- 6
set.seed(1)
source("./R/functions.R")
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")
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", quiet = TRUE), crs = 4326)
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")
```
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)
res <- 100
threshold <- units::set_units(1, km^2)
gridList <- list()
routesList <- list()
indexVec <- 1:nrow(WI_schools)
jj <- 1;
bad.school.vec <- c()
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)
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)
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")
}
}
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")
saveRDS(routesList, "./R/data/routesList.rds")
```
## Read List Data
```{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")
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))
}
```
## Plot List Data
### Median Non-Cycleway Duration
#### 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 = FALSE}
register_stadiamaps(key = substr(read_file(file = "api_keys/stadia_api_key"), 1, 36))
k <- 306
zoom.level <- 15
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 longest is `r WI_schools[k,] |> pull(SCHOOL)`.
#### Shortest
```{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))
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 = FALSE, results = "show", warning = TRUE, error = TRUE, message = FALSE}
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)
```
### 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
## Session Info
```{r sessionInfo, eval = TRUE, echo = FALSE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
sessionInfo()
```

192
archive/route_to_school.Rmd Normal file
View File

@ -0,0 +1,192 @@
---
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 = 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")
```
# External sources configurations
## Open Source Routing Machine (OSRM)
```{r osrm, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
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}
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 <- 4 # 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_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
```{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 = 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]),
right = as.double(bbox[3]),
top = as.double(bbox[4]))
zoom.level <- 12
basemap <- get_stadiamap(bbox = bbox, zoom = zoom.level, maptype = "stamen_toner_lite")
```
## Total Trip Time Map
```{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)
total.time.vec <- routes %>% pull(total.time)
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 = 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")
gg1
```
## Routes Map
```{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
```
# Route Characteristics
## Compute Percent of Trip on Cycleway
```{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)
}
grid <- cbind(grid, T.cycleway = x.vec)
grid <- cbind( grid, not.cycleway = (grid$total.time - grid$T.cycleway)/60)
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
```
# Archive
```{r chunklast, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
date()
sessionInfo()
```