Choosing a Projected CRS for Transit

For a feed that fits within a single UTM zone — a city or a metropolitan network — derive the zone number from the feed’s centroid longitude and use that UTM CRS, which gives you metre units with negligible distortion. For a feed that spans several zones, such as a national rail or intercity bus network, a single UTM zone distorts distance and area toward its edges, so switch to a regional equal-area coordinate reference system instead. The choice follows mechanically from the feed’s bounding box.

Root Cause Analysis

Every metric spatial operation on a GTFS feed — a walking-catchment buffer, a stop-density surface, a nearest-stop join against parcels — needs coordinates in metres. WGS84 delivers degrees, and a degree is not a fixed distance: it shrinks from about 111 km of east-west extent at the equator toward zero at the poles. Choosing a projected CRS is the act of picking how to flatten the curved surface into metres, and no flat projection preserves every property at once. You trade among conformality, equal area, and true distance, and the right trade depends entirely on the geographic extent of the feed.

Universal Transverse Mercator solves the local case elegantly. It divides the globe into sixty six-degree-wide zones, each with its own conformal projection tuned to keep distortion under about one part in a thousand within the zone. For a feed that lives inside one zone, UTM is close to ideal: shapes are preserved, distances are accurate, and the units are metres. The catch is the six-degree width. Push a feed past the zone boundary and the scale factor grows quickly; by two or three zones out, buffer radii and area calculations are visibly wrong at the edges. A national network draped across a single UTM zone will report stop densities that sag on one side of the country and swell on the other.

That is where equal-area projections earn their place. A continental equal-area CRS — Lambert Azimuthal Equal-Area for Europe, Albers Equal-Area for the conterminous United States, and similar national grids elsewhere — sacrifices local shape fidelity to guarantee that a square kilometre is a square kilometre everywhere in its domain. For any analysis that sums or compares area or density across a wide feed, that guarantee matters more than the slight shape distortion. The decision therefore reduces to one measurable question: does the feed fit inside a single UTM zone, or not? This page automates that question from the feed’s bounding box, and the same reasoning underpins choosing coordinate reference systems for transit data generally, including the spatial-analysis workflows in the route-geometry guide.

Production-Ready Python Implementation

The script reads stops.txt, computes the bounding box, measures how many UTM zones the longitude span crosses, and returns a CRS decision: a centroid-derived UTM EPSG for single-zone feeds, or a configured regional equal-area CRS for multi-zone feeds. It then confirms the chosen CRS reports metre units before any analysis proceeds. It uses geopandas and pyproj.

import logging
import math
from pathlib import Path

import pandas as pd
import geopandas as gpd
from pyproj import CRS

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

# Regional equal-area fallbacks, keyed by an approximate longitude window.
# Extend this table with national grids relevant to your feeds.
EQUAL_AREA_BY_REGION = {
    "europe": (-31.0, 45.0, 3035),      # ETRS89-LAEA Europe
    "conus": (-125.0, -66.0, 5070),     # NAD83 Albers, conterminous US
}
GLOBAL_EQUAL_AREA = 6933                # World Cylindrical Equal Area fallback


def _utm_zone(lon: float) -> int:
    """UTM zone number (1-60) containing a longitude in degrees."""
    return int(math.floor((lon + 180.0) / 6.0) + 1)


def _utm_epsg(lon: float, lat: float) -> int:
    """EPSG code for the WGS84 UTM zone containing (lon, lat)."""
    zone = _utm_zone(lon)
    return (32600 if lat >= 0 else 32700) + zone


def _regional_equal_area(centroid_lon: float) -> int:
    """Pick a regional equal-area EPSG by centroid longitude, else global."""
    for name, (lon_min, lon_max, epsg) in EQUAL_AREA_BY_REGION.items():
        if lon_min <= centroid_lon <= lon_max:
            logging.info("Selected %s equal-area CRS EPSG:%d", name, epsg)
            return epsg
    logging.info("No regional match; using global equal-area EPSG:%d", GLOBAL_EQUAL_AREA)
    return GLOBAL_EQUAL_AREA


