Performing Point-in-Polygon Joins for Store Catchments

In retail site selection automation, assigning transaction points, prospect centroids, or competitor locations to predefined trade area boundaries is a fundamental spatial operation. Performing Point-in-Polygon Joins for Store Catchments serves as the geometric backbone of location intelligence pipelines, directly feeding downstream revenue forecasting, lease evaluation, and network optimization models. This operation anchors the broader Demographic Data Integration & Spatial Joins workflow, where raw coordinate streams are systematically enriched with aggregated socioeconomic indicators. Production deployments require strict attention to coordinate reference system (CRS) alignment, topology validation, and predicate selection to prevent silent data loss or misattribution.

Configuration & Execution Parameters #

The spatial join evaluates whether a coordinate pair intersects the topological boundary of a polygonal catchment. Configuration choices dictate pipeline reliability:

• CRS: Always enforce a projected CRS (e.g., a regional UTM zone like EPSG:32617 for eastern US) before executing joins to eliminate angular distortion and ensure accurate spatial indexing.
• Predicate: Use predicate="within" to strictly match points inside polygon boundaries. When catchments overlap—common in multi-format retail portfolios—switch to predicate="intersects" and implement a deterministic tie-breaker (e.g., shortest Euclidean distance to the anchor store centroid) to prevent duplicate attribution.
• Fallback: Apply sjoin_nearest as a fallback for points that fall just outside a catchment due to GPS drift or boundary precision limits.

flowchart TD
    P["Transaction points"] --> AL["Align to projected CRS"]
    C["Catchment polygons<br/>make_valid"] --> AL
    AL --> SJ["sjoin predicate = within"]
    SJ --> U{"Unmatched points?"}
    U -->|"yes"| NN["sjoin_nearest fallback<br/>max_distance = 5 km"]
    U -->|"no"| AGG["Aggregate per catchment<br/>count · revenue · avg ticket"]
    NN --> AGG

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:32617"
) -> pd.DataFrame:
    """
    Joins transaction points to retail catchment polygons.
    Applies nearest-neighbor fallback for unmatched edge cases.
    """
    # 1. Load and repair catchment topologies
    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 identical projected CRS for spatial indexing
    if catchments.crs is None:
        raise ValueError("Catchment dataset lacks CRS definition. Assign before join.")
    catchments = catchments.to_crs(target_crs)
    points_gdf = points_gdf.to_crs(target_crs)

    # 4. Execute spatial join with explicit predicate
    # Reference: https://geopandas.org/en/stable/docs/reference/api/geopandas.sjoin.html
    joined = gpd.sjoin(points_gdf, catchments, how="left", predicate="within")

    # 5. Nearest-neighbor fallback for boundary/edge cases
    unmatched_mask = joined["index_right"].isna()
    if unmatched_mask.any():
        logging.warning(
            "%d points fell outside catchments. Applying nearest-neighbor fallback.",
            unmatched_mask.sum()
        )
        # sjoin_nearest operates on the original unmatched points (no "index_right" column 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

Debugging & Topology Management #

Boundary precision and sliver geometries frequently corrupt join outputs. Points from mobile GPS logs, POS systems, or third-party APIs often carry sub-meter drift, causing legitimate customer locations to fall outside drive-time isochrones. When working with manually digitized trade zones, validate polygon topology before joining to eliminate self-intersections and sliver artifacts that trigger false negatives. For detailed remediation workflows, consult Fixing sliver polygons in spatial join operations.

Always audit join coverage rates post-execution. A sudden drop below 95% typically indicates CRS mismatch, corrupted GeoJSON, or upstream coordinate inversion. Log geometry failure counts for upstream data engineering review.

Downstream Integration & Automation Triggers #

The output of a successful join is the primary key for demographic enrichment and predictive modeling. Once points are mapped to catchments, pipelines aggregate transaction volumes and attach socioeconomic profiles pulled via automated census feeds. This workflow aligns directly with Syncing US Census ACS Data via API, where block-group-level indicators are spatially aggregated to match catchment footprints. Following aggregation, analysts apply demographic weighting to isolate high-propensity segments and normalize for household size or income brackets, as detailed in Weighting Demographic Variables for Target Audiences.

To operationalize this process, embed the join routine in a scheduled orchestration framework (Apache Airflow, Prefect, or GitHub Actions) with the following automation triggers and validation gates:

1. Pre-flight Validation: Verify CRS consistency, geometry validity, and record counts. Fail fast if geometry validity drops below 99%.
2. Join Coverage Threshold: Trigger an alert if the null-match rate exceeds 5%. Route unmatched points to a quarantine table for manual review or nearest-neighbor reassignment.
3. Performance Monitoring: Log spatial index build times and join execution duration. Degradation beyond baseline indicates dataset bloat or missing spatial indexing.
4. Downstream Handoff: Upon successful aggregation, publish the catchment_metrics table to the data warehouse and trigger the demographic weighting pipeline.

By standardizing predicate selection, enforcing topology validation, and embedding automated quality gates, retail planners and location intelligence teams ensure spatial assignments remain deterministic, auditable, and production-ready across all market rollouts.