1. Home
  2. GTFS Feed Architecture Fundamentals
  3. Spatial Analysis And Route Geometry
  4. Geopandas Spatial Joins For Transit Stops

GeoPandas Spatial Joins for Transit Stops

Load stops.txt into a GeoDataFrame in EPSG:4326, project both the stops and your polygon layer to the same metric CRS, then call gpd.sjoin with predicate="within" for point-in-polygon matches or predicate="intersects" when you have buffered stops into walk catchments. The join attaches polygon attributes to each stop; group the result by the zone key to count stops and sum service per area. The projection step is mandatory — buffers and boundaries are only metric in a projected CRS.

Root Cause Analysis

A spatial join answers a question no shared identifier can: which zone does this stop fall in? GTFS stops carry no zone_id for census tracts, fare districts, or council boundaries — those live in separate polygon layers with no common key — so the relationship must be computed from geometry. gpd.sjoin does exactly that, pairing each stop with the polygons it satisfies a spatial predicate against.

Two things routinely break the join. The first is coordinate mismatch. sjoin compares raw coordinate numbers through a spatial index; it does not reconcile reference systems. If stops sit in EPSG:4326 (degrees) and the polygon layer sits in a projected grid (metres), the two live in different numeric universes and nothing matches — recent geopandas raises rather than returning silent emptiness, but the fix is the same: reproject both to one coordinate reference system. And it must be a metric one, because the moment you buffer a stop into a walk catchment the buffer radius is in the CRS’s units — a buffer(400) in EPSG:4326 means 400 degrees, which is nonsense. Choosing that CRS is covered in the projected-CRS selection guide.

The second is cardinality. sjoin emits one row per matching pair, not one row per stop. A point on a shared boundary, overlapping polygon layers, or a buffered catchment straddling several zones all produce multiple rows for one stop_id. Treating the raw join output as one-row-per-stop silently inflates every count you compute from it.

Production-Ready Python Implementation

The script loads stops, projects to an auto-derived UTM zone, joins against a polygon layer with a configurable predicate, resolves multi-matches, folds in trip counts from stop_times.txt, and aggregates by zone.

python 
import math
import logging
import pandas as pd
import geopandas as gpd

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 join_stops_to_zones(
    feed_dir: str,
    zones_path: str,
    zone_key: str = "zone_name",
    catchment_m: float | None = None,
) -> tuple[gpd.GeoDataFrame, pd.DataFrame]:
    """Attach polygon attributes to GTFS stops and aggregate by zone."""
    # --- Load stops with explicit dtypes ---
    stops = pd.read_csv(
        f"{feed_dir}/stops.txt",
        dtype={"stop_id": str, "stop_name": str, "stop_lat": float, "stop_lon": float},
        usecols=["stop_id", "stop_name", "stop_lat", "stop_lon"],
    )
    stops_gdf = gpd.GeoDataFrame(
        stops,
        geometry=gpd.points_from_xy(stops["stop_lon"], stops["stop_lat"]),
        crs="EPSG:4326",
    )

    # --- Derive a metric CRS and project BOTH layers to it ---
    target_epsg = utm_epsg(float(stops["stop_lon"].mean()), float(stops["stop_lat"].mean()))
    stops_gdf = stops_gdf.to_crs(target_epsg)
    zones = gpd.read_file(zones_path).to_crs(target_epsg)
    if zone_key not in zones.columns:
        raise KeyError(f"Zone layer has no column '{zone_key}'; got {list(zones.columns)}")
    zones = zones[[zone_key, "geometry"]]
    logging.info("Both layers in EPSG:%d | %d stops, %d zones",
                 target_epsg, len(stops_gdf), len(zones))

    # --- Optional: buffer stops into walk catchments (metres) ---
    if catchment_m:
        left = stops_gdf.copy()
        left["geometry"] = left.geometry.buffer(catchment_m)
        predicate = "intersects"   # areas: catch every zone the buffer touches
    else:
        left = stops_gdf
        predicate = "within"       # points: which zone contains the stop

    # --- The spatial join ---
    joined = gpd.sjoin(left, zones, how="left", predicate=predicate)
    logging.info("sjoin('%s') produced %d rows for %d stops",
                 predicate, len(joined), stops_gdf["stop_id"].nunique())

    # --- Resolve cardinality: keep first zone per stop for a clean 1:1 ---
    unmatched = joined[joined[zone_key].isna()]["stop_id"].nunique()
    if unmatched:
        logging.warning("%d stops matched no zone", unmatched)
    stops_zoned = (
        joined.sort_values(["stop_id", zone_key])
        .drop_duplicates(subset="stop_id", keep="first")
        .drop(columns=["index_right"], errors="ignore")
    )

    # --- Fold in trips-per-stop from stop_times.txt ---
    stop_times = pd.read_csv(
        f"{feed_dir}/stop_times.txt",
        dtype={"trip_id": str, "stop_id": str, "stop_sequence": "Int64"},
        usecols=["trip_id", "stop_id"],
    )
    trips_per_stop = (
        stop_times.groupby("stop_id", sort=False)["trip_id"].nunique()
        .rename("trip_count").reset_index()
    )
    stops_zoned = stops_zoned.merge(trips_per_stop, on="stop_id", how="left")
    stops_zoned["trip_count"] = stops_zoned["trip_count"].fillna(0).astype(int)

    # --- Aggregate by zone ---
    by_zone = (
        stops_zoned.groupby(zone_key, dropna=False)
        .agg(stop_count=("stop_id", "nunique"), total_trips=("trip_count", "sum"))
        .reset_index()
        .sort_values("stop_count", ascending=False)
    )
    return stops_zoned, by_zone


