Coordinate Reference Systems for Transit Data

Misaligned coordinate reference system (CRS) assumptions are one of the most common sources of silent spatial corruption in GTFS pipelines. The GTFS Feed Architecture & Fundamentals pillar treats the schema constraints that govern every feed file; this page goes deeper on the spatial layer: how to validate that stops.txt coordinates genuinely conform to WGS84, how to project them into metric planar systems for analysis, and how to round-trip them back without detectable drift. The workflow here is production-ready for transit analysts and Python GIS engineers dealing with feeds from a single agency up to multi-region batches.

The Spatial Problem in GTFS Pipelines

GTFS mandates a single geographic coordinate system for data exchange, but real-world transit analytics — walking catchment buffers, spatial joins with zoning or demographic layers, network distance routing — require projected metric coordinates. The gap between these two needs introduces risk at every stage: unvalidated inputs may already contain out-of-bounds values, projection grid files may be missing from the processing environment, and round-trip floating-point loss can shift stop positions by meters, enough to break routing engines that snap stops to road segments.

The challenge is compounded by the fact that GTFS schema validation tools such as the MobilityData canonical validator check column presence and type but do not verify that stop_lat/stop_lon values actually fall within WGS84 bounds or that they have sufficient decimal precision. Spatial errors therefore survive schema checks and propagate silently into production.

Prerequisites

Before running the code in this guide, ensure the following are in place:

• Python 3.9+ with pandas, geopandas, pyproj, shapely, and numpy installed:

pip install pandas geopandas pyproj shapely numpy

• A validated GTFS static feed with stops.txt, stop_times.txt, and optionally shapes.txt. If you are building or auditing a pipeline from scratch, review Understanding GTFS Static Feed Structure to confirm your ingestion layer aligns with specification requirements before applying spatial transformations.
• Access to PROJ’s datum-shift grids (proj-data package) if transforming between non-WGS84 datums.
• Knowledge of the geographic region covered by the feed, so you can select an appropriate projected CRS.

Concept and Spec Background

WGS84 and EPSG:4326

The GTFS specification requires all latitude and longitude values in stops.txt — and in shapes.txt for route geometry — to be expressed in WGS84 (EPSG:4326). WGS84 is a geographic (angular) coordinate system: coordinates are in decimal degrees, axes are latitude/longitude, and the reference ellipsoid is the 1984 World Geodetic System. It is globally consistent and universally supported by web mapping platforms, trip planners, and routing engines.

The trade-off is geometric distortion. Euclidean distance calculations in degree-based coordinates produce meaningless results at anything but equatorial latitudes. A 400-metre pedestrian buffer around a stop cannot be expressed as a degree offset without knowing the local scale factor — which varies by latitude and projection.

When to Use Projected Coordinates

Metric-based projected CRS options for transit work include:

The governing principle: project only within an isolated analysis layer; always revert to EPSG:4326 before writing GTFS output.

GTFS Files Involved

Step-by-Step Implementation

Step 1: Validate Input Coordinates and Build a GeoDataFrame

Silent failures in transit pipelines often originate from unvalidated input. Before projecting, verify that your stops.txt coordinates fall within WGS84 bounds and that the GeoDataFrame is constructed with an explicit CRS declaration.

import pandas as pd
import numpy as np
import geopandas as gpd
from shapely.geometry import Point

# Load stops with explicit dtypes to prevent silent coercion
stops = pd.read_csv(
    "gtfs/stops.txt",
    dtype={
        "stop_id": str,
        "stop_name": str,
        "stop_lat": float,
        "stop_lon": float,
        "location_type": "Int64",  # nullable integer (optional column)
        "parent_station": str,
    },
)

# Validate WGS84 coordinate bounds
lat_ok = stops["stop_lat"].between(-90.0, 90.0)
lon_ok = stops["stop_lon"].between(-180.0, 180.0)
invalid = stops[~(lat_ok & lon_ok)]

if not invalid.empty:
    raise ValueError(
        f"{len(invalid)} stops have coordinates outside WGS84 bounds:\n"
        f"{invalid[['stop_id', 'stop_lat', 'stop_lon']].to_string()}"
    )

