在RStudio中用r5r绘制公交可达性：出行时间与地图问题排查求助

核心问题 #

我在开展蒙特利尔到中央车站的公交可达性可视化工作时，生成旅行时间矩阵（TTM）环节出现多处异常：

• 初始hex_grid_points有1123个，但生成TTM后仅剩余934个，原因不明。
• 将结果叠加到蒙特利尔地图时，出行时间与到中央车站的距离完全不匹配，且中央车站未按预期显示为红点。

排查尝试 #

我尝试调整参数排查问题：

• 把max_trip_duration增大到300（此数值不符合实际）后，TTM中的点数变为1050个，但仍未达到预期的1123个。
• 调整参数后出现新问题：相邻1km大小的六边形出行时间差异极大——一个需要100分钟，相邻的仅需10分钟，而步行在两者间仅需10-15分钟即可通行。
• 我参照相关教程操作，生成的地图出行时间差异显著，尽管研究日期略有不同，但结果不应差距如此之大。

代码片段 #

library(pacman)
pacman::p_load(dplyr, ggplot2, sf, terra, here, readr, stargazer, giscoR, tidyverse,r5r,osmextract,osmdata,skimr,data.table,geosphere, geojsonsf, hexbin, viridis)

options(java.parameters = "-Xmx5G")

#Creating Single Polygon to represent bounds of sample area
montreal_divisions <-st_read(file.path(".","Montreal Data_transitfeeds","Montreal-limites-administratives-agglomeration-WGS84.geojson"))
montreal_boundary <- st_union(montreal_divisions)
montrealbbox <- st_bbox(montreal_boundary)

#Creation of Hexgrid
#CRS = 4326 , converting degrees to meters in order to set size 
desired_distance <- 1000 #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

hex_points <- st_make_grid(
  montreal_boundary, 
  cellsize = c(hex_size, hex_size),
  square = FALSE
  )

sf_use_s2(FALSE)
hex_grid <- st_intersection(st_as_sf(hex_points), montreal_boundary)
#Add ID column
hex_grid$id <-seq_len(nrow(hex_grid))
#Set ID column as first column
hex_grid <- hex_grid[, c("id", setdiff(names(hex_grid), "id"))]

write.csv(hex_grid, "montreal_hexgrid_og.csv", row.names = FALSE)

#Convert hex_grid to dataframe of centroids
hex_grid_centroids <- st_centroid(hex_grid)
hex_grid_df <- st_coordinates(hex_grid_centroids) %>%
  as.data.frame() %>%
  rename(lon=X, lat=Y)
write.csv(hex_grid_df, "montreal_hexgrid_df.csv", row.names = FALSE)

##Hard coding coordinates of destination = Central Station Montreal (source = 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,"+proj=longlat +datum=WGS84")

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))

#Build transport network
data_path <- file.path("C:\\")  #Containing file with GTFS.zip and osm.pbf files
r5r_core <- setup_r5(data_path = data_path)


#Travel Time Matrix

points <- st_read(file.path(data_path,"montreal_hexgrid_df.csv")) 
points$id <- seq_len(nrow(points))
#convert chr to int
points$lon <- as.numeric(points$lon)
points$lat <- as.numeric(points$lat)


mode<- c("WALK", "TRANSIT") 
max_walk_time <- 30 #in minutes
max_trip_duration <- 120 #in minutes
departure_datetime <- as.POSIXct("28-11-2023 08:00:00", 
                                 format = "%d-%m-%Y %H:%M:%S")

#Calculating travel time matrix
ttm <- travel_time_matrix(r5r_core = r5r_core,
                          origins = points,
                          destinations = station_central_sf,
                          mode = mode,
                          departure_datetime = departure_datetime,
                          max_walk_time = max_walk_time,
                          max_trip_duration = max_trip_duration,
                          progress = TRUE,
                          time_window = 60)

write.csv(ttm,file.path(data_path,"CentralStationDestTTM.csv"),row.names=FALSE)


#Producing a couple of graphs for troubleshooting
ggplot(ttm, aes(x = from_id, y = to_id, fill = travel_time_p50)) +
  geom_tile() +
  scale_fill_gradient(low = "lightblue", high = "darkblue") +
  labs(x = "Origin", y = "Destination", title = "Travel Time Matrix") +
  theme_minimal()

##Convert ttm to a data frame
ttm_df <- as.data.frame(ttm)

##Plot travel times as a bar graph
ggplot(ttm_df, aes(x = from_id, y = travel_time_p50)) +
  geom_bar(stat = "identity", fill = "skyblue") +
  labs(x = "Origin", y = "Travel Time (minutes)", title = "Travel Time to Central Station") +
  theme_minimal()

#Finally attempting to Graph the matrix values

points <- points[, c("id", setdiff(names(points), "id"))] ##Reorder the columns with 'id' as the first column
ttm <- ttm %>% mutate(from_id=row_number())
points <- points %>% mutate (id=row_number())
ttm_joined <- left_join(ttm,points,by=c("from_id" = "id"))

joined_ttm_sf <- st_as_sf(ttm_joined,coords=c("lon","lat"),crs=4326) #Convert to sf for mapping

#Join data with hex_grid 
joined_hex_grid <- left_join(hex_grid,ttm_joined,by = c("id" = "from_id"))
                    

breaks <- seq(0, max(joined_hex_grid$travel_time_p50, na.rm = TRUE), by = 10) #Define breaks for the color scale

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",breaks=breaks) +     # Color scale
  labs(title = "Travel Time Matrix to Central Station", fill = "Travel Time (minutes)") +
  theme_minimal()

使用的数据源 #

• STM公交GTFS数据（2023年12月19日）
• 蒙特利尔行政边界WGS84格式数据
• Geofabrik魁北克OSM数据
• STM公交线和站点Shapefile（暂未使用）

内容的提问来源于stack exchange，提问作者Audit_power28