def choose_crs(feed_dir: str) -> CRS:
    """
    Choose a projected CRS from a GTFS feed's bounding box.

    Single-zone feeds get the centroid UTM zone; multi-zone feeds get a
    regional (or global) equal-area CRS. The returned CRS always uses metres.
    """
    feed_path = Path(feed_dir)
    stops = pd.read_csv(
        feed_path / "stops.txt",
        dtype={"stop_id": str, "stop_lat": float, "stop_lon": float},
        usecols=["stop_id", "stop_lat", "stop_lon"],
    )

    # Validate WGS84 bounds before trusting the extent.
    valid = stops["stop_lat"].between(-90, 90) & stops["stop_lon"].between(-180, 180)
    if not valid.all():
        raise ValueError(f"{(~valid).sum()} stops fall outside WGS84 bounds")

    lon_min, lon_max = stops["stop_lon"].min(), stops["stop_lon"].max()
    centroid_lon = float(stops["stop_lon"].mean())
    centroid_lat = float(stops["stop_lat"].mean())

    zone_span = _utm_zone(lon_max) - _utm_zone(lon_min) + 1
    logging.info(
        "Feed spans %d UTM zone(s): lon %.4f to %.4f", zone_span, lon_min, lon_max
    )

    if zone_span <= 1:
        epsg = _utm_epsg(centroid_lon, centroid_lat)
        logging.info("Single-zone feed -> UTM EPSG:%d", epsg)
    else:
        epsg = _regional_equal_area(centroid_lon)

    crs = CRS.from_epsg(epsg)

    # Guard: the chosen CRS must be projected and metric.
    unit = crs.axis_info[0].unit_name
    if crs.is_geographic or unit != "metre":
        raise ValueError(f"Chosen CRS EPSG:{epsg} is not metre-based (unit={unit})")

    logging.info("Chosen CRS: EPSG:%d — %s (%s)", epsg, crs.name, unit)
    return crs


if __name__ == "__main__":
    crs = choose_crs("/data/gtfs/national_rail")
    # Reproject stops into the chosen CRS for downstream metric analysis.
    stops = pd.read_csv(
        "/data/gtfs/national_rail/stops.txt",
        dtype={"stop_id": str, "stop_lat": float, "stop_lon": float},
    )
    gdf = gpd.GeoDataFrame(
        stops,
        geometry=gpd.points_from_xy(stops["stop_lon"], stops["stop_lat"]),
        crs="EPSG:4326",
    ).to_crs(crs)
    print(f"Projected {len(gdf)} stops into {crs.name}")

Step-by-Step Walkthrough

Bounds validation comes first. A single out-of-range coordinate — a stop accidentally placed at longitude 0, 0 or with swapped lat/lon — would blow up the bounding box and push the zone-span count to a nonsensical value. Rejecting invalid coordinates before measuring extent keeps the CRS decision trustworthy.

Zone span, not raw width, drives the branch. _utm_zone converts each edge longitude to a zone number, and the count of zones crossed is what decides the branch. A metro that straddles a zone boundary by a few kilometres still counts as spanning two zones, and the code correctly prefers an equal-area CRS in that borderline case rather than picking a UTM zone that clips half the network.

UTM EPSG is computed, not looked up. For single-zone feeds, _utm_epsg builds the EPSG from the centroid: 32600 plus the zone number in the northern hemisphere, 32700 plus the zone in the southern. This works for any region on Earth without a hard-coded table, so the same function serves feeds from Berlin, Santiago, or Perth.

Equal-area selection is table-driven. For multi-zone feeds, _regional_equal_area matches the centroid longitude against a small table of continental grids and falls back to a global equal-area projection when nothing matches. Extending the table with a national grid — a country-specific Albers or Lambert definition — is the intended way to tune the pipeline for your feeds.

