Snapping GTFS Stops to Shapes in Python

Project both stops and shapes to a metric CRS, build one LineString per shape_id, then use LineString.project() to get each stop’s distance along its route shape and LineString.interpolate() to get the exact point on the line. The perpendicular distance between the original stop and that snapped point is a recomputed, trustworthy shape_dist_traveled plus an off-shape offset you can threshold for data quality. Never run this on raw WGS84 coordinates — the numbers come back in degrees.

Root Cause Analysis

Snapping is the operation of finding, for each stop, the nearest position along the route line it belongs to. Two things make it non-trivial in GTFS: the geometry is not pre-built, and the coordinates are angular.

The geometry gap is that a stop and a shape share no key. stop_times.txt links a trip_id to a stop_id, and trips.txt links that trip_id to a shape_id — so a stop’s shape is only reachable by joining through the trip. Two stops on the same physical corner can belong to different shapes depending on which trips serve them, so snapping is per-(trip, stop), resolved to per-(shape, stop) once you have the mapping. The stops.txt and stop_times.txt relationship guide covers those keys in depth.

The coordinate gap is more insidious. shapely’s project and interpolate operate in whatever units the geometry carries. On EPSG:4326 those units are degrees, so project returns a “distance” of, say, 0.037 — a degree measure that is not metres, not comparable across latitudes, and useless as shape_dist_traveled. The fix is to project both layers to a shared metric coordinate reference system — a local UTM zone or national grid chosen per the projected-CRS selection guide — before any linear referencing call.

Production-Ready Python Implementation

The script loads the four relevant files, projects to an automatically derived UTM zone, builds a LineString per shape_id, and snaps every distinct (shape_id, stop_id) pair. It emits distance-along-shape and off-shape offset, then flags stops beyond a configurable threshold.

import math
import logging
import numpy as np
import pandas as pd
import geopandas as gpd
from shapely.geometry import LineString, Point

logging.basicConfig(level=logging.INFO, format="%(levelname)s: %(message)s")


def utm_epsg(lon: float, lat: float) -> int:
    """EPSG code for the UTM zone containing (lon, lat)."""
    zone = math.floor((lon + 180) / 6) + 1
    return (32600 if lat >= 0 else 32700) + zone


def snap_stops_to_shapes(feed_dir: str, max_offset_m: float = 40.0) -> pd.DataFrame:
    """Snap each stop onto its trip's shape and return distance-along + offset."""
    # --- Load with explicit dtypes ---
    stops = pd.read_csv(
        f"{feed_dir}/stops.txt",
        dtype={"stop_id": str, "stop_lat": float, "stop_lon": float},
        usecols=["stop_id", "stop_lat", "stop_lon"],
    )
    trips = pd.read_csv(
        f"{feed_dir}/trips.txt",
        dtype={"trip_id": str, "route_id": str, "shape_id": str},
        usecols=["trip_id", "route_id", "shape_id"],
    )
    stop_times = pd.read_csv(
        f"{feed_dir}/stop_times.txt",
        dtype={"trip_id": str, "stop_id": str, "stop_sequence": "Int64"},
        usecols=["trip_id", "stop_id", "stop_sequence"],
    )
    shapes = pd.read_csv(
        f"{feed_dir}/shapes.txt",
        dtype={
            "shape_id": str,
            "shape_pt_lat": float,
            "shape_pt_lon": float,
            "shape_pt_sequence": int,
        },
        usecols=["shape_id", "shape_pt_lat", "shape_pt_lon", "shape_pt_sequence"],
    )

    # --- Derive the metric CRS from the feed centroid ---
    target_epsg = utm_epsg(float(stops["stop_lon"].mean()), float(stops["stop_lat"].mean()))
    logging.info("Projecting stops and shapes to EPSG:%d", target_epsg)

    # --- Distinct (shape_id, stop_id) pairs reachable via trips ---
    st_shape = stop_times.merge(trips[["trip_id", "shape_id"]], on="trip_id", how="inner")
    pairs = (
        st_shape.dropna(subset=["shape_id"])[["shape_id", "stop_id"]]
        .drop_duplicates()
        .reset_index(drop=True)
    )
    logging.info("Snapping %d distinct (shape_id, stop_id) pairs", len(pairs))

    # --- Project stops once; index by stop_id ---
    stops_gdf = gpd.GeoDataFrame(
        stops,
        geometry=gpd.points_from_xy(stops["stop_lon"], stops["stop_lat"]),
        crs="EPSG:4326",
    ).to_crs(target_epsg)
    stop_geom = dict(zip(stops_gdf["stop_id"], stops_gdf.geometry))

    # --- Build one projected LineString per shape_id ---
    shapes = shapes.sort_values(["shape_id", "shape_pt_sequence"], kind="stable")
    shape_lines_wgs = (
        shapes.groupby("shape_id", sort=False)
        .apply(lambda g: LineString(zip(g["shape_pt_lon"], g["shape_pt_lat"])))
        .rename("geometry")
    )
    shape_lines = (
        gpd.GeoDataFrame(shape_lines_wgs, geometry="geometry", crs="EPSG:4326")
        .to_crs(target_epsg)
    )
    line_geom = dict(zip(shape_lines.index, shape_lines.geometry))

    # --- Linear referencing per pair ---
    records = []
    for shape_id, stop_id in pairs.itertuples(index=False):
        line = line_geom.get(shape_id)
        point = stop_geom.get(stop_id)
        if line is None or point is None or len(line.coords) < 2:
            continue
        dist_along = line.project(point)          # metres from shape start
        snapped: Point = line.interpolate(dist_along)
        offset = point.distance(snapped)          # perpendicular metres
        records.append(
            {
                "shape_id": shape_id,
                "stop_id": stop_id,
                "shape_dist_traveled_m": round(dist_along, 3),
                "off_shape_offset_m": round(offset, 3),
                "snapped_lon_m": round(snapped.x, 3),
                "snapped_lat_m": round(snapped.y, 3),
            }
        )

    result = pd.DataFrame.from_records(records)
    flagged = result[result["off_shape_offset_m"] > max_offset_m]
    logging.info(
        "Snapped %d pairs | %d exceed %.0f m off-shape",
        len(result), len(flagged), max_offset_m,
    )
    return result


