From d1b1fd0253cec41c965730a3b312eb083d844bb5 Mon Sep 17 00:00:00 2001 From: syounkin Date: Wed, 13 Nov 2024 14:00:07 -0600 Subject: [PATCH 1/8] 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 Date: Wed, 13 Nov 2024 14:24:15 -0600 Subject: [PATCH 2/8] 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 Date: Wed, 13 Nov 2024 17:15:51 -0600 Subject: [PATCH 3/8] 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 Date: Wed, 20 Nov 2024 09:51:08 -0600 Subject: [PATCH 4/8] 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 Date: Thu, 21 Nov 2024 11:32:27 -0600 Subject: [PATCH 5/8] 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 Date: Thu, 21 Nov 2024 14:48:22 -0600 Subject: [PATCH 6/8] 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 Date: Thu, 21 Nov 2024 15:13:37 -0600 Subject: [PATCH 7/8] 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 Date: Thu, 21 Nov 2024 17:10:46 -0600 Subject: [PATCH 8/8] 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