Measuring Stop Spacing and Route Length
Project once into a metric CRS, build one LineString per shape, snap each stop to a position along that line, and take consecutive differences. Measuring straight-line distances between stops instead cuts every corner and understates a winding route by ten to twenty per cent. Then report the spacing as a distribution — percentiles, not a mean — because almost every real route is dense in the centre and sparse at the edges, and the average falls in the gap between the two. The geometry this is built on is covered in spatial analysis and route geometry.
Root Cause Analysis
Stop spacing and route length are the two most requested numbers in transit analysis, and both are routinely computed in a way that makes them wrong.
Measuring in degrees. The most common error, and the one that varies by city. Euclidean distance on latitude and longitude understates east–west distance by the cosine of the latitude — 26% at 40 degrees, 36% at 50, half at 60. The number looks plausible everywhere and is only correct near the equator. Projecting first is the fix, and it has to happen before any arithmetic.
Measuring between stops rather than along the route. A straight line from one stop to the next ignores the road. On a grid network the error is small; on a route that follows a river or loops through an industrial estate it is large, and it always understates, so route lengths computed this way are systematically short.
Averaging a bimodal distribution. A typical urban route has stops every 200 metres in the centre and every 800 metres on the outer sections. The mean lands around 400 metres, which is a spacing the route does not use anywhere. Any planning conclusion drawn from it is drawn from a number describing no part of the network.
Treating a route as one line. A trunk route may carry a dozen shape_id values — two directions, short turns, peak variants, diversions. “The length of route 47” is not well defined until a pattern is chosen, and different tools choosing differently is why two reports on the same feed disagree.
Production-Ready Python Implementation
"""Measure GTFS route length and stop spacing along the route shape."""
from __future__ import annotations
import logging
from dataclasses import dataclass, field
import numpy as np
import pandas as pd
from shapely.geometry import LineString, Point
log = logging.getLogger("gtfs.geometry.measure")
# A stop further than this from its route's shape is not on that pattern.
MAX_SNAP_OFFSET_M = 120.0
@dataclass
class RouteMeasurement:
route_id: str
shape_id: str
variants: int # how many distinct shapes this route has
length_m: float
stop_count: int
spacings_m: list[float] = field(default_factory=list)
off_shape_stops: int = 0
@property
def percentiles(self) -> dict[str, float]:
if not self.spacings_m:
return {}
a = np.asarray(self.spacings_m)
return {"p10": float(np.percentile(a, 10)), "p50": float(np.percentile(a, 50)),
"p90": float(np.percentile(a, 90)), "min": float(a.min()),
"max": float(a.max())}
def describe(self) -> str:
p = self.percentiles
if not p:
return f"route {self.route_id}: no usable spacing"
return (f"route {self.route_id} ({self.length_m / 1000:.1f} km, "
f"{self.stop_count} stops, {self.variants} pattern(s)): spacing "
f"p10 {p['p10']:.0f} m, median {p['p50']:.0f} m, p90 {p['p90']:.0f} m")
def build_lines(projected_shapes: pd.DataFrame) -> dict[str, LineString]:
"""One LineString per shape_id, vertices in sequence order."""
lines: dict[str, LineString] = {}
ordered = projected_shapes.sort_values(["shape_id", "shape_pt_sequence"])
for shape_id, group in ordered.groupby("shape_id"):
coords = list(zip(group["x"], group["y"]))
if len(coords) < 2:
log.warning("shape %s has %d vertex/vertices — skipped", shape_id, len(coords))
continue
lines[str(shape_id)] = LineString(coords)
return lines
def representative_shape(trips: pd.DataFrame,
lines: dict[str, LineString]) -> dict[str, tuple[str, int]]:
"""route_id -> (longest shape_id, how many distinct shapes that route has)."""
chosen: dict[str, tuple[str, int]] = {}
for route_id, group in trips.dropna(subset=["shape_id"]).groupby("route_id"):
candidates = [s for s in group["shape_id"].astype(str).unique() if s in lines]
if not candidates:
continue
longest = max(candidates, key=lambda s: lines[s].length)
chosen[str(route_id)] = (longest, len(candidates))
return chosen
def measure_route(route_id: str, shape_id: str, variants: int, line: LineString,
stop_points: list[tuple[str, Point]]) -> RouteMeasurement:
"""stop_points: the route's stops in stop_sequence order, already projected."""
positions: list[float] = []
off_shape = 0
for _, point in stop_points:
offset = line.distance(point)
if offset > MAX_SNAP_OFFSET_M:
off_shape += 1
continue
positions.append(line.project(point)) # distance ALONG the line
# Sort defensively: a loop route can snap out of order, and negative
# spacings would corrupt every percentile below.
positions.sort()
spacings = [b - a for a, b in zip(positions, positions[1:]) if b > a]
measurement = RouteMeasurement(
route_id=route_id, shape_id=shape_id, variants=variants,
length_m=float(line.length), stop_count=len(positions),
spacings_m=spacings, off_shape_stops=off_shape)
if off_shape:
log.warning("route %s: %d stop(s) more than %.0f m from shape %s — they belong "
"to a different pattern", route_id, off_shape, MAX_SNAP_OFFSET_M,
shape_id)
return measurement
def measure_feed(projected_stops: pd.DataFrame, projected_shapes: pd.DataFrame,
trips: pd.DataFrame, stop_times: pd.DataFrame
) -> list[RouteMeasurement]:
lines = build_lines(projected_shapes)
chosen = representative_shape(trips, lines)
xy = {str(r.stop_id): Point(r.x, r.y)
for r in projected_stops.itertuples(index=False)}
# One representative trip per route supplies the stop order.
trip_of_route = {}
for route_id, (shape_id, _) in chosen.items():
candidates = trips[(trips["route_id"] == route_id)
& (trips["shape_id"].astype(str) == shape_id)]
if not candidates.empty:
trip_of_route[route_id] = str(candidates.iloc[0]["trip_id"])
calls = (stop_times.sort_values(["trip_id", "stop_sequence"])
.groupby("trip_id")["stop_id"].apply(list).to_dict())
out: list[RouteMeasurement] = []
for route_id, (shape_id, variants) in chosen.items():
trip_id = trip_of_route.get(route_id)
if trip_id is None or trip_id not in calls:
continue
points = [(s, xy[str(s)]) for s in calls[trip_id] if str(s) in xy]
if len(points) < 2:
continue
out.append(measure_route(route_id, shape_id, variants, lines[shape_id], points))
log.info("measured %d route(s) of %d with usable geometry", len(out), len(chosen))
return out
Step-by-Step Walkthrough
Everything arrives already projected. measure_feed takes projected_stops and projected_shapes with x and y columns rather than doing the projection itself. That keeps the CRS decision in one place and makes it impossible for this module to accidentally measure degrees — there is no latitude column to reach for.
line.project(point) is the whole idea. Shapely’s project returns the distance along the line to the nearest point on it, which is precisely “how far into the route is this stop”. Consecutive differences of those positions give spacing along the road rather than across it, and no separate distance calculation is needed.
line.distance(point) gates it. A stop more than 120 metres from the line is not on this pattern — it belongs to a branch or a variant that the representative shape does not cover. Snapping it anyway would place it at whichever end of the line happened to be nearest and produce a wild spacing. Counting those separately keeps the measurement honest about what it covered.
Positions are sorted before differencing. On a loop route, project can return positions out of stop order where the line passes near itself. Sorting means every spacing is positive; the if b > a guard then drops the zero-length pairs that a repeated stop produces.
The representative shape is the longest. Short turns and branches are subsets of the full pattern, so the longest shape is the one that covers the most of the route. variants is carried alongside so the reader knows whether the number describes a simple route or one of twelve patterns.
Percentiles are computed on demand, not stored. RouteMeasurement keeps the raw spacings, which means a caller can compute any statistic — a different percentile, a histogram, a comparison against a planning standard — without re-running the measurement.
Verification and Output
def verify(measurement: RouteMeasurement, line: LineString) -> None:
assert measurement.length_m > 0, "zero-length route geometry"
assert all(s > 0 for s in measurement.spacings_m), "non-positive spacing"
total = sum(measurement.spacings_m)
assert total <= measurement.length_m + 1.0, (
f"stops span {total:.0f} m along a {measurement.length_m:.0f} m line — "
"positions were not taken along the shape")
if measurement.stop_count >= 2:
assert len(measurement.spacings_m) <= measurement.stop_count - 1
straight = LineString([line.coords[0], line.coords[-1]]).length
assert measurement.length_m >= straight - 1.0, (
"the route is shorter than the straight line between its ends, which is "
"geometrically impossible")
The last assertion is a cheap sanity check on the whole projection chain: a path between two points cannot be shorter than the straight line between them, and if it is, the coordinates are not in the space you think they are.
Output across a mixed network:
INFO gtfs.geometry.measure: measured 204 route(s) of 210 with usable geometry
route 1 (11.4 km, 42 stops, 6 pattern(s)): spacing p10 168 m, median 246 m, p90 611 m
route 77 (18.9 km, 31 stops, 4 pattern(s)): spacing p10 284 m, median 512 m, p90 1840 m
route RL (34.2 km, 17 stops, 2 pattern(s)): spacing p10 940 m, median 1720 m, p90 4100 m
Three lines that describe three genuinely different things — an urban trunk bus, a suburban route, a rail line — and the p10-to-p90 range on route 77 shows the bimodality directly. Its mean spacing would be around 640 metres, a value that occurs nowhere on the route.
Gotchas and Edge Cases
- Routes with no
shape_id. They are skipped entirely, and the count in the log says how many. Falling back to straight lines between stops would produce a number in the same column as the measured ones, which invites comparison between two quantities that are not the same. shape_dist_traveleddisagreeing with the measurement. The agency’s values are optional, unit-unspecified and occasionally reset mid-shape. Compare them against the measured length; where they agree within a few per cent, they are safe to use and much cheaper than measuring.- Circular routes. The first and last stop are the same place, so the final spacing is either zero or the whole loop. Neither is meaningful; detect a closed line and drop the wrap-around pair.
- Stops served twice on one trip. Their two positions along the line are different, which is correct, and sorting preserves both. This is the one case where sorting could reorder relative to
stop_sequence— and along-the-line order is the right order for spacing. - Very dense shapes. A shape with a vertex every metre makes
projectslower without making it more accurate. Simplifying the line to a one-metre tolerance before measuring is safe and can halve the runtime on rail geometries.
Frequently Asked Questions
Should route length be measured along the shape or between stops?
Along the shape. Straight lines between consecutive stops cut every corner, and on a winding route that understates the length by ten to twenty per cent. The shape is the path the vehicle actually takes, which is what a length is supposed to describe.
Which shape represents a route with a dozen variants?
Usually the longest, because short turns and branches are subsets of the full pattern. Report how many variants existed alongside the measurement, so a reader knows the figure describes one pattern rather than the whole route.
Why report percentiles instead of an average stop spacing?
Because the distribution is bimodal on almost every real route: dense stops in the centre, sparse ones on the outer sections. The mean falls in the gap between the two clusters and describes nothing that exists.
Can I trust shape_dist_traveled instead of measuring?
Only after checking it. It is optional, its unit is unspecified, and a minority of feeds reset it mid-shape. Measure it yourself and compare; where the two agree, the agency’s values are fine to use.
Related
- Extracting Route Geometry from GTFS Shapes — building the lines this measures
- Snapping GTFS Stops to Shapes in Python — the projection step, in detail
- Up: Spatial Analysis and Route Geometry — the geometry model behind both
- Section: GTFS Feed Architecture & Fundamentals · Home