Performing Point-in-Polygon Joins for Store Catchments

Assigning transaction points, prospect centroids, and competitor locations to predefined trade area boundaries is the geometric operation that converts raw coordinate streams into catchment-level intelligence inside the Demographic Data Integration & Spatial Joins pipeline.

This page covers the spatial join that sits at the heart of retail site selection automation: matching each point to the catchment polygon that contains it, then aggregating transactions, revenue, and demographics per catchment. The output becomes the primary key for every downstream model — revenue forecasting, lease evaluation, and network optimization all depend on points being attributed deterministically and without silent loss. Getting it right requires disciplined coordinate reference system (CRS) handling, topology validation, and explicit predicate selection.

Concept and theory: containment as a topological predicate #

A point-in-polygon join is not a relational key match — it resolves a geometric relationship between two layers. For each point geometry, the engine asks whether the point lies within the area enclosed by a polygon’s exterior ring (and outside any interior rings, or holes). Classical algorithms answer this with a ray-casting test (counting boundary crossings of a half-line extended from the point) or the winding-number method. Modern geospatial stacks delegate this to the GEOS library, which both GeoPandas and PostGIS wrap, so the predicate semantics are identical across the Python and database paths.

Two predicate distinctions matter for catchments:

• within vs contains: these are inverse relations. A point is within a polygon; a polygon contains a point. Choosing the wrong direction in a join silently returns zero matches.
• Boundary handling (within vs intersects): within excludes points that fall exactly on a shared edge, while intersects includes them. With adjacent, non-overlapping catchments — the normal case for a tiled trade area surface — a point on a shared border is a tie that must be broken deterministically.

The cost of a naive join is quadratic: testing every point against every polygon is O(n·m). Production pipelines avoid this with a bounding-box spatial index (an R-tree). The index reduces the containment test to a small candidate set per point — typically one or two polygons — before the exact GEOS predicate runs. This is the same indexing discipline described in the parent spatial joins overview, and it is non-negotiable once point counts pass a few thousand.

A subtle prerequisite underlies all of it: both layers must share an identical, projected CRS. Containment is unambiguous, but the spatial index is built on planar bounding boxes, and a mismatched or geographic CRS produces incorrect candidate sets and angular distance distortion. Establish CRS alignment before any geometric operation runs.

Architecture overview #

The catchment join is a five-stage flow: align both layers to a projected CRS and repair geometry, build a spatial index on the polygons, execute the within join, fall back to a nearest-neighbor match for unmatched edge points, and aggregate metrics per catchment.

Stages are deterministic and idempotent: re-running the flow on the same inputs produces byte-identical metrics, which is what makes the join safe to schedule on a fixed cadence.

Configuration parameters #

The table below lists the settings that govern join correctness and the retail-specific defaults that production pipelines should adopt. Treat any deviation from these defaults as a decision that needs documenting in the run config.

The combination of a projected target_crs and predicate="within" is the safe baseline. Overlapping catchments — common across multi-format portfolios where a supercenter and an express format share territory — are the main reason to relax it, and they always demand an explicit tie-breaker so a single transaction is never counted twice.

Step-by-step Python implementation #

The implementation loads catchments, repairs their topology, builds point geometries from raw longitude/latitude, enforces a shared projected CRS, runs the within join, applies a nearest-neighbor fallback for stray points, and aggregates per catchment. Every coordinate operation asserts a CRS before it runs, in line with the data foundations standard for store coordinate validation.

import geopandas as gpd
import pandas as pd
from shapely.validation import make_valid
import logging

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


def execute_catchment_join(
    points_path: str,
    catchments_path: str,
    target_crs: str = "EPSG:5070",
) -> pd.DataFrame:
    """
    Joins transaction points to retail catchment polygons.
    Applies a nearest-neighbor fallback for unmatched edge cases.
    """
    # 1. Load and repair catchment topologies (removes self-intersections / slivers)
    catchments = gpd.read_file(catchments_path)
    catchments["geometry"] = catchments["geometry"].apply(make_valid)

    # 2. Ingest coordinates and construct point geometries
    points_df = pd.read_csv(points_path)
    points_df["geometry"] = gpd.points_from_xy(
        points_df["longitude"], points_df["latitude"]
    )
    points_gdf = gpd.GeoDataFrame(points_df, geometry="geometry", crs="EPSG:4326")

    # 3. Enforce an identical projected CRS for indexing and the predicate
    if catchments.crs is None:
        raise ValueError("Catchment dataset lacks a CRS definition. Assign before join.")
    catchments = catchments.to_crs(target_crs)
    points_gdf = points_gdf.to_crs(target_crs)

    # 4. Execute the spatial join with an explicit predicate.
    #    how="left" keeps every point so unmatched rows stay visible.
    joined = gpd.sjoin(points_gdf, catchments, how="left", predicate="within")

    # 5. Nearest-neighbor fallback for boundary / GPS-drift edge cases
    unmatched_mask = joined["index_right"].isna()
    if unmatched_mask.any():
        logging.warning(
            "%d points fell outside all catchments. Applying nearest-neighbor fallback.",
            unmatched_mask.sum(),
        )
        # sjoin_nearest runs on the original unmatched points (no index_right yet)
        unmatched_points = points_gdf[unmatched_mask.values]
        nearest = gpd.sjoin_nearest(
            unmatched_points, catchments, how="left", max_distance=5000
        )
        joined.loc[unmatched_mask, "index_right"] = nearest["index_right"].values

    # 6. Aggregate transaction metrics per catchment
    catchment_metrics = (
        joined.groupby("index_right")
        .agg(
            transaction_count=("transaction_id", "count"),
            total_revenue=("revenue_usd", "sum"),
            avg_ticket=("revenue_usd", "mean"),
        )
        .reset_index()
        .rename(columns={"index_right": "catchment_id"})
    )

    return catchment_metrics

