Spatial Analysis and Route Geometry

Every map view, on-street distance, and walk-catchment metric a transit agency publishes depends on turning shapes.txt point sequences into real geometry and relating stops to those lines. The GTFS Feed Architecture & Fundamentals overview covers the schema constraints that make a feed loadable; this guide covers the spatial layer built on top of it: assembling a LineString per route from ordered shape points, snapping stops onto those lines with linear referencing, and joining stops to polygon layers with geopandas. The workflow is production-grade for analysts producing service-area reports and for engineers wiring GTFS into a spatial warehouse.

The Spatial Problem in GTFS Geometry

GTFS stores geometry as flat rows, not as objects. A route’s path lives in shapes.txt as thousands of (shape_pt_lat, shape_pt_lon) points tagged with a shape_id and an ordering key, shape_pt_sequence. Nothing in the file format guarantees those rows arrive sorted, that a route maps to exactly one shape, or that the optional shape_dist_traveled column is populated in a unit you can trust. Turning that into usable geometry is an assembly problem, and every step has a failure mode that survives schema validation.

Compounding this, the coordinates are angular. GTFS mandates WGS84 for all latitude and longitude values, which is correct for data exchange but wrong for measurement: LineString.length on degree coordinates returns a number that is neither metres nor miles and changes with latitude. Any pipeline that reports route lengths, snaps stops by distance, or buffers walk catchments must first project out of WGS84 into a metric system, do the analysis, and — if the output is a GTFS feed again — revert. The three topics in this section each isolate one stage of that pipeline.

Prerequisites

Before running the code across this section, confirm the following:

• Python 3.9+ with pandas, geopandas, shapely (2.x), pyproj, and numpy installed:

pip install pandas geopandas shapely pyproj numpy

• A GTFS static feed containing shapes.txt, trips.txt, stops.txt, and stop_times.txt. If any of these are unfamiliar, the stops.txt and stop_times.txt relationship guide establishes the join keys this section relies on.
• A chosen metric CRS for the feed’s region. The choosing a projected CRS for transit analysis guide walks through selecting a UTM zone or national grid from the feed’s bounding box.
• Shapely 2.x specifically — the vectorized constructors and STRtree behaviour referenced here differ on the 1.x series.

Concept and Spec Background

How shapes.txt encodes geometry

shapes.txt is a long table, not a geometry column. Each physical path is a set of rows sharing a shape_id; the rows are ordered by an integer shape_pt_sequence that must strictly increase along the path but need not be contiguous. A single LineString is the ordered sequence of (shape_pt_lon, shape_pt_lat) vertices for one shape_id. Routes reference geometry indirectly: trips.txt carries a shape_id, and a route_id typically owns many trips across directions, patterns, and time periods — so one route resolves to several shapes, not one.

Why shape_dist_traveled is not enough

The spec allows shape_dist_traveled on both shapes.txt and stop_times.txt so consumers can locate a stop along a shape without geometry. In practice the column is unreliable: it is optional, agencies populate it in metres, kilometres, miles, or feed-internal units with no declaration, and the values on stop_times.txt frequently do not align with the values on shapes.txt for the same trip. Rather than trust it, the robust approach is to build real geometry, project it, and recompute distance-along-line with shapely linear referencing — which is exactly what the snapping workflow below does.

The projected-CRS requirement

Linear referencing, length, and spatial-join distance predicates all assume planar metric coordinates. Running them on EPSG:4326 does not raise an error; it silently returns garbage, because the numbers are degrees. The rule for this entire section mirrors the coordinate reference systems guidance: load geometry in EPSG:4326, project to a metric CRS for every measurement, and only convert back to WGS84 if you are writing GTFS output.

Step-by-Step Implementation

This section assembles the end-to-end pipeline at a high level; each ### step below has a dedicated deep-dive guide with a single complete, copy-pasteable script. Read the steps here for the shape of the workflow, then follow the linked guide for the runnable implementation and its edge cases.

Step 1: Assemble a LineString per shape and pick a representative per route

The first job is to convert shapes.txt rows into geometry and decide which shape represents a route. Sorting is non-negotiable: an unsorted assembly produces a zig-zag LineString that self-intersects and reports a nonsensical length.

import pandas as pd
import geopandas as gpd
from shapely.geometry import LineString

# Explicit dtypes: shape_id must stay a string, sequence a sortable integer
shapes = pd.read_csv(
    "gtfs/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"],
)

# Order first, then group — never group an unsorted frame
shapes = shapes.sort_values(["shape_id", "shape_pt_sequence"], kind="stable")

