Build Transit Travel-Time Matrices with Free Tools

If we want to measure how long it actually takes transit-dependent residents to reach grocery stores, we need to calculate travel times at scale. A design with 408 tracts and 6,613 store records defines about 2.7 million candidate pairs. Those input populations and the completed-run claim are under reconciliation. The method below shows how such a matrix can be calculated locally without a paid per-request routing API.

Why does this matter? Many food access questions require transit time calculations, and commercial routing APIs may impose fees or quotas on large matrices. Local routing provides another option, but it requires assembling, freezing, and validating multiple data sources. Destination coordinates remain a separate input and may come from a source with its own key, license, or fee.

This post documents how to acquire GTFS transit data, build a pedestrian network from OpenStreetMap, and route with r5py. Public-material status for the project-specific result: article only.

---

The Tool Stack

Four free, open-source components make this possible:

The Tool Stack

Component	Tool	Purpose
Transit schedules	GTFS feeds	Bus and rail routes, stops, schedules
Street network	OpenStreetMap	Walking paths between transit and destinations
Routing engine	r5py (Conveyal R5)	Calculate multimodal travel times
Destination data	Project-specific point file	Grocery store coordinates; source terms and costs vary

GTFS (General Transit Feed Specification) is the standard format for transit data. Most US transit agencies publish free GTFS feeds that include routes, stops, and schedules. ^[3]

OpenStreetMap provides free street network data worldwide. For transit routing, we need pedestrian paths: sidewalks, crosswalks, and building entrances.

r5py is a Python wrapper for Conveyal's R5 routing engine, designed for rapid accessibility analysis. ^[4] It calculates travel times accounting for walking, waiting, riding, and transferring.

---

Step 1: Acquire GTFS Data

Transit agencies publish GTFS feeds on their own schedules. California's Cal-ITP dataset can help with statewide feed discovery and analytics, but its official data page says the aggregated tables are not for trip-planner ingestion. A routing network needs a complete, internally consistent GTFS ZIP from the relevant agency or another documented routing-feed endpoint. ^[5]

Do not append /stops.csv to the Cal-ITP catalog page; that URL is metadata, not a file endpoint. For this Santa Clara County example, VTA's official host provides a mutable agency ZIP. Freeze the downloaded bytes and record their hash, retrieval time, resolved URL, and response metadata. This page does not claim that today's feed reproduces the original article run:

import hashlib
import json
import requests
from datetime import datetime, timezone
from pathlib import Path

VTA_GTFS_URL = "https://gtfs.vta.org/gtfs_vta.zip"
vta_path = Path("data/gtfs/vta_gtfs.zip")
vta_path.parent.mkdir(parents=True, exist_ok=True)

response = requests.get(VTA_GTFS_URL, timeout=60)
response.raise_for_status()
gtfs_bytes = response.content
vta_path.write_bytes(gtfs_bytes)

gtfs_manifest = {
    "requested_url": VTA_GTFS_URL,
    "resolved_url": response.url,
    "retrieved_at_utc": datetime.now(timezone.utc).isoformat(),
    "sha256": hashlib.sha256(gtfs_bytes).hexdigest(),
    "etag": response.headers.get("ETag"),
    "last_modified": response.headers.get("Last-Modified"),
}
Path("data/gtfs/vta_gtfs.manifest.json").write_text(
    json.dumps(gtfs_manifest, indent=2) + "\n",
    encoding="utf-8",
)

Key files in a GTFS feed:

• stops.txt: Stop locations (latitude, longitude)
• routes.txt: Route names and types
• trips.txt: Individual scheduled trips
• stop_times.txt: Arrival/departure times at each stop
• calendar.txt: Service patterns (weekday, weekend)

Validate the frozen archive before building the routing network. This helper is deliberately fail-fast. Removing rows from a DataFrame in memory would not change the ZIP that r5py reads, and deleting stops without updating stop_times.txt would break referential integrity. Step 5 separately checks that the pre-specified routing date has at least one active service ID in this exact archive.

