diff --git a/Makefile b/Makefile index eb42337..6e05756 100644 --- a/Makefile +++ b/Makefile @@ -2,7 +2,6 @@ all: data containers cycle data: osrm-data brouter-data containers: osrm-container brouter-container -WI-cycle: WI-schools-cycle cycle: cycle_brouter walk: route_analysis.Rmd @@ -14,12 +13,6 @@ 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 diff --git a/R/data/.gitignore b/R/data/.gitignore deleted file mode 100644 index 5e7d273..0000000 --- a/R/data/.gitignore +++ /dev/null @@ -1,4 +0,0 @@ -# Ignore everything in this directory -* -# Except this file -!.gitignore diff --git a/R/functions.R b/R/functions.R index b6dc7d6..1e9d050 100644 --- a/R/functions.R +++ b/R/functions.R @@ -32,60 +32,3 @@ 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) - } - -} diff --git a/WI-schools-cycle.Rmd b/WI-schools-cycle.Rmd deleted file mode 100644 index b7971a5..0000000 --- a/WI-schools-cycle.Rmd +++ /dev/null @@ -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() -``` diff --git a/archive/route_to_school.Rmd b/archive/route_to_school.Rmd deleted file mode 100644 index c4db7e2..0000000 --- a/archive/route_to_school.Rmd +++ /dev/null @@ -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() -```