# --- Entry point ---
# stops_zoned, by_zone = join_stops_to_zones(
#     "/data/gtfs/king_county", "zones/council_districts.gpkg",
#     zone_key="district", catchment_m=None,
# )
# print(by_zone.to_string(index=False))

Step-by-Step Walkthrough

Both layers, one CRS. stops_gdf.to_crs(target_epsg) and gpd.read_file(zones_path).to_crs(target_epsg) land both inputs in the same metric grid. This is the line that makes the join valid; without it, geopandas either raises on the CRS mismatch or, on older versions, returns an empty match set.

Buffer decides the predicate. When catchment_m is set, each stop becomes a circular polygon of that radius in metres and the predicate switches to intersects, so the join captures every zone the walk catchment overlaps. With no buffer, stops stay as points and within answers the strict containment question. The two paths are deliberately mutually exclusive — a buffered point is no longer “within” anything in the useful sense.

left join keeps unmatched stops. how="left" preserves stops that fall outside every polygon (a stop just past a city boundary, say) with a null zone rather than dropping them. The count of nulls is logged, because silently losing stops is worse than surfacing that they are outside the layer.

Resolving multiple matches. sjoin can emit several rows per stop. Sorting by (stop_id, zone_key) and drop_duplicates(keep="first") collapses to one deterministic zone per stop. For overlapping layers where the largest-overlap zone is the right answer, replace this with an area-weighted resolution instead of first-match.

Service-weighting the stops. Counting stops per zone treats a once-a-day flag stop the same as a busy interchange. Merging trip_count from stop_times.txt and summing it per zone gives a service-weighted view — the total_trips column — which is usually the metric planners actually want.

Verification and Output

Guard the join with cardinality and coverage assertions before reporting:

python 
stops_zoned, by_zone = join_stops_to_zones(
    "/data/gtfs/king_county", "zones/council_districts.gpkg", zone_key="district",
)

# The resolved frame must be exactly one row per stop
assert stops_zoned["stop_id"].is_unique, "sjoin cardinality not resolved to 1:1"

# Stop counts across zones (plus unmatched) must reconcile to the feed total
feed_stops = pd.read_csv("/data/gtfs/king_county/stops.txt", dtype={"stop_id": str})
assert by_zone["stop_count"].sum() == feed_stops["stop_id"].nunique(), (
    "Zone stop counts do not sum to the number of stops in the feed"
)

# No negative or NaN service totals
assert (by_zone["total_trips"] >= 0).all(), "Negative trip total after aggregation"

print(by_zone.head(10).to_string(index=False))

The output is two frames: stops_zoned, one row per stop with its zone and trip count, and by_zone, one row per zone with stop_count and total_trips. Because the join uses a left join, the zone totals plus any unmatched (null-zone) stops always reconcile to the feed’s stop count — that reconciliation assertion is the single best guard against silently dropping stops in a CRS mismatch.

Gotchas and Edge Cases

  • Buffer radius in the wrong unit. buffer(400) means 400 units of the current CRS. In EPSG:4326 that is 400 degrees — a catchment larger than a hemisphere. Always project to metres before buffering, and sanity-check that a stop’s buffered area is on the order of half a square kilometre, not the whole globe.

  • Boundary-straddling points under intersects. A stop whose coordinate lands exactly on a shared polygon edge matches both adjacent polygons with intersects. Use within for strict containment, or snap-and-resolve to the polygon of largest overlap, so a single stop does not double-count across neighbouring zones.

  • index_right collisions on chained joins. gpd.sjoin writes the right layer’s index into an index_right column. If you run a second join on the result without dropping it, geopandas raises a column-collision error. Drop index_right after each join, as the script does.

  • Reused stop_id in merged multi-agency feeds. When stops from several agencies share an identifier space, deduplicating the join on stop_id alone silently merges distinct physical stops. Prefix identifiers per agency first; the same precaution applies to the stops.txt and stop_times.txt join keys throughout the pipeline.

Frequently Asked Questions

Why must both layers share the same CRS before a geopandas spatial join?

gpd.sjoin compares raw coordinate values, so if stops are in EPSG:4326 and polygons are in a projected CRS the geometries occupy different numeric spaces and nothing matches. geopandas raises on mismatched CRS in recent versions; always reproject both layers to one metric CRS first.

When should I use predicate='within' versus 'intersects'?

Use within for a point-in-polygon question — which zone contains this stop. Use intersects when the left geometry has area, such as a buffered walk catchment, and you want every polygon it touches. A point exactly on a shared boundary can match two polygons under intersects.

Why does my spatial join return more rows than I have stops?

gpd.sjoin emits one row per matching pair, so a stop that falls in overlapping polygons or a buffered stop that touches several zones produces multiple rows. Deduplicate on stop_id, or resolve to a single best match, before treating the result as one row per stop.