def validate_gtfs(gtfs_path):
    import zipfile
    import pandas as pd
    from pathlib import Path

    required = {
        "agency.txt", "stops.txt", "routes.txt",
        "trips.txt", "stop_times.txt",
    }
    with zipfile.ZipFile(gtfs_path) as zf:
        names = set(zf.namelist())
        missing = required - names
        if missing:
            raise ValueError(f"Missing required GTFS files: {sorted(missing)}")
        if not ({"calendar.txt", "calendar_dates.txt"} & names):
            raise ValueError("Feed needs calendar.txt or calendar_dates.txt")

        stops = pd.read_csv(
            zf.open("stops.txt"),
            usecols=["stop_id", "stop_lat", "stop_lon"],
            dtype={"stop_id": "string"},
        )
        trips = pd.read_csv(
            zf.open("trips.txt"), usecols=["trip_id"],
            dtype={"trip_id": "string"},
        )
        stop_times = pd.read_csv(
            zf.open("stop_times.txt"), usecols=["trip_id", "stop_id"],
            dtype={"trip_id": "string", "stop_id": "string"},
        )

    lat = pd.to_numeric(stops["stop_lat"], errors="coerce")
    lon = pd.to_numeric(stops["stop_lon"], errors="coerce")
    bad_coordinates = (
        lat.isna() | lon.isna() |
        ~lat.between(-90, 90) | ~lon.between(-180, 180) |
        ((lat == 0) & (lon == 0))
    )
    if bad_coordinates.any():
        raise ValueError(f"Invalid coordinates for {bad_coordinates.sum()} stops")

    unknown_stops = set(stop_times["stop_id"].dropna()) - set(stops["stop_id"].dropna())
    unknown_trips = set(stop_times["trip_id"].dropna()) - set(trips["trip_id"].dropna())
    if unknown_stops or unknown_trips:
        raise ValueError(
            f"Broken references: {len(unknown_stops)} stop IDs, "
            f"{len(unknown_trips)} trip IDs"
        )

    return Path(gtfs_path)  # validated original feed, unchanged

validated_gtfs = validate_gtfs(vta_path)

---

Step 2: Download OpenStreetMap Data

r5py needs a pedestrian network to route walking segments. OpenStreetMap provides this. Geofabrik's latest URL is mutable, so a reproducible run must freeze the bytes and a manifest. Verify the provider's matching checksum before promoting the temporary download:

import hashlib
import json
from datetime import datetime, timezone
from pathlib import Path
import requests

OSM_URL = "https://download.geofabrik.de/north-america/us/california-latest.osm.pbf"
OSM_MD5_URL = OSM_URL + ".md5"
osm_path = Path("data/osm/california.osm.pbf")
temporary_path = osm_path.with_name(osm_path.name + ".part")
osm_path.parent.mkdir(parents=True, exist_ok=True)

checksum_response = requests.get(OSM_MD5_URL, timeout=60)
checksum_response.raise_for_status()
expected_md5 = checksum_response.text.split()[0].lower()

md5 = hashlib.md5()
sha256 = hashlib.sha256()
with requests.get(OSM_URL, stream=True, timeout=(30, 300)) as response:
    response.raise_for_status()
    with temporary_path.open("wb") as output:
        for chunk in response.iter_content(chunk_size=1024 * 1024):
            if not chunk:
                continue
            output.write(chunk)
            md5.update(chunk)
            sha256.update(chunk)
    resolved_url = response.url
    etag = response.headers.get("ETag")
    last_modified = response.headers.get("Last-Modified")

if md5.hexdigest().lower() != expected_md5:
    temporary_path.unlink(missing_ok=True)
    raise ValueError("Geofabrik checksum mismatch; file not promoted")

temporary_path.replace(osm_path)
osm_manifest = {
    "requested_url": OSM_URL,
    "resolved_url": resolved_url,
    "retrieved_at_utc": datetime.now(timezone.utc).isoformat(),
    "provider_md5": expected_md5,
    "sha256": sha256.hexdigest(),
    "etag": etag,
    "last_modified": last_modified,
}
Path("data/osm/california.osm.manifest.json").write_text(
    json.dumps(osm_manifest, indent=2) + "\n",
    encoding="utf-8",
)

For smaller regions, extract a subset:

