Extracting Route Geometry from GTFS Shapes

Sort shapes.txt by (shape_id, shape_pt_sequence), assemble the ordered points into one LineString per shape_id, then join through trips.txt to map each shape to its route_id. Because a route owns many shapes, project to a metric CRS, measure each shape’s length, and keep the longest shape per route as its single representative geometry. The result is one clean LineString per route, ready for mapping, length reporting, or as the input to stop snapping.

Root Cause Analysis

The reason “extract route geometry” is a task and not a column lookup is that GTFS never stores route geometry. It stores shape points, and it stores them one level of indirection away from routes.

shapes.txt holds the vertices — rows of (shape_pt_lat, shape_pt_lon) grouped by shape_id and ordered by shape_pt_sequence. But shapes.txt has no route_id. The link runs route_id → trips.txt → shape_id, and it is many-to-many in practice: a single route runs inbound and outbound trips, peak express and all-day local patterns, and occasional short-turns, and each of those can reference a distinct shape_id. So “the geometry of route 47” is really a set of candidate shapes, and extracting one geometry means choosing among them.

Two failure points sit inside that assembly. First, ordering: shape_pt_sequence is only guaranteed to increase, not to be contiguous or pre-sorted, and if it is read as a string it sorts lexically so that "10" lands before "2". An unsorted assembly produces a LineString that jumps back and forth, self-intersects, and reports a meaningless length. Second, measurement: comparing shape lengths to pick the longest requires metric units, which means projecting out of WGS84 first, exactly as the coordinate reference systems guidance requires.

Production-Ready Python Implementation

This script reads shapes.txt and trips.txt, assembles a LineString per shape_id, maps shapes to routes, projects to an auto-derived UTM zone, and selects the longest shape per route as the representative geometry — returning a GeoDataFrame keyed by route_id.

import math
import logging
import numpy as np
import pandas as pd
import geopandas as gpd
from shapely import linestrings  # shapely 2.x vectorized constructor

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 extract_route_geometry(feed_dir: str) -> gpd.GeoDataFrame:
    """Return one representative LineString per route_id (longest shape)."""
    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,   # force integer so ordering is numeric
        },
        usecols=["shape_id", "shape_pt_lat", "shape_pt_lon", "shape_pt_sequence"],
    )
    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"],
    )

    # --- 1. Order rows: shape_id groups, numeric sequence within each ---
    shapes = shapes.sort_values(["shape_id", "shape_pt_sequence"], kind="stable")

    # --- 2. Assemble one LineString per shape_id (vectorized) ---
    codes, shape_ids = pd.factorize(shapes["shape_id"], sort=False)
    coords = np.column_stack([shapes["shape_pt_lon"].to_numpy(),
                              shapes["shape_pt_lat"].to_numpy()])
    geoms = linestrings(coords, indices=codes)
    shape_lines = gpd.GeoDataFrame(
        {"shape_id": shape_ids}, geometry=geoms, crs="EPSG:4326"
    )

    # Drop degenerate shapes (a single vertex cannot form a line)
    valid = shape_lines.geometry.apply(lambda ls: len(ls.coords) >= 2)
    if (~valid).any():
        logging.warning("Dropping %d single-point shapes", int((~valid).sum()))
    shape_lines = shape_lines[valid].reset_index(drop=True)

    # --- 3. Map shape_id -> route_id via trips.txt ---
    shape_to_route = (
        trips.dropna(subset=["shape_id"])[["shape_id", "route_id"]].drop_duplicates()
    )
    shape_lines = shape_lines.merge(shape_to_route, on="shape_id", how="inner")

    # --- 4. Project to a metric CRS so length is in metres ---
    centroid = shape_lines.geometry.union_all().centroid
    target_epsg = utm_epsg(centroid.x, centroid.y)
    logging.info("Projecting %d shapes to EPSG:%d", len(shape_lines), target_epsg)
    shape_lines = shape_lines.to_crs(target_epsg)
    shape_lines["length_m"] = shape_lines.geometry.length

    # --- 5. Longest shape per route is the representative geometry ---
    idx = shape_lines.groupby("route_id", sort=False)["length_m"].idxmax()
    representative = (
        shape_lines.loc[idx, ["route_id", "shape_id", "length_m", "geometry"]]
        .reset_index(drop=True)
    )
    logging.info(
        "Selected representative geometry for %d routes (from %d shapes)",
        representative["route_id"].nunique(), len(shape_lines),
    )
    return gpd.GeoDataFrame(representative, geometry="geometry", crs=target_epsg)