lines = (
    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(lines, geometry="geometry", crs="EPSG:4326")

Because one route_id owns many shape_id values, you then join through trips.txt and choose one representative shape per route — commonly the longest, measured after projection. The full script, including the trips join and the longest-shape selection, is in extracting route geometry from GTFS shapes.

Step 2: Snap stops onto the route line with linear referencing

Once you have a projected LineString, shapely gives you two functions that are the heart of every stop-to-shape workflow: LineString.project(point) returns the distance along the line to the point nearest a stop, and LineString.interpolate(distance) returns the coordinate at that distance. Subtracting the interpolated point from the original stop yields the perpendicular off-shape offset — a direct data-quality signal.

# shape_line and stops must already be in the SAME metric CRS
def snap_stop(shape_line, stop_geom):
    dist_along = shape_line.project(stop_geom)          # metres from shape start
    snapped = shape_line.interpolate(dist_along)        # point on the line
    offset = stop_geom.distance(snapped)                # perpendicular metres
    return dist_along, offset

dist_along is your recomputed, trustworthy shape_dist_traveled; offset flags stops that sit implausibly far from their route. The complete per-feed script — including projection, per-shape STRtree lookup, and offset thresholds — is in snapping GTFS stops to shapes in Python.

Step 3: Join stops to polygon layers with geopandas

With stops as a projected GeoDataFrame, gpd.sjoin attaches attributes from any polygon layer — census tracts, fare zones, service districts — using a spatial predicate rather than a shared key. The predicate choice (within for point-in-polygon, intersects for buffered catchments) determines the result set.

# stops and zones must share the SAME projected CRS before sjoin
stops_in_zones = gpd.sjoin(
    stops_projected,
    zones_projected[["zone_name", "geometry"]],
    how="left",
    predicate="within",
)

Buffering each stop into a walk catchment before the join turns a point-in-polygon question into a coverage question. The complete script — projection, buffering, sjoin, and per-zone aggregation — is in geopandas spatial joins for transit stops.

Validation and Verification

Geometry assembly fails quietly, so gate each stage with assertions rather than eyeballing a map:

from shapely.validation import explain_validity

# 1. Every assembled line must have at least two distinct vertices
too_short = shape_lines[shape_lines.geometry.apply(lambda ls: len(ls.coords) < 2)]
assert too_short.empty, f"{len(too_short)} shapes collapsed to a single point"

# 2. Projected route lengths must be plausible (metres, not degrees)
projected = shape_lines.to_crs(32617)  # example UTM zone
lengths = projected.geometry.length
assert lengths.min() > 50, "Sub-50 m route: likely an unsorted or degenerate shape"
assert lengths.max() < 500_000, "500 km route: likely a WGS84 length leaking through"

# 3. Geometry validity (self-intersections from mis-ordering)
invalid = projected[~projected.geometry.is_valid]
for shape_id, geom in zip(invalid.index, invalid.geometry):
    print(shape_id, explain_validity(geom))

The length bounds are the single most effective guard in the whole pipeline: a length in the tens of thousands when you expected metres, or a fraction of one when you expected kilometres, almost always means a CRS or sort mistake upstream.

Failure Modes and Edge Cases

• Unsorted shape rows. shapes.txt is frequently exported grouped by shape but not ordered by shape_pt_sequence, or ordered lexically so that sequence 10 precedes 2. Always cast shape_pt_sequence to integer and sort before assembly, or the LineString will fold back on itself.

• Routes with multiple shapes. A route_id spans inbound and outbound directions plus express and short-turn patterns, each with its own shape_id. There is no single “the” geometry for a route — you pick a representative (longest is a common, defensible default) and document the choice.

• Trusting shape_dist_traveled blindly. The optional distance column may be in miles on one feed and metres on the next, and its stop_times.txt values often disagree with its shapes.txt values. Recompute from projected geometry unless you have verified the unit against known landmarks.

• Stops far from their shape. A large perpendicular offset after snapping usually means the stop is matched to the wrong shape, the shape omits a deviation the stop sits on, or the stop coordinate is simply wrong. Treat offset as a data-quality metric, not noise to discard.

• Feeds without shapes.txt. shapes.txt is optional. When it is missing, you cannot snap or measure on-street distance; you fall back to an ordered straight-line path through the stops of each trip, which is acceptable for a rough map but not for distance metrics.

Performance and Scale Notes

National feeds carry shapes.txt files with tens of millions of point rows — the largest table in the archive. Read only the four columns you need with usecols, enforce the dtype map so pandas does not widen shape_pt_sequence to float64, and prefer shapely 2.x vectorized constructors over per-group Python loops:

import numpy as np
from shapely import linestrings  # shapely 2.x vectorized constructor

shapes = shapes.sort_values(["shape_id", "shape_pt_sequence"], kind="stable")
# Group boundaries without a Python-level groupby.apply
codes, uniques = pd.factorize(shapes["shape_id"], sort=False)
coords = np.column_stack([shapes["shape_pt_lon"], shapes["shape_pt_lat"]])
geoms = linestrings(coords, indices=codes)   # one LineString per group, vectorized
shape_lines = gpd.GeoDataFrame({"shape_id": uniques}, geometry=geoms, crs="EPSG:4326")

For spatial joins against large polygon layers, geopandas builds an R-tree automatically, but touching .sindex once before a batch of joins avoids rebuilding it per call. When geometry is a stage in a larger extract-transform-load run, cache the projected GeoDataFrame to GeoParquet so downstream steps skip re-reading and re-projecting. The broader techniques — chunked reads, category dtypes, and Parquet round-trips — are covered in memory-efficient processing for large feeds.

Frequently Asked Questions

Related

• Snapping GTFS Stops to Shapes in Python — linear referencing to compute distance-along-shape and off-shape offset per stop
• Extracting Route Geometry from GTFS Shapes — assemble LineStrings and choose one representative geometry per route
• GeoPandas Spatial Joins for Transit Stops — attach polygon attributes and aggregate stops by zone
• Coordinate Reference Systems for Transit Data — the projection rules every measurement in this section depends on
• Mastering stops.txt and stop_times.txt Relationships — the join keys linking geometry back to schedule data

Up: GTFS Feed Architecture & Fundamentals | Home