# Using osmium-tool to extract a bounding box
osmium extract -b -122.5,37.0,-121.5,37.8 data/osm/california.osm.pbf -o data/osm/bay_area.osm.pbf

A regional extract reduces network-build memory and runtime relative to the statewide file. Record the exact bounding box, osmium version, input hash, output hash, and measured build time in the run record rather than treating the original article's size and timing estimates as benchmarks.

---

Step 3: Set Up r5py

This tutorial pins r5py==1.1.7. The official installation guide for that release line requires a Java Development Kit version 21 or later. A pip install does not install Java, while the conda-forge package can install an appropriate OpenJDK dependency. For a pip environment, verify Java before running Python:

java -version  # must report JDK 21 or later
python -m pip install "r5py==1.1.7"

Build the transport network:

import r5py

# Build network from GTFS and OSM
transport_network = r5py.TransportNetwork(
    osm_pbf="data/osm/bay_area.osm.pbf",
    gtfs=[str(validated_gtfs)]  # Can include multiple validated feeds
)

Build the network once and reuse this transport_network object for every batch in the same Python process. This article does not claim a persistent cross-process network cache or a verified build-time benchmark.

---

Step 4: Prepare Origins and Destinations

Origins are census tract centroids (where people live). Destinations are grocery stores (where people need to go).

import pandas as pd
import geopandas as gpd

# Load a point layer of census tract centroids
tracts = gpd.read_file("data/census_tract_centroids.geojson")
if tracts.crs is None:
    raise ValueError("Origin CRS is missing")
if not tracts.geometry.geom_type.eq("Point").all():
    raise ValueError("Origin file must contain centroid points")
tracts["GEOID"] = tracts["GEOID"].astype(str)
origins = tracts[["GEOID", "geometry"]].to_crs("EPSG:4326")
origins["id"] = origins["GEOID"]

# Load grocery store locations
stores = pd.read_csv("data/grocery_stores.csv")
destinations = gpd.GeoDataFrame(
    stores,
    geometry=gpd.points_from_xy(stores.longitude, stores.latitude),
    crs="EPSG:4326"
)
destinations["id"] = destinations["store_id"].astype(str)

Critical: Both origins and destinations need an id column and geometry column in a GeoDataFrame.

---

Step 5: Calculate Travel Times

Current r5py constructs the result with r5py.TravelTimeMatrix(...); there is no separate compute_travel_times() call. First verify that the exact frozen feed schedules service on the pre-specified analysis date. The date comes from the run environment rather than from this page, because a mutable current feed cannot be assumed to cover the original November 2024 example:

import os
import zipfile
from datetime import date, datetime, time, timedelta

def active_service_ids(gtfs_path, service_date):
    yyyymmdd = service_date.strftime("%Y%m%d")
    weekday = service_date.strftime("%A").lower()
    active = set()

    with zipfile.ZipFile(gtfs_path) as zf:
        names = set(zf.namelist())
        trips = pd.read_csv(
            zf.open("trips.txt"),
            usecols=["service_id"],
            dtype={"service_id": "string"},
        )
        trip_service_ids = set(trips["service_id"].dropna())
        if "calendar.txt" in names:
            calendar = pd.read_csv(zf.open("calendar.txt"), dtype="string")
            scheduled = calendar[
                (calendar["start_date"] <= yyyymmdd)
                & (calendar["end_date"] >= yyyymmdd)
                & (calendar[weekday] == "1")
            ]
            active.update(scheduled["service_id"].dropna())

        if "calendar_dates.txt" in names:
            exceptions = pd.read_csv(
                zf.open("calendar_dates.txt"),
                dtype={"service_id": "string", "date": "string", "exception_type": "string"},
            )
            exceptions = exceptions[exceptions["date"] == yyyymmdd]
            active.update(
                exceptions.loc[exceptions["exception_type"] == "1", "service_id"].dropna()
            )
            active.difference_update(
                exceptions.loc[exceptions["exception_type"] == "2", "service_id"].dropna()
            )

    return active & trip_service_ids