# Require at least 5 decimal places (~1.1 m precision at equator)
insufficient_precision = stops[
    (stops["stop_lat"].round(5) == stops["stop_lat"].round(2)) |
    (stops["stop_lon"].round(5) == stops["stop_lon"].round(2))
]
if not insufficient_precision.empty:
    print(
        f"WARNING: {len(insufficient_precision)} stops have fewer than 5 decimal "
        "places of precision — trip planners may mis-snap them to road segments."
    )

# Build GeoDataFrame with CRS declared at construction time.
# Explicit declaration prevents pyproj from guessing axis order.
gdf_stops = gpd.GeoDataFrame(
    stops,
    geometry=gpd.points_from_xy(stops["stop_lon"], stops["stop_lat"]),
    crs="EPSG:4326",
)

print(f"Loaded {len(gdf_stops)} stops | CRS: {gdf_stops.crs.to_epsg()}")

The explicit dtype map prevents stop_lat and stop_lon being read as strings in feeds where missing values appear as empty fields — a common real-agency quirk that causes all downstream floating-point operations to fail silently with NaN propagation.

Step 2: Select and Apply a Projected CRS

Use the geographic centroid of the feed to derive the appropriate UTM zone automatically, rather than hard-coding a zone that may be wrong for feeds from other regions.

from pyproj import CRS
import math

def utm_epsg_from_lonlat(lon: float, lat: float) -> int:
    """Return the EPSG code for the UTM zone containing (lon, lat)."""
    zone_number = math.floor((lon + 180) / 6) + 1
    if lat >= 0:
        return 32600 + zone_number   # Northern hemisphere
    return 32700 + zone_number       # Southern hemisphere

# Derive centroid from the feed's bounding box
centroid_lon = float(gdf_stops["stop_lon"].mean())
centroid_lat = float(gdf_stops["stop_lat"].mean())
utm_epsg = utm_epsg_from_lonlat(centroid_lon, centroid_lat)

target_crs = CRS.from_epsg(utm_epsg)
print(f"Selected CRS: EPSG:{utm_epsg} — {target_crs.name}")
print(f"Linear unit: {target_crs.axis_info[0].unit_name}")  # must be 'metre'

# Project — geopandas internally uses pyproj with always_xy=True since v0.7
gdf_projected = gdf_stops.to_crs(target_crs)

# Sanity-check: projected coordinates should be in the hundreds-of-thousands range
x_range = gdf_projected.geometry.x.agg(["min", "max"])
y_range = gdf_projected.geometry.y.agg(["min", "max"])
print(f"Easting range:  {x_range['min']:.0f} – {x_range['max']:.0f} m")
print(f"Northing range: {y_range['min']:.0f} – {y_range['max']:.0f} m")

Always confirm that axis_info[0].unit_name returns 'metre' before any buffer or distance operation. If the unit is 'degree', the CRS declaration was not applied correctly.

Step 3: Spatial Operations in Projected Space

With metric coordinates in place, buffer operations and spatial joins produce geometrically accurate results. A typical workflow combines walking catchments with stop schedule frequency data from stop_times.txt.

# Create 400 m walking-catchment buffers around each stop
gdf_buffers = gdf_projected.copy()
gdf_buffers["geometry"] = gdf_buffers["geometry"].buffer(400)

# Load stop_times with explicit dtypes
stop_times = pd.read_csv(
    "gtfs/stop_times.txt",
    dtype={
        "trip_id": str,
        "stop_id": str,
        "stop_sequence": int,
        "arrival_time": str,    # kept as str; GTFS times can exceed 24:00
        "departure_time": str,
    },
    usecols=["trip_id", "stop_id", "stop_sequence"],
)

# Count distinct trips serving each stop
trips_per_stop = (
    stop_times
    .groupby("stop_id", sort=False)["trip_id"]
    .nunique()
    .reset_index(name="trip_count")
)

# Join frequency data to spatial buffers
gdf_analysis = gdf_buffers.merge(
    trips_per_stop,
    on="stop_id",
    how="left",
)
gdf_analysis["trip_count"] = gdf_analysis["trip_count"].fillna(0).astype(int)

# Referential integrity check: every stop_id in stop_times must exist in stops
orphaned_stops = set(stop_times["stop_id"]) - set(gdf_stops["stop_id"])
if orphaned_stops:
    print(f"WARNING: {len(orphaned_stops)} stop_ids in stop_times.txt have no "
          f"matching entry in stops.txt: {list(orphaned_stops)[:5]} …")