# --- Entry point ---
# snapped = snap_stops_to_shapes("/data/gtfs/mbta", max_offset_m=40.0)
# snapped.sort_values("off_shape_offset_m", ascending=False).head(20)

Step-by-Step Walkthrough

Deriving the CRS from the data. utm_epsg turns the feed’s mean longitude and latitude into a UTM zone code (32600 + zone north, 32700 + zone south). Hard-coding a zone breaks the moment you run a feed from another region; deriving it keeps the script portable while still yielding metre units.

Resolving stops to shapes through trips. A stop’s geometry is only reachable through its trips, so the script merges stop_times.txt with trips.txt on trip_id, then reduces to the distinct (shape_id, stop_id) pairs. This is what makes snapping correct for feeds where the same physical stop serves several routes with different shapes.

Projecting before referencing. Both the stop points and the shape lines are converted with to_crs(target_epsg) before any project/interpolate call. This is the single most important line: skip it and every distance below is in degrees. The dictionaries stop_geom and line_geom cache the projected geometries so each is built once, not once per pair.

project then interpolate. line.project(point) returns the along-line distance in metres to the foot of the perpendicular from the stop. line.interpolate(dist_along) returns the coordinate at that distance — the snapped point. point.distance(snapped) is the perpendicular offset. The first value replaces the untrustworthy shape_dist_traveled; the second is a pure data-quality metric.

Thresholding the offset. Anything beyond max_offset_m is surfaced separately. For a dense urban feed, 30–50 m catches genuinely misplaced stops without drowning you in false positives; widen it for rural or rail feeds where shapes are coarser.

Verification and Output

Assert the invariants before trusting the output downstream:

snapped = snap_stops_to_shapes("/data/gtfs/mbta", max_offset_m=40.0)

# Distance-along-shape must be non-negative and finite
assert (snapped["shape_dist_traveled_m"] >= 0).all(), "Negative distance-along-shape"
assert np.isfinite(snapped["off_shape_offset_m"]).all(), "Non-finite offset present"

# Offsets in metres are single/double digits for well-formed urban feeds
median_offset = snapped["off_shape_offset_m"].median()
assert median_offset < 25, (
    f"Median off-shape offset {median_offset:.1f} m is high — "
    "check that projection ran (degrees leaking through) or shape matching is wrong"
)

worst = snapped.nlargest(5, "off_shape_offset_m")
print(worst[["shape_id", "stop_id", "shape_dist_traveled_m", "off_shape_offset_m"]]
      .to_string(index=False))

A healthy output shows shape_dist_traveled_m increasing monotonically along a trip’s stops and off_shape_offset_m in the low tens of metres. A median offset in the thousands is the classic symptom of forgetting the projection step — the numbers are degrees multiplied by the projected scale.

Gotchas and Edge Cases

• Self-intersecting shapes cause ambiguous projection. When a shape crosses itself (a loop route, or a shape that runs out and back on the same street), project returns the distance to the globally nearest point, which may be on the wrong pass. For loop routes, split the shape at the terminus or snap within a sliding distance window rather than against the whole line.

• Stops shared across shapes get one row per shape. Because snapping is per-(shape_id, stop_id), a stop served by five shapes produces five rows with five different shape_dist_traveled_m values — which is correct, since the same stop sits at a different distance along each route. Do not deduplicate on stop_id alone.

• Degenerate single-point shapes. A shape with one vertex (an export bug) makes a zero-length line on which project always returns 0. The len(line.coords) < 2 guard skips these; audit them separately, because they indicate a broken shapes.txt export.

• Times above 24:00:00 are irrelevant here but co-occur with bad feeds. Snapping ignores arrival_time, but feeds sloppy enough to ship overnight times as invalid are often the same feeds with unsorted shape points — sort shape_pt_sequence as an integer, never lexically, or the assembled line zig-zags and every offset inflates.

Frequently Asked Questions

Related

• Extracting Route Geometry from GTFS Shapes — build the LineStrings this snapping step consumes
• GeoPandas Spatial Joins for Transit Stops — take the snapped, projected stops into polygon joins
• Up: Spatial Analysis and Route Geometry — where linear referencing fits in the full geometry pipeline
• Section: GTFS Feed Architecture & Fundamentals — the schema foundation under every spatial operation · Home