Building transit network graphs
In gtfstools: General Transit Feed Specification (GTFS) Editing and Analysing Tools

# the queries in this vignette use {dplyr}, and the graphs use {igraph},
# {tidygraph}, {sfnetworks} and {ggraph}, which gtfstools only suggests
has_dplyr <- requireNamespace("dplyr", quietly = TRUE)
has_igraph <- has_dplyr && requireNamespace("igraph", quietly = TRUE)
has_viz <- has_dplyr &&
  requireNamespace("tidygraph", quietly = TRUE) &&
  requireNamespace("sfnetworks", quietly = TRUE) &&
  requireNamespace("ggraph", quietly = TRUE) &&
  requireNamespace("ggplot2", quietly = TRUE)

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = has_dplyr
)

Sys.setenv(OMP_THREAD_LIMIT = 2)
data.table::setDTthreads(2)

Many transport analyses, such as network centrality and transit assignment models, are carried out on a graph of the transit network rather than on its timetable. build_graph() builds such a graph from the trips of a GTFS feed that depart within a period of day. The graph has two types of nodes:

stop nodes, one per stop, the physical place where a passenger stands;
service nodes, one for each route and direction at each stop it serves.

obs. this way, a stop served by routes 1 and 2, both in one direction, gets one stop node and two service nodes, all at the same coordinates.

And four types of links:

in-vehicle links between the service nodes of consecutive stops, with the mean travel time between them;
boarding links from each stop to the service nodes from which trips depart, with the number of departures and the mean waiting time (half the headway);
alighting links from each service node back to its stop;
walking transfer links between nearby stops.

The graph can be returned as tables following the General Modeling Network Specification (GMNS), or as an {igraph} or {tidygraph} graph. This vignette shows how to build a transit network graph from GTFS feeds and how to analyse and plot the network with these packages like {sfnetworks} and {ggraph}. Please see the documentation of build_graph() for the details of how nodes and links are built.

Preparing the feed

We will use the sample feed of Porto Alegre, included in the package, and build the network of the morning peak, from 7am to 10am. build_graph() pools the trips of every service in the feed. In this example, we will remove services running on weekends with filter_by_weekday():

library(gtfstools)
library(dplyr)

# read GTFS feed
path <- system.file("extdata/poa_gtfs.zip", package = "gtfstools")
gtfs <- read_gtfs(path)

# drop weekends
gtfs <- filter_by_weekday(
  gtfs,
  weekday = c("saturday", "sunday"), 
  keep = FALSE
  )

Oftentimes, a GTFS feed comes with a few inconsistencies, such as duplicated trips, unused stops, etc. This particular feed lists only the times of the first and the last stop of each trip, and the links between stops without times would have no travel time, so they couldn't be used to calculate travel times through the network. Here we will use the fix_gtfs() function to fix several of these inconsistencies. We also want to ensure that no trips exceed 25 km/h and that all trips have a minimum dwell time of 30 seconds to allow passengers to board and alight.

gtfs <- fix_gtfs(gtfs, max_speed = 25, min_dwell_time = 30)

obs. in this case, fix_gtfs() throws two warnings that refer to a few late night trips whose times go back to "00:00:00" after midnight. These trips depart outside the time of day used below, so they don't affect the network.

Building a GMNS network

By default, build_graph() returns the GMNS node, link and config tables:

g_gmns <- build_graph(
  gtfs, 
  from = "07:00:00", 
  to = "10:00:00", 
  format = "gmns"
  )

head(g_gmns$node)

head(g_gmns$link)

g_gmns$link |>
  dplyr::count(link_type)

The tables can be written to the CSV files of a GMNS network with data.table::fwrite(), as shown in the examples of build_graph(), and then read by other modelling software.

Building and igraph and Calculating the centrality of stops

With format = "igraph", the network is returned as an {igraph} graph, whose links have the travel time as an attribute. We use it to calculate the shortest travel times between every pair of stops, which include the waiting time at the stops where vehicles are boarded:

library(igraph)

g_igraph <- build_graph(
  gtfs, 
  from = "07:00:00", 
  to = "10:00:00", 
  format = "igraph"
  )

# get stop nodes
stop_nodes <- igraph::V(g_igraph)[node_type == "stop"]

# get travel times between stops
travel_times <- igraph::distances(
  g_igraph,
  v = stop_nodes,
  to = stop_nodes,
  mode = "out",
  weights = igraph::E(g_igraph)$travel_time
)

# travel time from one stop to itself is NA
# an `Inf` value means {igraph} found no path for a given origin-destination pair.
diag(travel_times) <- NA

From these travel times, we calculate the closeness centrality of each stop, which tell us how close a node is to all other nodes in the network. To do this, we calculate harmonic mean of the inverse of the travel times between each stop and the other stops. Stops that can't be reached count as 0. We don't use the centrality functions of {igraph}, such as closeness(), because they would also take the service nodes as destinations, and because they require positive weights, while alighting links have a travel time of 0. This sample feed has a single direction of each route, so the closeness from a stop depends on how far it is from the end of its routes, and the closeness to it on how far it is from their start. Hence we average the closeness from and to each stop:

closeness <- data.frame(
  stop_id = stop_nodes$stop_id,
  stop_name = stop_nodes$stop_name,
  closeness_from = rowMeans(1 / travel_times, na.rm = TRUE),
  closeness_to = colMeans(1 / travel_times, na.rm = TRUE)
) |>
  dplyr::mutate(closeness = (closeness_from + closeness_to) / 2)

closeness |>
  dplyr::slice_max(closeness, n = 10)

Plotting the network with tidygraph, sfnetworks and ggraph

With format = "tidygraph", the network is returned as a {tidygraph} graph, whose nodes have the coordinates of their stops. These coordinates can be used to plot the network with {ggraph}:

library(ggplot2)
library(ggraph)
library(tidygraph)
library(sfnetworks)

a <- build_graph(gtfs, from = "08:00:00", to = "10:00:00", format = "tidygraph")
a

ggraph(a, layout = "manual", x = x_coord, y = y_coord) +
  geom_edge_link() +
  geom_node_point()

This plot treats longitudes and latitudes as plain coordinates. To work with the network as spatial data, we convert it into an {sfnetworks} network, in which nodes are points and links are lines. Since boarding and alighting links connect nodes at the same place, we keep only the in-vehicle and transfer links, and place the in-vehicle links last, so that they are drawn on top:

net <- a |>
  tidygraph::activate("edges") |>
  dplyr::filter(link_type %in% c("in_vehicle", "transfer")) |>
  dplyr::arrange(dplyr::desc(link_type)) |>
  sfnetworks::as_sfnetwork(
    coords = c("x_coord", "y_coord"),
    crs = 4326,
    remove = FALSE,
    edges_as_lines = TRUE,
    message = FALSE
  )

net

net can be used with the functions of {sfnetworks} and {sf}, for example to filter the network by a spatial extent or to join it with other spatial data. It can also be mapped with the "sf" layout of {ggraph}:

ggraph(net, layout = "sf") +
  geom_edge_sf(aes(colour = link_type)) +
  geom_node_sf(size = 0.5) +
  scale_edge_colour_manual(
    name = "Link type",
    values = c(in_vehicle = "#0072B2", transfer = "#D55E00"),
    labels = c(in_vehicle = "In-vehicle", transfer = "Walking transfer")
  ) +
  theme_void()