The metre-unit guard is non-negotiable. Before returning, the code asserts the chosen CRS is projected and reports metre as its axis unit. This catches the classic mistake of accidentally selecting a geographic CRS, which would silently make every downstream buffer radius meaningless. Getting the units right here is the same discipline that keeps schedule times correct when converting to UTC — the spatial and temporal layers each have one canonical unit that must be verified, never assumed.

Verification and Output

Prove that the chosen CRS actually produces sane metric coordinates for the feed.

import geopandas as gpd
import pandas as pd

crs = choose_crs("/data/gtfs/national_rail")

stops = pd.read_csv(
    "/data/gtfs/national_rail/stops.txt",
    dtype={"stop_id": str, "stop_lat": float, "stop_lon": float},
)
gdf = gpd.GeoDataFrame(
    stops,
    geometry=gpd.points_from_xy(stops["stop_lon"], stops["stop_lat"]),
    crs="EPSG:4326",
).to_crs(crs)

# Projected coordinates must be finite and metric-scaled, not tiny degree values
xs, ys = gdf.geometry.x, gdf.geometry.y
assert xs.notna().all() and ys.notna().all(), "Projection produced NaN coordinates"
assert xs.abs().max() > 1000, "Coordinates look like degrees, not metres — wrong CRS"

# A 400 m buffer should have an area near pi * 400^2 square metres
sample_area = gdf.geometry.iloc[0].buffer(400).area
assert 490_000 < sample_area < 510_000, f"Buffer area {sample_area:.0f} m^2 is distorted"

print(f"CRS EPSG:{crs.to_epsg()} verified: {len(gdf)} stops projected to metres")

The buffer-area check is the decisive test: a 400-metre buffer must enclose roughly 502,655 square metres regardless of where the feed sits. If a feed you expected to be single-zone trips the distortion assertion, its stops span more longitude than you thought, and the equal-area branch is the correct answer. The function returns a pyproj.CRS object you pass straight to GeoDataFrame.to_crs, so the decision and the projection stay in one place.

Gotchas and Edge Cases

• Feeds crossing the antimeridian. A network spanning the 180th meridian (parts of the Pacific, the Russian Far East) produces a bounding box that appears to wrap the entire globe, inflating the zone span to sixty. Detect a longitude range greater than 180 degrees and shift coordinates into a continuous frame before measuring extent, or the equal-area fallback will trigger spuriously.

• A handful of far-flung stops. One mislocated stop — an island depot or a test coordinate left in the feed — can drag the bounding box across several zones and force an equal-area CRS on an otherwise local network. Consider using a robust extent (for example the 1st–99th longitude percentile) rather than the raw min and max, or clean outliers first with the validation guide.

• National grids beat generic UTM at the edges. Where a country publishes an official grid — British National Grid, RGF93 Lambert-93, and their peers — it is usually tuned better than the raw UTM zone for that territory. If you know the feed’s country, prefer its national grid over the computed UTM zone by adding it to the region table.

• Always revert to WGS84 for output. A projected CRS is for analysis only. Any stops.txt or shapes.txt you write back must be in EPSG:4326, because trip planners and routing engines reject anything else. Keep projection confined to an isolated analysis layer and reproject to WGS84 before export.

Frequently Asked Questions

Related

• Coordinate Reference Systems for Transit Data — validating, projecting, and round-tripping GTFS coordinates end to end
• Spatial Analysis and Route Geometry — the buffer, join, and geometry workflows that consume the CRS chosen here
• Converting Local Transit Times to UTC in Python — the temporal counterpart to picking the right spatial reference frame
• Up: Coordinate Reference Systems for Transit Data — the spec background and projection round-trip this decision fits into
• Section: GTFS Feed Architecture & Fundamentals — the schema foundations beneath every spatial workflow · Home