I’m very new to RStudio, so apologies in advance for a messy, inefficient program. I’m working on plotting access to a single point in Montreal– the central station– and the plots I’m coming up with do not seem to be producing sensible output, so I feel like I’m missing some critical piece that would make this all flow together nicely. The current map produced at the end of this analysis appear to be very random–in my result, hexagons close to the station are coded as needing the max transit time to reach.

Any feedback on parts of this program that are really not necessary for this problem, anything that will help my output. I tried following the tutorial here:https://citygeographics.org/r5r-workshop/accessibility-example-travel-time-matrix/ but was not able to reproduce anything approaching this. Please see end of code for the data files used in the analysis.

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

I used the following sources of GTFS data:

• https://transitfeeds.com/p/societe-de-transport-de-montreal/39 #19Dec2023 File
• https://open.canada.ca/data/dataset/9797a946-9da8-41ec-8815-f6b276dec7e9 #WGS84 version for the montreal boundaries
• https://download.geofabrik.de/north-america/canada/quebec.html #Geofabrik Quebec OSM “latest” file
• https://www.stm.info/en/about/developers # I also downloaded the bus lines & stops shape file from here, but haven’t necessarily used that for anything in particular.