analysis_date = date.fromisoformat(os.environ["ANALYSIS_DATE"])
active_services = active_service_ids(validated_gtfs, analysis_date)
if not active_services:
    raise ValueError(
        f"The frozen GTFS has no active service on {analysis_date}; "
        "choose a pre-specified date covered by this feed or obtain the correct historical feed"
    )

departure = datetime.combine(analysis_date, time(9, 0))
routing_kwargs = {
    "departure": departure,
    "departure_time_window": timedelta(hours=2),
    "transport_modes": [r5py.TransportMode.TRANSIT],
    "access_modes": [r5py.TransportMode.WALK],
    "egress_modes": [r5py.TransportMode.WALK],
    "max_time": timedelta(minutes=60),
    "percentiles": [50],
}

travel_times = r5py.TravelTimeMatrix(
    transport_network,
    origins=origins,
    destinations=destinations,
    **routing_kwargs,
)

transport_modes=[TRANSIT] specifies the main transit leg; access_modes and egress_modes specify walking to and from transit. A direct all-walking alternative is a different choice set and should be computed as a separate matrix with transport_modes=[WALK] if the research question needs it. Do not add WALK to the main-mode list without documenting that change.

A completed constructor call returns a DataFrame with from_id, to_id, and travel_time. The article's project-specific approximately 2.7-million-pair and 45-minute execution claims are under reconciliation and are not presented here as verified outputs.

---

Step 6: Find Nearest Store by Transit

With travel times calculated, find the minimum for each origin:

# Group by origin and find minimum travel time
nearest_by_transit = (
    travel_times
    .groupby("from_id")
    .agg(
        min_transit_time=("travel_time", "min"),
        stores_within_30_min=("travel_time", lambda x: (x <= 30).sum()),
        stores_within_45_min=("travel_time", lambda x: (x <= 45).sum())
    )
    .reset_index()
)
nearest_by_transit["from_id"] = nearest_by_transit["from_id"].astype(str)

# Preserve every origin, including tracts with no reachable destination.
n_origins = len(tracts)
tracts_with_transit = tracts.merge(
    nearest_by_transit,
    left_on="GEOID",
    right_on="from_id",
    how="left",
    validate="one_to_one",
    indicator=True,
)
if len(tracts_with_transit) != n_origins:
    raise AssertionError("The merge changed the origin row count")

count_columns = ["stores_within_30_min", "stores_within_45_min"]
tracts_with_transit[count_columns] = (
    tracts_with_transit[count_columns].fillna(0).astype(int)
)
# Keep min_transit_time missing for origins with no reachable store.
# A returned row with an all-missing travel time is not reachable.
tracts_with_transit["reachable_store"] = (
    tracts_with_transit["min_transit_time"].notna()
)
tracts_with_transit = tracts_with_transit.drop(columns=["from_id", "_merge"])

An inner merge would drop origins that have no returned route. The left merge preserves them. Zero is appropriate for the count of reachable stores; it is not an appropriate travel time, so min_transit_time remains missing and the reachability flag makes the distinction explicit.

---

What Didn't Go as Expected

The implementation had three problems:

Problem 1: Memory Limits

Calculating all pairs at once exceeded available RAM. Solution: batch processing.

# Process in batches of 50 tracts
batch_size = 50
results = []

for i in range(0, len(origins), batch_size):
    batch_origins = origins.iloc[i:i+batch_size]

    batch_result = r5py.TravelTimeMatrix(
        transport_network,
        origins=batch_origins,
        destinations=destinations,
        **routing_kwargs,
    )

    results.append(batch_result)

# Combine all batches
travel_times = pd.concat(results, ignore_index=True)

Problem 2: Invalid GTFS Data

Some GTFS feeds contain errors such as coordinates at (0, 0), missing references, or inconsistent schedules. Run the fail-fast validator from Step 1 before constructing TransportNetwork. If it fails, obtain a corrected feed or run a maintained GTFS validator and repair the complete relational feed. Filtering stops in memory does not clean the ZIP and can orphan rows in stop_times.txt.

validated_gtfs = validate_gtfs(vta_path)
transport_network = r5py.TransportNetwork(
    osm_pbf="data/osm/bay_area.osm.pbf",
    gtfs=[str(validated_gtfs)],
)

