1. Home
  2. GTFS Feed Architecture Fundamentals
  3. Coordinate Reference Systems For Transit Data
  4. Reprojecting GTFS Shapes With Pyproj

Reprojecting GTFS Shapes with pyproj

Compute the feed’s own bounding box, derive the UTM zone from its centroid, build one Transformer with always_xy=True, and transform whole coordinate arrays in a single call rather than looping over points. Then verify the round trip returns the original coordinates before trusting anything measured in the projected space. The two errors that make everything downstream wrong are omitting always_xy — which silently swaps latitude and longitude — and reaching for Web Mercator, which distorts distance by the secant of the latitude. Choosing the projection itself is covered in coordinate reference systems for transit data.

East-west distance error from measuring in degrees How far a degree of longitude falls short of a degree of latitude, by latitude — the factor by which unprojected east-west distances are wrong. Latitude 60° 50 % Oslo, Helsinki, Saint Petersburg Latitude 50° 36 % London, Vancouver, Prague Latitude 40° 23 % New York, Madrid, Beijing Latitude 20° 6 % Mexico City, Mumbai Latitude 0° 0 % the only place degrees behave the same code, the same feed, a different city

Root Cause Analysis

GTFS stores coordinates in WGS84 decimal degrees, which is the correct choice for exchange and the wrong unit for every measurement anyone wants to make. A degree of latitude is about 111 km everywhere; a degree of longitude is 111 km at the equator and 71 km at 50 degrees north. Treating the pair as a planar coordinate system means east–west distances are wrong by a factor that depends on where the agency happens to be, which is why the same snapping code behaves acceptably in Quito and badly in Helsinki.

Reprojection fixes it, and introduces three failure modes of its own.

Axis order. EPSG:4326 defines its axes as latitude, then longitude. Almost every piece of software, and every GTFS column pair, works in longitude, latitude order. pyproj respects the authority definition unless told otherwise, so a transformer built without always_xy=True interprets the first value as latitude. The result is not an error — it is a set of coordinates in the wrong place, and if the feed happens to sit near the diagonal it can look almost plausible.

Per-point transformation. Transformer.transform accepts arrays. Called inside a loop over 214,000 shape vertices it does 214,000 individual transformations with the projection pipeline overhead paid each time. The array form does the same work in one call, roughly a hundred times faster.

Rebuilding the transformer. Constructing a Transformer involves resolving the CRS definitions and building a pipeline, which is expensive relative to the transformation itself. Building one inside a per-route or per-trip function is a common and very effective way to make a spatial pipeline slow.

Production-Ready Python Implementation

python 
"""Reproject GTFS coordinates into a metric CRS with pyproj."""
from __future__ import annotations

import logging
from dataclasses import dataclass
from functools import lru_cache

import numpy as np
import pandas as pd
from pyproj import CRS, Transformer

log = logging.getLogger("gtfs.reproject")

WGS84 = "EPSG:4326"
# Beyond this longitudinal extent a single UTM zone distorts too much to be honest.
MAX_UTM_SPAN_DEGREES = 6.0


@dataclass(frozen=True)
class Extent:
    min_lon: float
    min_lat: float
    max_lon: float
    max_lat: float

    @property
    def centroid(self) -> tuple[float, float]:
        return ((self.min_lon + self.max_lon) / 2, (self.min_lat + self.max_lat) / 2)

    @property
    def lon_span(self) -> float:
        return self.max_lon - self.min_lon


def feed_extent(*frames: tuple[pd.Series, pd.Series]) -> Extent:
    """Bounding box over any number of (lon, lat) column pairs."""
    lons = pd.concat([lon.astype("float64") for lon, _ in frames])
    lats = pd.concat([lat.astype("float64") for _, lat in frames])
    lons, lats = lons.dropna(), lats.dropna()
    if lons.empty or lats.empty:
        raise ValueError("no usable coordinates in this feed")
    return Extent(float(lons.min()), float(lats.min()),
                  float(lons.max()), float(lats.max()))


