ben/route_analysis
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Compare commits

base: ben:f594e4cbd6841522bc601e5692cdd529136f3115
Branches Tags
ben:main
..
compare: ben:2d1d432710ec4b6b6abb84b77096a57c8cc22625
Branches Tags
ben:main

No commits in common. "f594e4cbd6841522bc601e5692cdd529136f3115" and "2d1d432710ec4b6b6abb84b77096a57c8cc22625" have entirely different histories.
f594e4cbd6 ... 2d1d432710

5 changed files with 0 additions and 550 deletions
Show Stats Expand all files Collapse all files

7
Makefile
View File

@ -2,7 +2,6 @@ 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
@ -14,12 +13,6 @@ 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
View File

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

57
R/functions.R
View File

@ -32,60 +32,3 @@ 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
View File

@ -1,290 +0,0 @@
---
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
View File

@ -1,192 +0,0 @@
---
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()
```