Spatial joins and attribute merges must respect primary key relationships. Misaligned stop identifiers between spatial and temporal tables silently drop records rather than raising errors — referential integrity checks like the one above should be part of every pipeline run. For schema-level validation of key relationships across the entire feed, see GTFS Validation Rules and Common Schema Errors.

Step 4: Round-Trip Reversion and Drift Measurement

After analysis, revert all geometries to WGS84. GTFS consumers — trip planners, OpenTripPlanner, Valhalla, and mobility-as-a-service APIs — expect strict EPSG:4326 in every feed they ingest.

# Revert to WGS84
gdf_export = gdf_analysis.to_crs("EPSG:4326")

# Extract coordinates back into DataFrame columns
gdf_export = gdf_export.copy()
gdf_export["stop_lon"] = gdf_export.geometry.x
gdf_export["stop_lat"] = gdf_export.geometry.y

# Measure round-trip drift against original values
original = stops.set_index("stop_id")[["stop_lat", "stop_lon"]]
exported = gdf_export.set_index("stop_id")[["stop_lat", "stop_lon"]]
drift = np.sqrt(
    (exported["stop_lon"] - original["stop_lon"]) ** 2 +
    (exported["stop_lat"] - original["stop_lat"]) ** 2
)

max_drift = drift.max()
p99_drift = drift.quantile(0.99)
print(f"Max coordinate drift: {max_drift:.10f}°  ({max_drift * 111_320:.4f} m)")
print(f"P99 coordinate drift: {p99_drift:.10f}°  ({p99_drift * 111_320:.4f} m)")

# Acceptable threshold: < 1e-7° (~1 cm)
assert max_drift < 1e-7, (
    f"Round-trip drift {max_drift:.2e}° exceeds 1e-7° threshold. "
    "Check projection grid files and datum alignment."
)

# Build the GTFS-compliant export DataFrame
export_df = (
    gdf_export
    .drop(columns="geometry")
    .rename(columns={})   # column names already match stops.txt spec
)

The threshold of 1e-7 degrees corresponds to roughly 1 centimetre at mid-latitudes — well within the precision that any downstream routing engine can act on. Drift larger than 1e-5 degrees (~1 metre) indicates a datum shift file is missing or was silently bypassed.

Validation and Verification

Run the following checks as a gate in your CI pipeline after every feed ingest:

def validate_crs_pipeline(
    original_stops: pd.DataFrame,
    exported_stops: pd.DataFrame,
    max_drift_degrees: float = 1e-7,
) -> dict:
    """
    Compare original and exported stop coordinates.
    Returns a dict with pass/fail status and diagnostic metrics.
    """
    orig = original_stops.set_index("stop_id")[["stop_lat", "stop_lon"]]
    expo = exported_stops.set_index("stop_id")[["stop_lat", "stop_lon"]]

    # Ensure same set of stop_ids
    assert set(orig.index) == set(expo.index), "stop_id sets differ between input and output"

    drift = np.sqrt(
        (expo["stop_lon"] - orig["stop_lon"]) ** 2 +
        (expo["stop_lat"] - orig["stop_lat"]) ** 2
    )

    results = {
        "stop_count": len(orig),
        "max_drift_deg": float(drift.max()),
        "max_drift_cm": float(drift.max() * 111_320 * 100),
        "p99_drift_deg": float(drift.quantile(0.99)),
        "passes_threshold": bool(drift.max() < max_drift_degrees),
        "worst_stop_id": str(drift.idxmax()),
    }

    assert results["passes_threshold"], (
        f"CRS round-trip drift {results['max_drift_deg']:.2e}° at stop "
        f"{results['worst_stop_id']} exceeds {max_drift_degrees:.0e}° threshold."
    )

    return results

Failure Modes and Edge Cases

• Axis-order swaps (lat/lon vs lon/lat). Older versions of pyproj (before 2.2) defaulted to lat, lon axis ordering for geographic CRS. When constructing a Transformer manually, always pass always_xy=True. Using GeoDataFrame.to_crs() via geopandas 0.7+ applies always_xy=True internally, but any raw pyproj.Transformer call in surrounding code may not.

