Issues with travel_time_matrix including reasonable transit times
Audit-uCanB opened this issue · 2 comments
Brief description of the problem:
I thought this might be because I was using the wrong GTFS file, but after reading some posts on here the other night, I found the GTFS pro resource and decided to give that a try hoping it would correct my issue, however, I'm still get a plotted map that does not make any particular sense in terms of travel time. I've used the files for the gtfs and osm.pbf available here (they are too large to attach): https://gtfs.pro/en/canada/STM/stm
I used the montreal boundaries provided in geojson WGS84 here: https://open.canada.ca/data/dataset/9797a946-9da8-41ec-8815-f6b276dec7e9
Despite using an appropriate date, I have found that when I graph the travel_time_matrix in a map of Montreal Island, the values just don't make sense -- there will be a travel time of 120 minutes for example right next to the destination point for the series: the central station. There are relatively few entries in the travel_time_matrix for less than 30minutes (maybe just 10% of the origin points), and other examples I have found online show that travel times up to 30 minutes should be pretty reasonable, with hardly anything on the island reaching above 100 minutes. This seems reasonable when considering a highly metropolitan area. Over 15% of my ttm values are over 100 minutes in contrast. Additionally, when running the travel_time_matrix, many points are dropped on account of not being within the max_trip_duration of 240 minutes which also seems like an error.
Reproducible example here
options(java.parameters = "-Xmx5G")
library(pacman)
pacman::p_load(dplyr, ggplot2, sf, terra, here, readr, stargazer, giscoR, tidyverse,r5r,osmextract,osmdata,skimr,data.table,geosphere)
montreal_divisions <-st_read(file.path(".","GTFS_Pro","Montreal-limites-administratives-agglomeration-WGS84.geojson"))
#Create one boundary
montreal_boundary <- st_union(montreal_divisions)
plot(montreal_boundary)
##Setting the size of the hexagons in hexgrid
#reference_point <- st_centroid(montreal_boundary)
desired_distance <- 250 #meter To be changed subject to resolution desires
meters_to_degrees <- 1/111000 #1 degree to 111km or 111,000 meters
hex_size <- desired_distance*meters_to_degrees
##Generation of hexgrid within selected bounding box
hex_points <- st_make_grid(
montreal_boundary,
cellsize = c(hex_size, hex_size),
square = FALSE
)
sf_use_s2(FALSE)
hex_grid <-st_transform(hex_grid, st_crs(montreal_boundary))
write.csv(hex_grid, "montreal_hexgrid_og.csv", row.names = FALSE)
hex_grid_centroids <- st_centroid(hex_grid)
hex_grid_df <- st_coordinates(hex_grid_centroids) %>%
as.data.frame() %>%
rename(lon=X, lat=Y)
#Save grid as CSV
write.csv(hex_grid_df, "montreal_hexgrid_df.csv", row.names = FALSE)
##Hard coding central station coordinates from google
point_data <-data.frame(lon = -73.5665, lat = 45.5001 )
station_central_sf <- st_as_sf(point_data,coords=c("lon", "lat"))
station_central_sf <- st_set_crs(station_central_sf, st_crs(montreal_boundary))
station_central_sf$lat <-st_coordinates(station_central_sf)[,"Y"]
station_central_sf$lon <-st_coordinates(station_central_sf)[,"X"]
station_central_sf$id <- seq_len(nrow(station_central_sf))
data_path <- file.path("C://Users//Administrator//Documents//RStudio//Spatial Analysis//Research Project//GTFS_Pro")
r5r_core2 <- setup_r5(data_path = data_path, verbose = TRUE)
points <- st_read(file.path(data_path,"montreal_hexgrid_df.csv"))
points$id <- seq_len(nrow(points))
points$id <- as.character(points$id)
#convert chr to int
points$lon <- as.numeric(points$lon)
points$lat <- as.numeric(points$lat)
# Reorder the columns with 'id' as the first column
points <- points[, c("id", setdiff(names(points), "id"))]
hex_origins <- points
station_central_df <- as.data.frame(st_coordinates(station_central_sf))
colnames(station_central_df) <- c("lon", "lat")
station_central_df$id <- seq_len(nrow(station_central_df))
station_central_df <- station_central_df[, c("id", setdiff(names(station_central_df), "id"))]
station_central_df$id <- as.character(station_central_df$id)
mode<- c("WALK","TRANSIT")
max_walk_time <- 60 #in minutes
max_trip_duration <- 240 #in minutes
departure_datetime <- as.POSIXct("05-05-2024 08:00:00",
format = "%d-%m-%Y %H:%M:%S")
#Calculating travel time matrix
ttm <- travel_time_matrix(r5r_core = r5r_core2,
origins = hex_origins,
destinations = station_central_df,
mode = mode,
departure_datetime = departure_datetime,
max_walk_time = max_walk_time,
max_trip_duration = max_trip_duration,
percentile = 50,
progress = TRUE,
time_window = 120)
##joining ttm with origin point information to graph
ttm <- ttm %>% mutate(from_id=row_number())
hex_origins <- hex_origins %>% mutate (id=row_number())
ttm_joined <- left_join(ttm,hex_origins,by=c("from_id" = "id"))
#Convert to sf for mapping
joined_ttm_sf <- st_as_sf(ttm_joined,coords=c("lon","lat"),crs=st_crs(montreal_boundary))
#Join data with hex_grid
joined_hex_grid <- left_join(hex_grid,ttm_joined,by = c("id" = "from_id"))
ggplot() +
geom_sf(data = station_central_sf, color = "red", size = 3) + # Central station
geom_sf(data = joined_hex_grid, aes(fill = travel_time_p50)) + # Travel time matrix
scale_fill_gradient(low = "green", high = "red") + # Color scale
labs(title = "Travel Time Matrix to Central Station", fill = "Travel Time (minutes)") +
theme_minimal()
Situation report here
r5r::r5r_sitrep()
> r5r::r5r_sitrep()
$r5r_package_version
[1] ‘2.0’
$r5_jar_version
[1] "7.1"
$java_version
[1] "21.0.2"
$set_memory
[1] "-Xmx5G"
$session_info
R version 4.4.1 (2024-06-14 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 22631)
Matrix products: default
locale:
[1] LC_COLLATE=English_United States.utf8
[2] LC_CTYPE=English_United States.utf8
[3] LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.utf8
time zone: Europe/Berlin
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods
[7] base
other attached packages:
[1] rJava_1.0-11 viridis_0.6.5 viridisLite_0.4.2
[4] hexbin_1.28.3 geojsonsf_2.0.3 geosphere_1.5-18
[7] data.table_1.15.4 skimr_2.1.5 osmdata_0.2.5
[10] osmextract_0.5.1 r5r_2.0 lubridate_1.9.3
[13] forcats_1.0.0 stringr_1.5.1 purrr_1.0.2
[16] tidyr_1.3.1 tibble_3.2.1 tidyverse_2.0.0
[19] giscoR_0.4.2 stargazer_5.2.3 readr_2.1.5
[22] here_1.0.1 terra_1.7-71 sf_1.0-16
[25] ggplot2_3.5.1 dplyr_1.1.4 pacman_0.5.1
loaded via a namespace (and not attached):
[1] gtable_0.3.5 xfun_0.43 lattice_0.22-6
[4] tzdb_0.4.0 vctrs_0.6.5 tools_4.4.1
[7] generics_0.1.3 proxy_0.4-27 fansi_1.0.6
[10] pkgconfig_2.0.3 KernSmooth_2.23-24 checkmate_2.3.1
[13] lifecycle_1.0.4 farver_2.1.2 compiler_4.4.1
[16] munsell_0.5.1 repr_1.1.7 codetools_0.2-20
[19] htmltools_0.5.8.1 class_7.3-22 yaml_2.3.8
[22] pillar_1.9.0 classInt_0.4-10 wk_0.9.1
[25] tidyselect_1.2.1 digest_0.6.35 stringi_1.8.4
[28] labeling_0.4.3 rprojroot_2.0.4 fastmap_1.1.1
[31] grid_4.4.1 colorspace_2.1-0 cli_3.6.2
[34] magrittr_2.0.3 base64enc_0.1-3 utf8_1.2.4
[37] e1071_1.7-14 withr_3.0.0 backports_1.4.1
[40] scales_1.3.0 sp_2.1-4 timechange_0.3.0
[43] rmarkdown_2.26 gridExtra_2.3 hms_1.1.3
[46] evaluate_0.23 knitr_1.46 s2_1.1.6
[49] rlang_1.1.3 Rcpp_1.0.12 glue_1.7.0
[52] DBI_1.2.2 rstudioapi_0.16.0 jsonlite_1.8.8
[55] R6_2.5.1 units_0.8-5
Hi @Audit-uCanB . From what I've been able to replicate your code, it seems the problem occurs because the id of hex_origins are not "compatible" with the id of the station_central_df object.
The station_central_df object has an id equal to 1. However, the location of station_central_df is not the same as the cell from hex_origins which also has id equal to 1. I believe this is messing up with the rest of the data analysis. So this is not an issue with {r5r}.
To avoid this problenm above, I would suggest using an H3 hexagonal grid. See below.
library(h3jsr)
i <- h3jsr::polygon_to_cells(geometry = montreal_boundary,
res = 9)
hexgrid <- h3jsr::cell_to_polygon(input = i)
closing this issue due to no response. I'm happy to reopen it if the problem persists