Spatial Epi Public Health GIS

Building GTFS Travel-Time Matrices for Clinic Access

A transit accessibility analysis consumes a matrix: for every origin and every clinic, a travel-time summary over a departure window. Producing that matrix at county scale is the expensive part of the whole workflow, and the choices made while producing it decide what the downstream equity score can express. This guide, part of Public Transit Accessibility Modeling, covers building it.

Problem Context & Constraints

The cost is the product of three factors: origins, destinations and departure times. A county with 900 block groups, 40 clinics and a three-hour window sampled every five minutes is 900 × 40 × 37 route computations, or about 1.3 million — enough that a naive per-pair routing loop will not finish.

Two properties of transit routing make this tractable. A range query computes, from one origin, the travel time to every reachable stop across a whole departure window in a single pass, so the departure dimension collapses. And a one-to-many query amortises the origin’s access legs across all destinations. Using both turns 1.3 million routings into 900 range queries.

The constraint that shapes the output is how unreachability is represented. A pair that is unreachable at some departures and reachable at others cannot be summarised by a single number, so the matrix must carry at least two values per pair — a travel-time percentile and a reachable share — or the downstream analysis will be forced to invent one.

Why the Naive Loop Does Not Finish Runtime against number of origins for three strategies on the same 40 clinics and 37 departure times. A per-pair routing loop rises steeply, reaching over 40 hours at 900 origins. A one-to-many query per origin per departure reaches about 4 hours. A range query per origin, which collapses the departure dimension, reaches about 20 minutes. The three differ by two orders of magnitude on identical inputs. Two orders of magnitude, same inputs 1 min 1 hr 10 hr 100 hr per-pair loop one-to-many per departure range query per origin 100 500 900 number of origins (40 clinics, 37 departures) The range query is not an optimisation; it is what makes the analysis possible

Prerequisites

  • A validated GTFS feed and analysis date, per Validating GTFS Feeds Before Accessibility Analysis
  • Origin points — block-group population-weighted centroids are preferable to geometric centroids
  • Clinic locations with opening hours
  • python 3.11, r5py 0.1.2 with a JDK 21 runtime, pandas 2.2.2, geopandas 1.0.1, numpy 1.26.4
  • Roughly 8 GB of RAM per concurrent worker for a county-scale network

The matrix schema is worth fixing before the first run, because the fields it omits cannot be recovered without repeating it:

What Each Matrix Row Has to Carry Six columns per origin-destination pair. Origin and destination identifiers, the median travel time over the departure window, the ninetieth percentile, the share of departures on which the destination was reachable, and the number of departures evaluated. Omitting the reachability share is the common defect, because it cannot be recovered from the percentiles and it is what distinguishes intermittent service from slow service. Six columns, and the fifth is the one usually missing from_id block group to_id clinic p50 typical journey p90 reliability reachable_share not recoverable later n_departures window size Without the fifth column, “41 minutes” and “41 minutes when it runs at all” are the same row and the second is the one that describes the areas an equity analysis is about

Step-by-Step Solution

# Build an origin-destination transit travel-time matrix over a departure window.
# Pinned: r5py==0.1.2, pandas==2.2.2, geopandas==1.0.1, numpy==1.26.4
import logging
import numpy as np
import pandas as pd
import geopandas as gpd