# --- Entry point ---
# routes_geom = extract_route_geometry("/data/gtfs/trimet")
# routes_geom.to_parquet("cache/route_geometry.parquet")

Step-by-Step Walkthrough

Forcing a numeric sequence. shape_pt_sequence is declared int in the dtype map so the sort is numeric, not lexical. This one declaration prevents the most common route-geometry bug: a shape assembled in the order 1, 10, 11, 2, 20, 3 …, which produces a line that criss-crosses itself.

Sort before assembly. sort_values(["shape_id", "shape_pt_sequence"], kind="stable") groups the rows by shape and orders vertices within each group. kind="stable" preserves input order among equal keys, which matters if a feed has duplicate sequence numbers you have not yet cleaned.

Vectorized LineString construction. pd.factorize turns shape_id into integer group codes, and shapely 2.x’s linestrings(coords, indices=codes) builds every line in one C-level call — far faster than a Python groupby().apply loop on a national feed with tens of millions of shape points.

Mapping shapes to routes. trips.txt is reduced to distinct (shape_id, route_id) pairs and inner-joined onto the assembled lines. The inner join drops any shape that no trip references — orphaned geometry that would otherwise inflate your route set.

Projecting to measure length. The union centroid picks a UTM zone; to_crs reprojects; geometry.length is then in metres. Only after projection is comparing lengths across shapes valid — the whole longest-shape selection depends on it.

Longest shape per route. groupby("route_id")["length_m"].idxmax() returns, for each route, the index of its longest shape. Selecting those rows yields exactly one representative LineString per route. The longest shape is chosen deliberately: it captures the full end-to-end pattern rather than a short-turn or a peak-only express variant.

Verification and Output

Confirm the geometry is well-formed and the representative set is complete:

routes_geom = extract_route_geometry("/data/gtfs/trimet")

# Exactly one representative geometry per route_id
assert routes_geom["route_id"].is_unique, "Duplicate route_id in representative set"

# No degenerate or invalid geometry survived
assert routes_geom.geometry.is_valid.all(), "Invalid geometry in output"
assert (routes_geom.geometry.length > 0).all(), "Zero-length route geometry"

# Lengths are metric and plausible (metres, not degrees or millimetres)
lengths_km = routes_geom["length_m"] / 1000
assert lengths_km.between(0.1, 300).all(), (
    "Route length outside 0.1-300 km — check projection or shape sorting"
)

print(routes_geom.assign(length_km=lengths_km.round(2))
      [["route_id", "shape_id", "length_km"]]
      .sort_values("length_km", ascending=False)
      .head(10).to_string(index=False))

The output is a GeoDataFrame with one row per route_id, its chosen shape_id, a metric length_m, and the projected LineString. A well-formed feed shows route lengths clustered in single-to-low-double-digit kilometres for bus and larger for rail; a length in the hundreds of thousands means the geometry is still in degrees and the projection step did not run.

Gotchas and Edge Cases

• Lexical sequence sorting. If shape_pt_sequence is read as a string, "100" sorts before "20" and the LineString self-crosses. Always pin it to an integer dtype, and if a feed uses non-contiguous sequences (10, 20, 30 for later insertions), that is spec-legal — only the ordering matters, not the gaps.

• Direction collapses into one line. Picking the longest shape per route_id merges inbound and outbound into a single representative. If your analysis is direction-sensitive, group by (route_id, direction_id) from trips.txt instead, keeping one representative per direction.

• Routes with no shape at all. Some feeds ship trips.txt rows with an empty shape_id, or omit shapes.txt entirely. Those routes produce no geometry here; handle them by falling back to a straight-line path through each trip’s ordered stop coordinates, joined via stops.txt and stop_times.txt.

• Duplicate points and zero-length segments. Consecutive identical vertices are common and harmless for assembly, but they can trip up downstream simplification. If you plan to simplify geometry, drop consecutive duplicates before measuring so segment counts and lengths stay honest.

Frequently Asked Questions

Related

• Snapping GTFS Stops to Shapes in Python — feed these LineStrings into linear referencing to place stops along the route
• GeoPandas Spatial Joins for Transit Stops — join the projected geometry against zones and service areas
• Up: Spatial Analysis and Route Geometry — how geometry assembly opens the full spatial pipeline
• Section: GTFS Feed Architecture & Fundamentals — the schema rules governing shapes.txt and trips.txt · Home