For teams that keep catchments in the database rather than in files, the equivalent PostGIS form is a single indexed query — the geometry column carries a GiST index and ST_Within runs the same GEOS predicate:

SELECT c.catchment_id,
       count(*)        AS transaction_count,
       sum(p.revenue_usd) AS total_revenue
FROM   transactions p
JOIN   catchments   c
  ON   ST_Within(p.geom, c.geom)   -- both geoms in EPSG:5070, GiST-indexed
GROUP  BY c.catchment_id;

Keep the Python and SQL paths interchangeable: both must use the same projected CRS and the same predicate so their outputs reconcile during validation.

Edge cases and failure modes #

The join rarely fails loudly. The dangerous failures are silent — points vanish or land in the wrong catchment without raising an exception.

• CRS mismatch: joining a geographic layer against a projected one (or two different projected CRSs) yields empty or wrong matches with no error. The CRS assertion in step 3 turns this into a hard failure.
• Inverted coordinates: a longitude/latitude swap places points in the wrong hemisphere; coverage collapses to near zero. Bounding-box sanity checks catch it immediately.
• Sliver and self-intersecting polygons: hand-digitized or clipped catchments produce micro-artifacts that trigger GEOS TopologyException or false negatives. The make_valid() repair handles most, but persistent cases need targeted remediation — see fixing sliver polygons in spatial join operations.
• GPS drift on boundaries: POS, mobile, and third-party feeds carry sub-meter to tens-of-meters drift, pushing legitimate customers just outside a catchment. The bounded sjoin_nearest fallback recovers these without dragging far-flung outliers in.
• Overlapping catchments: predicate="within" returns multiple rows per point, inflating counts. Use intersects with a deterministic tie-breaker (shortest distance to the anchor store centroid, or a documented priority rank).

Performance and scaling #

The single biggest lever is the spatial index, which GeoPandas builds automatically for sjoin. Beyond that, throughput on national portfolios depends on memory discipline:

• Batch by region: partition points by state or metro FIPS and join against the matching catchment subset. This keeps the candidate R-tree small and lets batches run in parallel.
• Tune batch size to memory: a join materializes the matched pairs in memory, so very large left frames spike resident set size. The memory-management patterns in reducing memory overhead for 10,000-point batch routing apply directly — chunk the input and concatenate metrics, rather than joining one monolithic frame.
• Simplify polygon vertices: high-vertex catchments (drive-time isochrones can carry thousands of vertices) slow the exact predicate. A small-tolerance simplify() on the polygons cuts predicate cost without moving boundaries materially.
• Push the join into the database for large static catchments: a GiST-indexed PostGIS query streams results without loading the full point set into Python, which matters once transaction tables reach tens of millions of rows.

Log the spatial-index build time and the join duration each run; a regression against baseline usually signals dataset bloat or a dropped index, not a code change.

Validation and QA gates #

Run these checks before publishing metrics to any downstream consumer. They are deterministic and belong in the orchestration layer (Apache Airflow, Prefect, or GitHub Actions) as fail-fast gates:

1. Pre-flight validity: confirm both layers carry a CRS, geometry validity is at or above 99%, and record counts match expectations. Fail the run otherwise.
2. Coverage threshold: compute the null-match rate after the join. If it exceeds 5%, alert and route unmatched points to a quarantine table for manual review or nearest-neighbor reassignment rather than dropping them.
3. One-to-one attribution: assert that no transaction id appears under more than one catchment — the canary for an unbroken intersects tie.
4. Bounding-box sanity: verify the joined points fall inside the expected national or regional envelope, catching coordinate inversion and unit errors.
5. Performance baseline: record index build and join duration; flag degradation beyond the rolling baseline.

Coverage and topology gates are necessary but not sufficient — they confirm the join ran cleanly, not that it attributed points to the correct catchments. Confirm attribution accuracy against surveyed reality using the framework in validating spatial join accuracy with ground truth.

Integration notes #

The catchment_metrics table is the join key for the rest of the demographic pipeline. Once points are mapped to catchments, the next stage attaches socioeconomic profiles drawn from automated census feeds: block-group estimates are aggregated to match each catchment footprint, the workflow detailed in syncing US Census ACS data via API. Where a catchment straddles partially covered block groups, area-proportional interpolation handles the join of ACS 5-year estimates to custom trade area polygons.

From there, analysts apply business-logic weights to isolate high-propensity segments and normalize for household size or income, as covered in weighting demographic variables for target audiences. On successful aggregation, publish catchment_metrics to the warehouse and trigger the weighting stage — at which point the join has done its job: every point is attributed deterministically, every metric is auditable, and the catchment surface is ready for site scoring.

Related #

• Fixing sliver polygons in spatial join operations
• Syncing US Census ACS data via API
• Weighting demographic variables for target audiences
• Validating spatial join accuracy with ground truth

← Back to Demographic Data Integration & Spatial Joins