• Missing datum-shift grid files. Transforming between datums (NAD27 → WGS84, OSGB36 → WGS84, etc.) requires PROJ grid shift files (.tif via the proj-data package). Without them, PROJ falls back to a 3-parameter Helmert transform that can introduce errors of 10–100 metres. Check for proj-data in your environment and validate against a known benchmark coordinate.

• Feeds that cross UTM zone boundaries. A single feed covering a large region (e.g. a national rail network or a multi-state transit authority) may span multiple UTM zones. Projecting the entire feed into a single UTM zone introduces scale distortion toward the edges. Use a continental equal-area projection (e.g. EPSG:3035 for Europe, EPSG:5070 for the conterminous US) for whole-network analysis.

• Coordinate clipping by routing engines. Some trip planners and GTFS importers enforce strict bounding boxes. After round-trip export, verify that no stop coordinate has been clamped or rounded by intermediate serialisation. A stop at stop_lat=47.000000 where the original was 47.000012 indicates somewhere a float was rounded to 6 decimal places — safe for most consumers but worth logging.

• Timezone vs. CRS confusion in combined pipelines. Timezone handling and schedule normalisation and CRS projection are independent operations. Schedule timestamps do not carry spatial information; coordinate transformations do not alter clock times. Pipelines that interleave both operations in the same transform step are a maintenance hazard — keep them in separate, testable functions.

• stop_lat/stop_lon written in wrong column order. When extracting projected geometry back to columns, geometry.x is longitude (easting) and geometry.y is latitude (northing) — the reverse of the GTFS column order. A one-line swap (stop_lat = geometry.y, stop_lon = geometry.x) in the wrong direction will silently produce feeds with coordinates on the wrong continent.

Performance and Scale Notes

For feeds above 500 MB (large regional or national operators), the following strategies reduce memory pressure:

import pyarrow as pa
import pyarrow.parquet as pq

# Chunked read for large stops.txt (rare but possible in aggregated feeds)
chunk_iter = pd.read_csv(
    "gtfs/stops.txt",
    dtype={"stop_id": str, "stop_lat": float, "stop_lon": float},
    chunksize=50_000,
)

validated_chunks = []
for chunk in chunk_iter:
    lat_ok = chunk["stop_lat"].between(-90.0, 90.0)
    lon_ok = chunk["stop_lon"].between(-180.0, 180.0)
    if not (lat_ok.all() and lon_ok.all()):
        raise ValueError("Out-of-bounds coordinates detected in chunk")
    validated_chunks.append(chunk)

stops_full = pd.concat(validated_chunks, ignore_index=True)

# After projection, cache the GeoDataFrame as GeoParquet for reuse
# (avoids re-reading and re-projecting on every pipeline run)
gdf_projected.to_parquet("cache/stops_projected.parquet")

# Reload without re-projecting:
# gdf_projected = gpd.read_parquet("cache/stops_projected.parquet")

For spatial joins against large polygon layers (parcel data, census tracts), use a spatial index. geopandas builds an R-tree index automatically when you call sjoin, but explicitly calling .sindex in advance avoids rebuilding it for repeated joins in a batch pipeline:

# Pre-build spatial index before batched sjoin operations
_ = gdf_projected.sindex   # triggers index construction and caches it

# Now batch sjoin calls reuse the same index
for zone_gdf in zone_batch:
    result = gpd.sjoin(gdf_projected, zone_gdf, how="left", predicate="within")

For memory-efficient processing of large GTFS feeds more broadly — including chunked reading of stop_times.txt and trips.txt — see the dedicated guide on that topic.

Related

• Mastering stops.txt and stop_times.txt Relationships — referential integrity and join patterns between the spatial and schedule layers
• GTFS Validation Rules and Common Schema Errors — schema-level checks that run before spatial validation
• Timezone Handling and Schedule Normalisation — the temporal counterpart to spatial projection in GTFS pipelines
• Memory-Efficient Processing for Large Feeds — chunked reading and Parquet caching for feeds that exceed available RAM
• Understanding GTFS Static Feed Structure — the full file-by-file schema context before applying spatial transformations

Up: GTFS Feed Architecture & Fundamentals | Home