def utm_epsg(extent: Extent) -> str:
    """The UTM zone containing the feed's centroid."""
    lon, lat = extent.centroid
    zone = int((lon + 180) // 6) + 1
    # 326xx is northern hemisphere, 327xx southern.
    code = 32600 + zone if lat >= 0 else 32700 + zone
    if extent.lon_span > MAX_UTM_SPAN_DEGREES:
        log.warning("feed spans %.1f degrees of longitude, more than one UTM zone — "
                    "EPSG:%d will distort at the edges; consider a national grid or "
                    "an equal-area projection", extent.lon_span, code)
    return f"EPSG:{code}"


@lru_cache(maxsize=32)
def _transformer(source: str, target: str) -> Transformer:
    """Cached: building a Transformer is far more expensive than using one."""
    # always_xy keeps coordinates in (lon, lat) order, which is what every GTFS
    # column pair uses. Without it, EPSG:4326's authority order silently swaps them.
    return Transformer.from_crs(source, target, always_xy=True)


def project(lon: pd.Series, lat: pd.Series, target: str,
            source: str = WGS84) -> tuple[np.ndarray, np.ndarray]:
    """Transform whole arrays at once — never point by point."""
    transformer = _transformer(source, target)
    x, y = transformer.transform(lon.astype("float64").to_numpy(),
                                 lat.astype("float64").to_numpy())
    return x, y


def project_frame(frame: pd.DataFrame, lon_column: str, lat_column: str,
                  target: str) -> pd.DataFrame:
    out = frame.copy()
    out["x"], out["y"] = project(frame[lon_column], frame[lat_column], target)
    return out


def choose_and_project(stops: pd.DataFrame,
                       shapes: pd.DataFrame | None = None
                       ) -> tuple[str, pd.DataFrame, pd.DataFrame | None]:
    pairs = [(stops["stop_lon"], stops["stop_lat"])]
    if shapes is not None and not shapes.empty:
        pairs.append((shapes["shape_pt_lon"], shapes["shape_pt_lat"]))

    extent = feed_extent(*pairs)
    target = utm_epsg(extent)
    name = CRS.from_user_input(target).name
    log.info("feed extent %.3f,%.3f to %.3f,%.3f — projecting to %s (%s)",
             extent.min_lon, extent.min_lat, extent.max_lon, extent.max_lat,
             target, name)

    projected_stops = project_frame(stops, "stop_lon", "stop_lat", target)
    projected_shapes = (project_frame(shapes, "shape_pt_lon", "shape_pt_lat", target)
                        if shapes is not None and not shapes.empty else None)
    return target, projected_stops, projected_shapes
Reprojecting once, at the boundary Measure the feed's own extent, derive the zone from its centroid, build one cached transformer, and transform whole arrays. Extent from the coordinates Zone UTM from the centroid Transformer built once, cached Transform arrays, not points always_xy keeps longitude first, which is the order every GTFS column pair uses

Step-by-Step Walkthrough

The zone comes from the feed, not from configuration. utm_epsg derives it from the centroid of the actual coordinates, so the same code works for a feed in Boston, Bogotá or Brisbane without anyone editing a constant. The 32600/32700 offsets encode the hemisphere, which is why the latitude sign is checked rather than assumed.

A feed too wide for one zone gets a warning, not an error. UTM zones are six degrees wide, and a regional feed can exceed that. The warning names the alternative — a national grid or an equal-area projection — rather than silently picking one, because that decision depends on what the measurements are for.

always_xy=True carries a comment explaining what it prevents. This is the single most consequential argument in the module. It is also invisible in its effect: without it nothing raises, the coordinates are simply somewhere else.

_transformer is lru_cached on the CRS pair. Every call site can ask for a transformer without thinking about reuse, and the expensive construction happens once per pair. A cache size of 32 covers any realistic set of source and target combinations in one process.

transform receives NumPy arrays. Series.to_numpy() hands pyproj a contiguous block it can transform in one pass. This is where the two-orders-of-magnitude difference lives, and it is why project takes Series rather than scalars — the signature makes the fast path the obvious one.

project_frame copies rather than mutating. The original degree columns are kept alongside the projected ones, so anything that needs to write coordinates back out — a PostGIS load, a re-export — still has the values as published.

Verification and Output

python 
def verify_round_trip(lon: pd.Series, lat: pd.Series, target: str,
                      tolerance_m: float = 0.001) -> None:
    """Project and unproject; the result must land back where it started."""
    x, y = project(lon, lat, target)
    back_lon, back_lat = _transformer(target, WGS84).transform(x, y)

    # Compare in metres, not degrees, so the tolerance means the same everywhere.
    dlat_m = (np.asarray(back_lat) - lat.to_numpy()) * 111_320.0
    dlon_m = ((np.asarray(back_lon) - lon.to_numpy()) * 111_320.0
              * np.cos(np.radians(lat.to_numpy())))
    worst = float(np.max(np.hypot(dlat_m, dlon_m)))
    assert worst < tolerance_m, (
        f"round trip moved a point by {worst * 1000:.1f} mm — check always_xy and "
        "that the source CRS really is WGS84")


def verify_extent(projected: pd.DataFrame, target: str) -> None:
    x, y = projected["x"], projected["y"]
    # A UTM easting is always within 500 km of the zone's central meridian.
    assert (x.between(-500_000, 1_500_000)).all(), (
        "projected easting is outside any plausible UTM range — coordinates were "
        "probably passed in the wrong axis order")
    assert (y.between(-10_000_000, 10_000_000)).all(), "implausible northing"
    span_km = float(np.hypot(x.max() - x.min(), y.max() - y.min())) / 1000
    log.info("projected extent spans %.1f km diagonally in %s", span_km, target)

The easting range check is a cheap, decisive test for the axis-order bug. A UTM easting is bounded by construction; a latitude passed in as a longitude produces a value far outside it, so the assertion fires immediately rather than letting a mirrored network reach a measurement.

Normal output:

text 
INFO gtfs.reproject: feed extent -71.191,42.227 to -70.986,42.450 — projecting to EPSG:32619 (WGS 84 / UTM zone 19N)
INFO gtfs.reproject: projected extent spans 30.4 km diagonally in EPSG:32619

A regional feed that has outgrown UTM:

text 
WARN gtfs.reproject: feed spans 9.4 degrees of longitude, more than one UTM zone — EPSG:32632 will distort at the edges; consider a national grid or an equal-area projection
When UTM stops being the right answer Three extents and the projection each calls for; the warning fires when a feed outgrows a single zone. How wide is the feed's extent? under 6° of longitude One UTM zone metres, minimal distortion one country A national grid designed for that extent continental Equal-area preserve area, accept distance error

Gotchas and Edge Cases

  • Feeds crossing the antimeridian. The bounding box wraps, the centroid lands on the wrong side of the world, and the chosen zone is nonsense. Rare enough to handle by exception rather than by general code, but worth an explicit check on any feed whose longitude span exceeds 180 degrees.
  • Coordinates of exactly 0, 0. A missing coordinate written as zero rather than left blank puts a stop in the Gulf of Guinea and drags the feed extent across the world. Filter implausible coordinates before computing the extent, not after.
  • shape_dist_traveled after reprojection. It is published in whatever unit the agency chose and is unaffected by reprojection. Do not recompute it from projected coordinates unless you intend to replace the agency’s values entirely — see route geometry extraction.
  • Mixing projected and unprojected frames in a join. A spatial join across mismatched CRS returns an empty result rather than an error. Carrying the target EPSG code alongside every projected frame, as choose_and_project returns it, makes the mismatch checkable.
  • pyproj versions before 2.0. The Transformer API and always_xy do not exist there, and the older pyproj.transform function has the opposite axis convention. Pin the version.

Frequently Asked Questions

Why not just use Web Mercator?

Because it distorts distance by the secant of the latitude — about 1.6x at 50 degrees north. Stop spacing, buffer radii and snapping distances computed in it are all wrong by that factor. Web Mercator is a display projection and nothing else.

What does always_xy actually do?

It forces the transformer to take and return coordinates in longitude, latitude order rather than the authority-defined order, which for EPSG:4326 is latitude first. Omitting it is the most common pyproj bug and it silently swaps every coordinate.

Should I transform point by point?

No. pyproj transforms whole arrays in one call, which is roughly two orders of magnitude faster than looping. On a shapes.txt with 214,000 vertices that is the difference between milliseconds and half a minute.

What if the feed spans more than one UTM zone?

Pick the zone containing the centroid and accept the distortion at the edges, or move to a national grid or an equal-area projection covering the whole extent. A feed spanning more than about six degrees of longitude has outgrown UTM.