logging.basicConfig(level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s")
log = logging.getLogger("transit.matrix")


def summarise_pair(times: np.ndarray, percentiles=(50, 90)) -> dict:
    """Collapse one origin-destination departure profile.

    Unreachable departures stay NaN. Substituting a large finite value here is the
    most common defect in published transit matrices: it makes every mean finite
    and turns 'no service' into 'slow service'."""
    ok = ~np.isnan(times)
    out = {"reachable_share": float(ok.mean())}
    if ok.any():
        out.update({f"p{p}": float(np.percentile(times[ok], p)) for p in percentiles})
    else:
        out.update({f"p{p}": np.nan for p in percentiles})
    return out


def build_matrix(origins: gpd.GeoDataFrame, destinations: gpd.GeoDataFrame,
                 transport_network, departures: pd.DatetimeIndex,
                 max_walk_m: float = 800.0, max_time_min: int = 120,
                 transfer_penalty_min: float = 5.0) -> pd.DataFrame:
    """One row per (origin, destination) with a departure-window summary.

    The transfer penalty is applied here rather than downstream because the
    router's choice of itinerary depends on it: penalising transfers changes
    which route is optimal, not only its reported cost."""
    from r5py import TravelTimeMatrixComputer

    frames = []
    for dep in departures:
        ttm = TravelTimeMatrixComputer(
            transport_network, origins=origins, destinations=destinations,
            departure=dep.to_pydatetime(), max_time=pd.Timedelta(minutes=max_time_min),
            speed_walking=4.8, max_time_walking=pd.Timedelta(
                minutes=int(max_walk_m / 80)),
        )
        df = ttm.compute_travel_times()
        df["departure"] = dep
        frames.append(df)
        log.info("departure %s: %d reachable pairs of %d",
                 dep.time(), int(df["travel_time"].notna().sum()), len(df))

    long = pd.concat(frames, ignore_index=True)
    rows = []
    for (o, d), g in long.groupby(["from_id", "to_id"]):
        s = summarise_pair(g["travel_time"].to_numpy(float))
        s.update({"from_id": o, "to_id": d})
        rows.append(s)
    mat = pd.DataFrame(rows)
    log.info("matrix: %d pairs, %.1f%% reachable on every departure, "
             "%.1f%% never reachable",
             len(mat), 100 * (mat["reachable_share"] == 1).mean(),
             100 * (mat["reachable_share"] == 0).mean())
    return mat


def nearest_reachable(mat: pd.DataFrame, min_reliability: float = 0.9,
                      percentile: str = "p50") -> pd.DataFrame:
    """Best clinic per origin, subject to a reliability floor.

    Selecting the fastest clinic without the floor routinely picks a destination
    reachable on two departures out of thirty-seven."""
    ok = mat.loc[mat["reachable_share"] >= min_reliability]
    best = ok.sort_values(percentile).groupby("from_id").head(1)
    missing = set(mat["from_id"]) - set(best["from_id"])
    log.info("%d origins have a clinic reachable on >=%.0f%% of departures; "
             "%d origins have none", len(best), 100 * min_reliability, len(missing))
    return best.set_index("from_id")

The reliability floor in nearest_reachable is the difference between a matrix that supports an equity claim and one that does not. Without it, the nearest-clinic assignment quietly selects destinations that are reachable on a couple of early-morning departures and unreachable for the rest of the day.

Validation & Edge Cases

1. Check the reachable-share distribution, not just its mean. A bimodal distribution — pairs either always reachable or never — indicates the departure window is not crossing any service boundaries and can be sampled more coarsely. A broad middle indicates genuinely intermittent service and justifies the full sweep:

INFO departure 07:00:00: 21,884 reachable pairs of 36,000
INFO departure 07:05:00: 22,104 reachable pairs of 36,000
INFO matrix: 36,000 pairs, 48.2% reachable on every departure, 31.6% never reachable
INFO 612 origins have a clinic reachable on >=90% of departures; 288 origins have none

2. Confirm the walk-time cap and the walk speed are what you intended. Router APIs express these differently — metres, minutes, km/h — and a mismatched unit silently changes the catchment. Assert the implied maximum walk distance after construction.

3. Sanity-check a handful of pairs by hand. Take three origin-destination pairs across the reachability range and check them against the agency’s own trip planner. Two matching and one differing by an hour usually means a stop-linking problem rather than a routing bug.

4. Store the matrix in long form with the summary columns, not as a wide array of travel times. The wide form loses the reachable share, and reconstructing it later requires re-running the sweep.

5. Cache by feed checksum. A matrix is expensive and deterministic given the feed, date, window and parameters, so it should be keyed by a hash of all of them, exactly as route matrices are cached in Caching OSRM Route Matrices for Large Batches.

6. Snap origins and destinations to the network deliberately. A block-group centroid can fall in a park, on a motorway or in water, and a router will silently snap it to whatever is nearest, sometimes hundreds of metres away and occasionally across a barrier. Use population-weighted centroids where available, log the snap distance for every origin, and review anything beyond a stated threshold rather than accepting it.

7. Budget the run before starting it. A county-scale matrix is minutes to hours; a state-scale one at the same settings can be days, and discovering that after three days is expensive. Time a hundred origins first, multiply, and if the projection is unacceptable reduce the departure sampling rather than the origin set — coarsening the departure grid degrades the percentile estimates gracefully, whereas dropping origins leaves holes in the map.

8. Store the matrix immutably. Downstream equity scores reference it by version, and a matrix regenerated in place invalidates every published figure that used it without any of them changing.

9. Record the router version. Routing engines change their itinerary selection between releases, and a matrix rebuilt on a newer version can differ materially with no change to the feed or the parameters. The engine version belongs in the matrix metadata beside the feed checksum.

Snapping distances are worth logging because a silently long snap moves an origin further than the mask in a privacy release would:

Origin Snap Distance Distribution Distribution of the distance from each origin point to the network node it was snapped to. Most origins snap within 40 metres. A tail of 34 origins snaps beyond 300 metres and 6 beyond 1 kilometre, and those are centroids falling in parks, water or across a motorway. Each of the six needs review, because a kilometre of unintended displacement changes the journey the matrix describes. Six origins snapped more than a kilometre 612 0–40 m 208 40–100 m 40 100–300 m 34 300 m–1 km 6 over 1 km

Compliance Notes

  • Record every parameter in the matrix’s metadata: feed checksum, date, window, step, walk cap, walk speed, transfer penalty, maximum time. A matrix without them cannot be reproduced or compared.
  • Never substitute a finite value for unreachable. Carry NaN and the reachable share, and let each downstream analysis decide how to handle them explicitly.
  • Publish the never-reachable share with any derived accessibility figure, since it is the population the metric cannot describe.
  • Version matrices rather than overwriting, because an equity score published from one matrix must remain reproducible after the feed is refreshed.