--- 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) runLTS <- FALSE logo <- image_read(path = "other/BFW_Logo_180_x_200_transparent_background.png") school_symbol <- image_read_svg(path = "other/school_FILL0_wght400_GRAD0_opsz24.svg") ``` ## School Location Data ```{r gpkg, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE} WI_schools <- st_transform(st_read(dsn = "data/Schools/Wisconsin_Public_Schools_-5986231931870160084.gpkg"), crs = 4326) WI_schools <- WI_schools %>% mutate(geom = SHAPE) ``` ## Addresses Data ```{r addresses, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE} addresses <- read_csv(file="data/addresses/Addresses_Students_EastHS_2024_GeocodeResults.csv") %>% filter(lat > 0) %>% st_as_sf(coords=c("lon","lat"), crs=4326, remove = FALSE) ``` (Remember that x = lon and y = lat.) ## Bike Level of Traffic Stress (LTS) ```{r bikelts, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE} bike_lts <- st_transform(st_read("data/bike_lts/bike_lts_DANE.geojson"), crs = 4326) # make lts attribute a factor bike_lts[["lts"]] <- as.factor(bike_lts$LTS_F) # remove segments with an LTS value of 9 bike_lts <- bike_lts %>% filter(lts != 9) # set color scale bike_lts_scale <- data.frame(code = c(1, 2, 3, 4, 9), color = c("#1a9641", "#a6d96a", "#fdae61", "#d7191c", "#d7191c")) ``` # 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" ``` ## 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 ```{r analysisPreamble, eval = TRUE, echo = FALSE, results = "show", warning = FALSE, error = TRUE, message = FALSE} radius <- 3 # miles levels <- c(1) res <- 100 threshold <- 1 ``` ## Subset Addresses Within `r radius` Miles ```{r cycleBoundary, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE} cycle_boundary_m <- radius*1609 school_focus <- data.frame(name = c("East High School"), NCES_CODE = c("550852000925")) #school_focus <- data.frame(name = c("IMAP"), NCES_CODE = c("550008203085")) cycle_boundary_poly <- st_transform(fill_holes(st_make_valid(osrmIsodistance( loc = WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE), # breaks = c(cycle_boundary_m), breaks = cycle_boundary_m*levels, res = res) ), units::set_units(threshold, km^2)), crs = 4326) addresses_near <- st_intersection(addresses, cycle_boundary_poly) ``` Notes: - _osrmIsoDistance_ is the primary function in the above chunk. - This function computes areas that are reachable within a given road distance from a point and returns the reachable regions as polygons. These areas of equal travel distance are called isodistances. - Input is a point represented as an sf object (extended data.frame-like objects with a simple feature list column) could be other classes, e.g., vector of coods, data.frame of lat tand long. etc. - Arguments to osrmIsodistances used here are breaks and res - breaks: a numeric vector of break values to define isodistance areas, in meters. - res: number of points used to compute isodistances, one side of the square grid, the total number of points will be res*res. Increase res to obtain more detailed isodistances. - _fill\_holes_ is also used with a threshold of `r threshold` km^2. - _st\_intersection_ is also used on sf objects (simple features?) ## Calculate Routes ```{r routes, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE} school_focus_location <- WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE) %>% select(LAT, LON) #calculate routes with the safety profile brouter_profile <- "safety" routes_safety <- list(NULL) for(i in addresses_near %>% arrange(number) %>% pull(number)) { query <- paste0( brouter_url, "?lonlats=", addresses_near %>% filter(number == i) %>% pull(lon), ",", addresses_near %>% filter(number == i) %>% pull(lat), "|", 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_safety[[i]] <- route_run message(paste0("done - ", i, " of ", max(addresses_near$number))) } # combine the list of routes into a data table and make sure its the right crs routes_safety <- st_transform(bind_rows(routes_safety), crs = 4326) # calculate routes with the "shortest" profile brouter_profile <- "shortest" routes_shortest <- list(NULL) for(i in addresses_near %>% arrange(number) %>% pull(number)) { query <- paste0( brouter_url, "?lonlats=", addresses_near %>% filter(number == i) %>% pull(lon), ",", addresses_near %>% filter(number == i) %>% pull(lat), "|", 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_shortest[[i]] <- route_run message(paste0("done - ", i, " of ", max(addresses_near$number))) } # combine the list of routes into a data table and make sure its the right crs routes_shortest <- st_transform(bind_rows(routes_shortest), crs = 4326) routes <- bind_rows(routes_safety, routes_shortest) %>% mutate(total.time = as.double(total.time), total.energy = as.double(total.energy), track.length = as.double(track.length)) ``` Notes: - this queries the brouter server to get routes ## plot efficiency loss ```{r plotefficiency, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE} ggplot(data = as.data.frame(routes) %>% select(student_number, name, track.length) %>% pivot_wider(names_from = name, values_from = track.length) %>% mutate(difference_percent = ((brouter_safety_0/brouter_shortest_0) - 1) * 100)) + geom_point(aes(x = brouter_shortest_0 / 1609, y = difference_percent)) + labs(x = "distance for shortest route", y = "Percent difference between safest and shortest routes") ``` ## Combine routes with Bike LTS ```{r ltscount, eval = runLTS, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE} # Count the routes that intersect or overlap with each segment of the bike_tls network. # The intersections have a buffer of 10 m bike_lts_buffer <- st_buffer(st_intersection(bike_lts, cycle_boundary_poly), 10) bike_lts_buffer["student_use"] <- unlist(lapply(st_intersects(bike_lts_buffer, routes_safety), length)) bike_lts_studentuse <- left_join(bike_lts, as.data.frame(bike_lts_buffer) %>% select(OBJECTID, student_use), by = "OBJECTID") %>% filter(student_use > 0) ``` Notes: - for each segment in bike_lts, this counts how many student’s calculated routes intersect with it (within a 10 m buffer) ```{r functions, eval = runLTS, echo = FALSE, results = "show", warning = FALSE, error = TRUE, message = FALSE} source("./R/functions.R") ``` ```{r routeslts, eval = runLTS, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE} # Start with routes_lts as a NULL list routes_lts <- list(NULL) # Pre-filter the bike_lts_buffer for relevant student use relevant_buffer <- bike_lts_buffer %>% filter(student_use > 0) routes_lts <- mclapply(addresses_near %>% arrange(number) %>% pull(number), getLTSForRoute, route_table = routes_safety, mc.cores = detectCores() / 2, mc.cleanup = TRUE, mc.preschedule = TRUE, mc.silent = FALSE) routes_lts <- bind_rows(routes_lts) ``` Notes: - for each student's route, this finds which bike_lts segment it intersects with and calculates a max and an average level of traffic stress (LTS). This takes a while, so it runs in parallel. There's probably a more efficient way to do this calculation. - see ./R/functions.R for definition of getLTSForRoute() ```{r addresslts, eval = runLTS, echo = FALSE, results = "show", warning = FALSE, error = TRUE, message = FALSE} # Join the route lts data with the addresses data addresses_near <- left_join(addresses_near, routes_lts %>% select(c("student_number", "lts_max", "lts_average", "lts_1_dist", "lts_2_dist", "lts_3_dist", "lts_4_dist")), join_by("number"=="student_number"), multiple = "any") # add supplemental analysis addresses_near <- addresses_near %>% mutate(lts_34_dist = lts_3_dist + lts_4_dist) ``` # Make Maps ## 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])) #get basemap basemap <- get_stadiamap(bbox = bbox, zoom = 15, maptype = "stamen_toner_lite") ``` Notes: - This chunk retrieves the base map from Stadia Maps (API key required) ## Generate map of addresses ```{r mapaddresses, eval = TRUE, echo = FALSE, results = "show", warning = FALSE, error = TRUE, message = FALSE} ggmap(basemap) + labs(title = paste0("Student homes at ", school_focus %>% pull(name)), x = NULL, y = NULL, color = NULL, fill = "How many students live there") + theme(axis.text=element_blank(), axis.ticks=element_blank(), plot.caption = element_text(color = "grey")) + geom_hex(data = addresses %>% extract(geometry, into = c('Lat', 'Lon'), '\\((.*),(.*)\\)', conv = T), aes(x = Lat, y = Lon), alpha = 0.7) + scale_fill_distiller(palette = "YlOrRd", direction = "reverse") + geom_sf(data = cycle_boundary_poly, inherit.aes = FALSE, aes(color = paste0(radius, " mile cycling boundary")), fill = NA, linewidth = 1) + scale_color_manual(values = "blue", name = NULL) + new_scale_color() + annotation_raster(school_symbol, # Position adjustments here using plot_box$max/min/range ymin = as.double((WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE) %>% pull(geom))[[1]])[2] - 0.001, ymax = as.double((WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE) %>% pull(geom))[[1]])[2] + 0.001, xmin = as.double((WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE) %>% pull(geom))[[1]])[1] - 0.0015, xmax = as.double((WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE) %>% pull(geom))[[1]])[1] + 0.0015) + geom_sf_label(data = WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE), inherit.aes = FALSE, mapping = aes(label = school_focus %>% pull(name)), nudge_y = 0.0015, label.size = 0.04, size = 2) ggsave(file = paste0("figures/", school_focus %>% pull(name), " Addresses_cycling.pdf"), title = paste0(school_focus %>% pull(name), " Addresses"), device = pdf, height = 8.5, width = 11, units = "in", create.dir = TRUE) ``` ## Generate map of routes ```{r maproutes, eval = TRUE, echo = FALSE, results = "show", warning = FALSE, error = TRUE, message = FALSE} # generate map ggmap(basemap) + labs(title = paste0("Cycling routes for students at ", school_focus %>% pull(name)), subtitle = paste0("only showing the safest routes for students within the ", radius, " mile cycling boundary"), x = NULL, y = NULL, color = NULL, linewidth = "Potential student cyclists") + theme(axis.text=element_blank(), axis.ticks=element_blank(), plot.caption = element_text(color = "grey")) + geom_sf(data = cycle_boundary_poly, inherit.aes = FALSE, aes(color = paste0(radius, " mile cycling boundary")), fill = NA, linewidth = 1) + scale_color_manual(values = "blue", name = NULL) + new_scale_color() + geom_sf(data = bike_lts_studentuse %>% filter(!is.na(student_use), student_use > 3), inherit.aes = FALSE, aes(linewidth = student_use), color = "mediumvioletred", fill = NA) + scale_linewidth_continuous(range = c(0, 3)) + annotation_raster(school_symbol, # Position adjustments here using plot_box$max/min/range ymin = as.double((WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE) %>% pull(geom))[[1]])[2] - 0.001, ymax = as.double((WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE) %>% pull(geom))[[1]])[2] + 0.001, xmin = as.double((WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE) %>% pull(geom))[[1]])[1] - 0.0015, xmax = as.double((WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE) %>% pull(geom))[[1]])[1] + 0.0015) + geom_sf_label(data = WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE), inherit.aes = FALSE, mapping = aes(label = school_focus %>% pull(name)), nudge_y = 0.0015, label.size = 0.04, size = 2) ggsave(file = paste0("figures/", school_focus %>% pull(name), " Routes_cycling.pdf"), title = paste0(school_focus %>% pull(name), " Cycling Routes"), device = pdf, height = 8.5, width = 11, units = "in", create.dir = TRUE) ``` ## Generate map of routes with LTS (1) ```{r maprouteslts, eval = runLTS, echo = FALSE, results = "show", warning = FALSE, error = TRUE, message = FALSE} # generate map ggmap(basemap) + labs(title = paste0("Cycling routes for students at ", school_focus %>% pull(name)), subtitle = "only showing routes within the cycling boundary", x = NULL, y = NULL, color = NULL, linewidth = "Potential student cyclists") + theme(axis.text=element_blank(), axis.ticks=element_blank(), plot.caption = element_text(color = "grey")) + geom_sf(data = cycle_boundary_poly, inherit.aes = FALSE, aes(color = paste0(radius, " mile cycling boundary")), fill = NA, linewidth = 1) + scale_color_manual(values = "blue", name = NULL) + new_scale_color() + geom_sf(data = bike_lts_studentuse %>% filter(!is.na(student_use), student_use > 0), inherit.aes = FALSE, aes(color = lts, linewidth = student_use)) + scale_color_manual(values = bike_lts_scale$color, name = "Bike Level of Traffic Stress") + scale_linewidth_continuous(range = c(0, 3)) + annotation_raster(school_symbol, # Position adjustments here using plot_box$max/min/range ymin = as.double((WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE) %>% pull(geom))[[1]])[2] - 0.001, ymax = as.double((WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE) %>% pull(geom))[[1]])[2] + 0.001, xmin = as.double((WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE) %>% pull(geom))[[1]])[1] - 0.0015, xmax = as.double((WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE) %>% pull(geom))[[1]])[1] + 0.0015) + geom_sf_label(data = WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE), inherit.aes = FALSE, mapping = aes(label = school_focus %>% pull(name)), nudge_y = 0.0015, label.size = 0.04, size = 2) ggsave(file = paste0("figures/", school_focus %>% pull(name), " Routes - Traffic Stress_cycling.pdf"), title = paste0(school_focus %>% pull(name), " Cycling Routes - Traffic Stress"), device = pdf, height = 8.5, width = 11, units = "in", create.dir = TRUE) ``` ## Generate map of routes with LTS (2) ```{r mapaddresseslts, eval = runLTS, echo = FALSE, results = "show", warning = FALSE, error = TRUE, message = FALSE} # generate map ggmap(basemap) + labs(title = paste0("Level of Traffic stress for biking for students at ", school_focus %>% pull(name)), subtitle = "only showing routes within the cycling boundary", x = NULL, y = NULL, linewidth = "Potential student cyclists") + theme(axis.text=element_blank(), axis.ticks=element_blank(), plot.caption = element_text(color = "grey")) + geom_sf(data = cycle_boundary_poly, inherit.aes = FALSE, aes(color = paste0(radius, " mile cycling boundary")), fill = NA, linewidth = 1) + scale_color_manual(values = "blue", name = NULL) + new_scale_color() + geom_sf(data = routes_lts %>% filter(route$student_use >= 3), inherit.aes = FALSE, aes(geometry = route$geometry, color = route$lts, linewidth = route$student_use)) + #scale_color_gradientn(colors = bike_lts_scale$color, name = "Length of high stress travel on route from that address", limits = c(1,4)) + scale_color_manual(values = bike_lts_scale$color, name = "Bike Level of Traffic Stress") + #scale_color_distiller(palette = "YlOrRd", direction = "reverse") + scale_linewidth_continuous(range = c(0, 3)) + annotation_raster(school_symbol, # Position adjustments here using plot_box$max/min/range ymin = as.double((WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE) %>% pull(geom))[[1]])[2] - 0.001, ymax = as.double((WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE) %>% pull(geom))[[1]])[2] + 0.001, xmin = as.double((WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE) %>% pull(geom))[[1]])[1] - 0.0015, xmax = as.double((WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE) %>% pull(geom))[[1]])[1] + 0.0015) + geom_sf_label(data = WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE), inherit.aes = FALSE, mapping = aes(label = school_focus %>% pull(name)), nudge_y = 0.0015, label.size = 0.04, size = 2) ggsave(file = paste0("figures/", school_focus %>% pull(name), " Routes - Traffic Stress_cycling_new.pdf"), title = paste0(school_focus %>% pull(name), " Student Addresses - Cycling Traffic Stress"), device = pdf, height = 8.5, width = 11, units = "in", create.dir = TRUE) ``` # Appendix ## Notes ### R Package sf - Simple features or simple feature access refers to a formal standard (ISO 19125-1:2004) that describes how objects in the real world can be represented in computers, with emphasis on the spatial geometry of these objects. It also describes how such objects can be stored in and retrieved from databases, and which geometrical operations should be defined for them. - The standard is widely implemented in spatial databases (such as PostGIS), commercial GIS (e.g., ESRI ArcGIS) and forms the vector data basis for libraries such as GDAL. A subset of simple features forms the GeoJSON standard. - R has well-supported classes for storing spatial data (sp) and interfacing to the above mentioned environments (rgdal, rgeos), but has so far lacked a complete implementation of simple features, making conversions at times convoluted, inefficient or incomplete. The package sf tries to fill this gap, and aims at succeeding sp in the long term. - See source [here](https://r-spatial.github.io/sf/articles/sf1.html) - all functions and methods in sf that operate on spatial data are prefixed by st_, which refers to spatial type; this makes them easily findable by command-line completion. - Tessellation st_make_grid() ### Sandbox #### Create a Grid Over Bikeable Area ```{r sandbox1, 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 Centroid to School ```{r sandbox2, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE} grid_pts <- st_centroid(grid) grid_coods <- st_coordinates(grid_pts) 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 ```{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)/60) total.energy.vec <- routes %>% pull(total.energy) grid <- cbind(grid, total.energy = as.numeric(total.energy.vec)) ggmap(basemap) + geom_sf(data = grid, aes(fill = total.time), inherit.aes = FALSE ) ``` #### Available Route Data ```{r sandbox4, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE} attributes(routes)$names ``` ##### 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() ``` # Archive ```{r archive1, eval = FALSE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE} # for(i in addresses_near %>% arrange(number) %>% pull(number)) { # lts_segments <- bike_lts_buffer$OBJECTID[st_intersects(bike_lts_buffer, routes %>% filter(student_number == i), sparse = FALSE)] # lts_max <- max(bike_lts_buffer %>% filter(OBJECTID %in% lts_segments) %>% pull(LTS_F), na.rm = TRUE) # lts_average <- mean(bike_lts_buffer %>% filter(OBJECTID %in% lts_segments) %>% pull(LTS_F), na.rm = TRUE) # routes_lts[[i]] <- data.frame("student_number" = c(i), "lts_max" = c(lts_max), "lts_average" = c(lts_average)) # message(paste0("done - ", i, " of ", max(addresses_near$number))) # } # routes_lts <- lapply(head(addresses_near %>% arrange(number) %>% pull(number)), # getLTSForRoute) # system.time(routes_lts <- lapply(head(addresses_near %>% arrange(number) %>% pull(number)), # getLTSForRoute)) ```