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

Merge pull request #18 from syounkin/sgy

Browse source 
Created New Script route_to_school.Rmd
This commit is contained in:
bvarick 2024-11-21 17:12:56 -06:00 • committed by GitHub
commit bfa8e3880f
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
5 changed files with 493 additions and 1 deletions
Download patch file Download diff file Expand all files Collapse all files

8
Makefile
View file

@ -2,7 +2,7 @@ all: data containers cycle
data: osrm-data brouter-data
containers: osrm-container brouter-container
cycle: cycle_brouter
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")'
@ -13,6 +13,12 @@ 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")'
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

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

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

57
R/functions.R
View file

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

233
WI-schools-cycle.Rmd Normal file
View file

@ -0,0 +1,233 @@
---
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 <- FALSE
```
## 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)
```
## 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)
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 = 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()
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
### 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
```{r chunklast, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
date()
sessionInfo()
```

192
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()
```