Problem 3: OSM Size Limits

A statewide OSM file can exceed the memory or runtime budget of a local environment. Extract the documented region needed for the analysis, then benchmark the build:

# Calculate bounding box from your data
min_lon = min(origins.geometry.x.min(), destinations.geometry.x.min()) - 0.1
max_lon = max(origins.geometry.x.max(), destinations.geometry.x.max()) + 0.1
min_lat = min(origins.geometry.y.min(), destinations.geometry.y.min()) - 0.1
max_lat = max(origins.geometry.y.max(), destinations.geometry.y.max()) + 0.1

# Extract using osmium
import subprocess
subprocess.run([
    "osmium", "extract",
    "-b", f"{min_lon},{min_lat},{max_lon},{max_lat}",
    "data/osm/california.osm.pbf",
    "-o", "data/osm/region.osm.pbf"
], check=True)

---

When to Use This Approach

Good fit:

• Research requiring many origin-destination pairs
• Budget constraints preclude commercial APIs
• Need for reproducibility (GTFS + OSM are public data)
• Batch processing is acceptable (not real-time queries)

Less suitable:

• Real-time routing for individual trips
• Need for traffic-aware car routing (GTFS is transit only)
• Very small analyses where API costs are negligible
• Regions where transit agencies don't publish GTFS feeds (some rural areas, demand-responsive services) or where OpenStreetMap pedestrian paths are incomplete (typically rural or newly developed areas)

---

Alternatives

Commercial APIs (Google Maps, Mapbox): Easier setup and real-time capabilities, but pricing and quota rules change. Price the exact origin-destination workload against the provider's current official calculator before choosing a service.

OSRM (Open Source Routing Machine): Free and fast, but car/bike/walk only. No transit support.

OpenTripPlanner: Full-featured transit router with real-time updates and trip planning APIs. Deployment requires running a Java server and configuring multiple components (graph building, API endpoints, frontend). This overhead makes sense for transit agency websites serving individual trip queries, but for batch research calculating millions of routes, r5py's single Python script is simpler.

Conveyal Analysis: Web-based interface to R5 with built-in isochrone mapping and accessibility visualization. Handles data loading and visualization without code. The limitation: you can't easily customize routing parameters, export raw travel time matrices, or integrate results into a larger analysis pipeline.

---

Limitations

Schedule-based, not real-time: GTFS represents planned service, not actual arrivals. Delays, cancellations, and disruptions aren't captured.

Average conditions: The departure time window produces median travel times. Some trips will be faster or slower than calculated.

Walking assumptions: r5py assumes constant walking speed. Elderly residents, those with mobility limitations, or those carrying groceries may walk slower.

Network completeness: OpenStreetMap varies in quality by region. Some pedestrian paths may be missing or incorrectly mapped.

---

Public Materials

Article only. The external repository named in the original version contains a title-only README, not the code and sample data previously described.

---

Notes

Evidence note. The original cost illustration was based on a candidate-pair cross-product, not a frozen API invoice or completed-run record.

Repository note. The external repository is not designated as evidence for this article because it does not contain the described package.

[3] The official GTFS overview and documentation describe the schedule format and its required relationships. ↩

[4] Fink, C., Klumpenhouwer, W., Saraiva, M., Pereira, R. H. M., & Tenkanen, H. (2022). r5py: Rapid realistic routing with R5 in Python. Zenodo. The routing engine's accessibility design is described by Conway, M. W., Byrd, A., & van Eggermond, M. (2018), Accounting for uncertainty and variation in accessibility metrics for public transport sketch planning, Journal of Transport and Land Use, 11(1), 541-558. ↩

[5] The official Cal-ITP GTFS ingest dataset page describes statewide analytics tables and explicitly says not to use them for trip-planner ingestion. Cal-ITP's GTFS data host supports feed discovery. For this example, VTA publishes a mutable agency ZIP from its official GTFS host; a reproducible study must publish the retrieved file's hash and service-date coverage. ↩

---

Tags: #FoodSecurity #TransitAnalysis #GTFS #r5py #OpenStreetMap #Methods #Tutorial #OpenSource

---