From 49eba8eb71e8dd26cbbe0369b0c8c979668f985b Mon Sep 17 00:00:00 2001 From: Ben Varick Date: Thu, 31 Oct 2024 16:39:12 -0500 Subject: [PATCH] removed unnecessary st_transform steps --- R/route_analysis.Rmd | 36 +- R/route_analysis.html | 1876 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 1895 insertions(+), 17 deletions(-) create mode 100644 R/route_analysis.html diff --git a/R/route_analysis.Rmd b/R/route_analysis.Rmd index 0a1729a..34c06c4 100644 --- a/R/route_analysis.Rmd +++ b/R/route_analysis.Rmd @@ -7,6 +7,8 @@ output: toc_float: collapsed: false smooth_scroll: true +editor_options: + chunk_output_type: console --- ```{r libs, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE} @@ -51,7 +53,7 @@ register_stadiamaps(key = substr(read_file(file = "api_keys/stadia_api_key"), 1, ## subset addresses within 1.5 miles walk_boundary_poly <- fill_holes(st_make_valid(osrmIsodistance( - loc = st_transform(WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE), crs = 4326), + loc = WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE), breaks = c(walk_boundary_m), res = 80) ), units::set_units(1, km^2)) @@ -94,7 +96,7 @@ 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") -bbox <- st_bbox(st_transform(st_buffer(walk_boundary_poly, dist = 500), crs = 4326)) +bbox <- st_bbox(st_buffer(walk_boundary_poly, dist = 500)) bbox <- c(left = as.double(bbox[1]), bottom = as.double(bbox[2]), right = as.double(bbox[3]), @@ -130,11 +132,11 @@ ggmap(basemap) + scale_linewidth_continuous(range = c(0, 3)) + annotation_raster(school_symbol, # Position adjustments here using plot_box$max/min/range - ymin = as.double((st_transform(WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE), crs = 4326) %>% pull(geom))[[1]])[2] - 0.001, - ymax = as.double((st_transform(WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE), crs = 4326) %>% pull(geom))[[1]])[2] + 0.001, - xmin = as.double((st_transform(WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE), crs = 4326) %>% pull(geom))[[1]])[1] - 0.0015, - xmax = as.double((st_transform(WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE), crs = 4326) %>% pull(geom))[[1]])[1] + 0.0015) + - geom_sf_label(data = st_transform(WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE), crs = 4326), + 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, @@ -178,11 +180,11 @@ ggmap(basemap) + scale_linewidth_continuous(range = c(0, 3)) + annotation_raster(school_symbol, # Position adjustments here using plot_box$max/min/range - ymin = as.double((st_transform(WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE), crs = 4326) %>% pull(geom))[[1]])[2] - 0.001, - ymax = as.double((st_transform(WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE), crs = 4326) %>% pull(geom))[[1]])[2] + 0.001, - xmin = as.double((st_transform(WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE), crs = 4326) %>% pull(geom))[[1]])[1] - 0.0015, - xmax = as.double((st_transform(WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE), crs = 4326) %>% pull(geom))[[1]])[1] + 0.0015) + - geom_sf_label(data = st_transform(WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE), crs = 4326), + 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, @@ -222,11 +224,11 @@ ggmap(basemap) + new_scale_color() + annotation_raster(school_symbol, # Position adjustments here using plot_box$max/min/range - ymin = as.double((st_transform(WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE), crs = 4326) %>% pull(geom))[[1]])[2] - 0.001, - ymax = as.double((st_transform(WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE), crs = 4326) %>% pull(geom))[[1]])[2] + 0.001, - xmin = as.double((st_transform(WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE), crs = 4326) %>% pull(geom))[[1]])[1] - 0.0015, - xmax = as.double((st_transform(WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE), crs = 4326) %>% pull(geom))[[1]])[1] + 0.0015) + - geom_sf_label(data = st_transform(WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE), crs = 4326), + 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, diff --git a/R/route_analysis.html b/R/route_analysis.html new file mode 100644 index 0000000..f91f004 --- /dev/null +++ b/R/route_analysis.html @@ -0,0 +1,1876 @@ + + + + + + + + + + + + + +Route Analysis + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ + + +
+
+
+
+
+ +
+ + + + + + + +
date()
+
## [1] "Thu Oct 31 16:37:36 2024"
+
rm(list=ls())
+library(tidyverse)
+library(ggmap)
+library(sf)
+library(osrm)
+library(smoothr)
+library(magick)
+library(ggnewscale)
+library(rsvg)
+fig.height <- 6
+set.seed(1)
+
+

Main R script

+
## school focus
+school_focus <- data.frame(name = c("East High School"), NCES_CODE = c("550852000925"))
+
+## walk boundary
+walk_boundary_m <- 1.5 * 1609
+
+## load school locations
+
+WI_schools <- st_transform(st_read(dsn = "data/Schools/Wisconsin_Public_Schools_-5986231931870160084.gpkg"), crs = 4326)
+
## Reading layer `WI_Public_Schools' from data source 
+##   `/home/ben/Bike Fed/Documents/data_analysis/route_analysis/data/Schools/Wisconsin_Public_Schools_-5986231931870160084.gpkg' 
+##   using driver `GPKG'
+## replacing null geometries with empty geometries
+## Simple feature collection with 2274 features and 32 fields (with 120 geometries empty)
+## Geometry type: POINT
+## Dimension:     XY
+## Bounding box:  xmin: 298957.7 ymin: 225953.4 xmax: 760547.6 ymax: 704946.8
+## Projected CRS: NAD83(HARN) / Wisconsin Transverse Mercator
+
WI_schools <- WI_schools %>% mutate(geom = SHAPE)
+
+## load addresses
+addresses <- read_csv(file="data/addresses/Addresses_Students_EastHS_2024_GeocodeResults.csv") %>%
+  filter(lat > 0) %>%
+  st_as_sf(coords=c("lon","lat"), crs=4326) # remember x=lon and y=lat
+
+## set osrm options
+options(osrm.server = "http://127.0.0.1:5000/")
+options(osrm.profile = "walk")
+
+register_stadiamaps(key = substr(read_file(file = "api_keys/stadia_api_key"), 1, 36))
+
+## subset addresses within 1.5 miles
+walk_boundary_poly <- fill_holes(st_make_valid(osrmIsodistance(
+  loc = WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE),
+  breaks = c(walk_boundary_m),
+  res = 80)
+), units::set_units(1, km^2))
+
+addresses_near <- st_intersection(addresses, walk_boundary_poly)
+
+## load bike tls
+bike_lts <- st_read("data/bike_lts/bike_lts_DANE.geojson")
+
## Reading layer `Bike_LTS' from data source 
+##   `/home/ben/Bike Fed/Documents/data_analysis/route_analysis/data/bike_lts/bike_lts_DANE.geojson' 
+##   using driver `GeoJSON'
+## Simple feature collection with 43355 features and 5 fields
+## Geometry type: MULTILINESTRING
+## Dimension:     XY
+## Bounding box:  xmin: -89.84212 ymin: 42.8404 xmax: -89.0082 ymax: 43.2942
+## Geodetic CRS:  WGS 84
+
bike_lts[["lts"]] <- as.factor(bike_lts$LTS_F)
+
+bike_lts_scale <- data.frame(code = c(1, 2, 3, 4, 9),
+                             color = c("#1a9641",
+                                       "#a6d96a",
+                                       "#fdae61",
+                                       "#d7191c",
+                                       "#d7191c"))
+
+## calculate routes
+routes <- list(NULL)
+for(i in addresses_near$number) {
+ routes[[i]] <- osrmRoute(
+      src = addresses_near %>% filter(number == i),
+      dst = WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE))
+  message(paste0("done - ", i, "of", max(addresses_near$number)))
+}
+routes <- bind_rows(routes)
+
+## combine routes
+# Count the routes that intersect or overlap with each segment of the bike_tls network.
+# The intersections have a buffer of 20m
+bike_lts_buffer <- st_buffer(st_intersection(bike_lts, walk_boundary_poly), 20)
+
+bike_lts_buffer["student_use"] <- unlist(lapply(st_intersects(bike_lts_buffer, routes), length))
+
+bike_lts <- st_join(bike_lts, bike_lts_buffer %>% select(OBJECTID, student_use))
+
+## make maps
+# load logo
+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")
+
+
+bbox <- st_bbox(st_buffer(walk_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")
+
+# generate map
+ggmap(basemap) +
+  labs(title = paste0("Walking routes for students at ",
+    school_focus %>% pull(name)),
+    subtitle = "only showing routes within the 1.5 walk boundary",
+    x = NULL,
+    y = NULL,
+    color = NULL,
+    linewidth = "Potential student walkers") +
+  theme(axis.text=element_blank(),
+        axis.ticks=element_blank(),
+        plot.caption = element_text(color = "grey")) +
+  geom_sf(data = walk_boundary_poly,
+          inherit.aes = FALSE,
+          aes(color = paste0(1.5, " mile walking boundary")),
+          fill = NA,
+          linewidth = 1) +
+  scale_color_manual(values = "blue", name = NULL) +
+  new_scale_color() +
+  geom_sf(data = bike_lts %>% 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.pdf"),
+       title = paste0(school_focus %>% pull(name), " Walking Routes"),
+       device = pdf,
+       height = 8.5,
+       width = 11,
+       units = "in",
+       create.dir = TRUE)
+
+# generate map
+ggmap(basemap) +
+  labs(title = paste0("Walking routes for students at ",
+                      school_focus %>% pull(name)),
+       subtitle = "only showing routes within the walk boundary",
+       x = NULL,
+       y = NULL,
+       color = NULL,
+       linewidth = "Potential student walkers") +
+  theme(axis.text=element_blank(),
+        axis.ticks=element_blank(),
+        plot.caption = element_text(color = "grey")) +
+  geom_sf(data = walk_boundary_poly,
+          inherit.aes = FALSE,
+          aes(color = paste0(1.5, " mile walking boundary")),
+          fill = NA,
+          linewidth = 1) +
+  scale_color_manual(values = "blue", name = NULL) +
+  new_scale_color() +
+  geom_sf(data = bike_lts %>% 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.pdf"),
+       title = paste0(school_focus %>% pull(name), " Walking Routes - Traffic Stress"),
+       device = pdf,
+       height = 8.5,
+       width = 11,
+       units = "in",
+       create.dir = TRUE)
+
+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 = walk_boundary_poly,
+          inherit.aes = FALSE,
+          aes(color = paste0(1.5, " mile walking 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.pdf"),
+       title = paste0(school_focus %>% pull(name), " Addresses"),
+       device = pdf,
+       height = 8.5,
+       width = 11,
+       units = "in",
+       create.dir = TRUE)
+
+
+

Appendix

+
date()
+
## [1] "Thu Oct 31 16:38:18 2024"
+
sessionInfo()
+
## R version 4.1.2 (2021-11-01)
+## Platform: x86_64-pc-linux-gnu (64-bit)
+## Running under: Pop!_OS 22.04 LTS
+## 
+## Matrix products: default
+## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
+## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
+## 
+## locale:
+##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
+##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
+##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
+##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
+##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
+## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
+## 
+## attached base packages:
+## [1] stats     graphics  grDevices utils     datasets  methods   base     
+## 
+## other attached packages:
+##  [1] rsvg_2.6.1       ggnewscale_0.5.0 magick_2.8.5     smoothr_1.0.1   
+##  [5] osrm_4.2.0       sf_1.0-18        ggmap_4.0.0      lubridate_1.9.3 
+##  [9] forcats_1.0.0    stringr_1.5.1    dplyr_1.1.4      purrr_1.0.2     
+## [13] readr_2.1.5      tidyr_1.3.1      tibble_3.2.1     ggplot2_3.5.1   
+## [17] tidyverse_2.0.0 
+## 
+## loaded via a namespace (and not attached):
+##  [1] httr_1.4.7            sass_0.4.9            bit64_4.5.2          
+##  [4] vroom_1.6.5           jsonlite_1.8.9        bslib_0.8.0          
+##  [7] highr_0.11            yaml_2.3.10           lattice_0.20-45      
+## [10] pillar_1.9.0          glue_1.8.0            digest_0.6.37        
+## [13] RColorBrewer_1.1-3    colorspace_2.1-1      htmltools_0.5.8.1    
+## [16] plyr_1.8.9            pkgconfig_2.0.3       s2_1.1.7             
+## [19] scales_1.3.0          terra_1.7-83          jpeg_0.1-10          
+## [22] tzdb_0.4.0            timechange_0.3.0      proxy_0.4-27         
+## [25] generics_0.1.3        farver_2.1.2          googlePolylines_0.8.5
+## [28] cachem_1.1.0          withr_3.0.2           hexbin_1.28.4        
+## [31] cli_3.6.3             magrittr_2.0.3        crayon_1.5.3         
+## [34] evaluate_1.0.1        fansi_1.0.6           class_7.3-20         
+## [37] tools_4.1.2           hms_1.1.3             RgoogleMaps_1.5.1    
+## [40] mapiso_0.3.0          lifecycle_1.0.4       munsell_0.5.1        
+## [43] isoband_0.2.7         compiler_4.1.2        jquerylib_0.1.4      
+## [46] e1071_1.7-16          rlang_1.1.4           classInt_0.4-10      
+## [49] units_0.8-5           grid_4.1.2            rstudioapi_0.17.1    
+## [52] bitops_1.0-9          labeling_0.4.3        rmarkdown_2.28       
+## [55] wk_0.9.4              gtable_0.3.6          codetools_0.2-18     
+## [58] DBI_1.2.3             curl_5.2.3            RcppSimdJson_0.1.12  
+## [61] R6_2.5.1              knitr_1.48            fastmap_1.2.0        
+## [64] bit_4.5.0             utf8_1.2.4            KernSmooth_2.23-20   
+## [67] stringi_1.8.4         parallel_4.1.2        Rcpp_1.0.13          
+## [70] vctrs_0.6.5           png_0.1-8             tidyselect_1.2.1     
+## [73] xfun_0.48
+
+ + + +
+
+ +
+ + + + + + + + + + + + + + + +