Spatial Overlay Analysis

When a planning commission adopts a new overlay district, the question every downstream system needs answered is deceptively simple: which parcels does this actually touch, and by how much? Getting it wrong is expensive in both directions. Miss a parcel that a floodplain or transit overlay clips by three percent and an underwriting model prices a property as if nothing changed; flag a parcel that only shares a boundary line with the new district and you generate a false notice that erodes trust in the whole feed. Spatial overlay analysis is the layer that turns a changed zoning polygon into a defensible, area-weighted list of affected parcels — and it fails in quiet, systematic ways when the predicate is loose, the coordinate reference system is geographic instead of projected, or the intersection leaves behind slivers thinner than a survey error. This topic sits inside the Spatial Impact Analysis & Zoning Change Detection framework and is the computational heart of it: it consumes changed geometries and emits the impacted-parcel table that everything else reasons over.

Prerequisites and operational context jump to heading

Overlay analysis is a settling pond: whatever contamination enters upstream concentrates here, because area math amplifies every geometric flaw. Before the patterns below produce trustworthy numbers, four things must already be true.

• Both layers share a projected CRS. Area computed on latitude/longitude degrees is meaningless — a “square degree” is not a unit of land. Districts and parcels must be reprojected to a projected system with real linear units before any .area call, following the CRS alignment strategies your organization standardizes on. For a single county, the local State Plane zone in feet or metres is ideal; for anything spanning multiple zones, an equal-area projection such as a US National Atlas Equal Area or a locally centred Albers keeps areas comparable.
• The parcel fabric is authoritative and clean. Overlay against a parcel layer riddled with self-intersections or duplicate geometries produces area fractions above 100% and negative slivers. Repair invalid geometry (Shapely make_valid) and deduplicate on the stable parcel identifier before the join.
• A spatial index exists on the parcel layer. A naive overlay compares every district against every parcel — O(n·m) — which is fatal at county scale. GeoPandas builds an STRtree lazily via .sindex, and for anything backed by a database the equivalent GiST index matters just as much; the trade-offs live in spatial database indexing and performance.
• The changed geometries are already isolated. Overlay does not detect what changed; it measures the consequence of a change someone else identified. The set of new or modified district polygons arrives from change detection and geometry diffing, so overlay never rescans the entire zoning map — only the delta.

If any of these is missing, fix it first. Running overlay on unprojected, unindexed, or unrepaired data does not fail loudly; it returns plausible-looking numbers that are wrong.

Architecture: the join, the intersection, and why they differ jump to heading

Two operations are easy to conflate and produce very different results: the spatial join and the spatial overlay. A spatial join (geopandas.sjoin) attaches attributes from one layer to another wherever a predicate holds, but it never changes geometry — a parcel that a district clips by one percent comes back with the district’s columns and its own full polygon intact. An overlay (geopandas.overlay) actually cuts geometry: with how="intersection" it returns the polygonal slivers of overlap, so a parcel clipped one percent comes back as a tiny fragment carrying both layers’ attributes. You need both. The join is a fast, index-backed filter that answers “which parcels are candidates?”; the overlay is the exact geometric operation that answers “how much of each candidate is inside?”

The predicate is where correctness is won or lost. intersects is the loosest test — it is true even when two polygons share only a single boundary point or edge, which is exactly the situation at every parcel that merely abuts a new district. If you treat intersects as “affected”, you will flag the entire ring of neighbours around a rezoned block. within and contains are the strict containment tests; overlaps requires interiors to intersect while neither fully contains the other. The correct pattern for impact analysis is to use intersects as the cheap candidate filter, then compute the actual intersection area and let a threshold — not the predicate — decide what counts as impacted. That way a shared edge contributes zero area and drops out naturally.

The third foundation of the architecture is the coordinate reference system, and it is worth being blunt about the failure it prevents. GeoPandas will happily run .overlay() and .area on data in EPSG:4326. The numbers come back as decimal fractions of a square degree, which vary with latitude and are useless for an area fraction. Reprojecting both layers to a projected CRS before the overlay — not after — is non-negotiable, because the intersection geometry must be computed in the same planar space the areas are measured in. Reproject once, up front, and keep every geometry column in that CRS for the whole pipeline.

Production implementation jump to heading

The function below is the core of a county-scale impact run. It reprojects both layers to a supplied projected CRS, repairs geometry, uses the spatial index to pre-filter candidate pairs, runs a true overlay intersection, computes each parcel’s affected-area fraction against its original area, drops slivers below a configurable threshold, and returns a tidy impacted-parcel frame with enough provenance to audit. It is defensive about the failure modes that actually occur in municipal data: geographic CRS input, invalid polygons, empty overlays, and division by zero on degenerate parcels.

from __future__ import annotations

import logging
from dataclasses import dataclass

import geopandas as gpd
import pandas as pd
from shapely import make_valid

logger = logging.getLogger("overlay_impact")


@dataclass(frozen=True)
class OverlayConfig:
    # A PROJECTED CRS with linear units. State Plane feet for one county,
    # an equal-area projection for multi-zone extents. Never a geographic CRS.
    working_epsg: int = 2263          # NY State Plane Long Island (US ft)
    sliver_frac: float = 0.005        # <0.5% of a parcel affected => ignore
    sliver_area_min: float = 5.0      # absolute floor in CRS units (ft^2)
    parcel_id_col: str = "parcel_id"
    district_id_col: str = "district_id"


def _prepare(layer: gpd.GeoDataFrame, epsg: int, name: str) -> gpd.GeoDataFrame:
    """Reproject to the working CRS and repair invalid geometry."""
    if layer.crs is None:
        raise ValueError(f"{name} has no CRS set; refusing to guess.")
    if layer.crs.is_geographic and layer.crs.to_epsg() == epsg:
        # working_epsg must be projected; a geographic target is a config bug.
        raise ValueError(f"working CRS {epsg} is geographic; area math needs a projected CRS.")

    projected = layer.to_crs(epsg=epsg)
    # Repair self-intersections / bowties so overlay does not raise or emit noise.
    invalid = ~projected.geometry.is_valid
    if invalid.any():
        logger.warning("%s: repairing %d invalid geometries", name, int(invalid.sum()))
        projected.loc[invalid, projected.geometry.name] = projected.loc[
            invalid, projected.geometry.name
        ].apply(make_valid)
    # Drop anything that is still empty or non-polygonal after repair.
    projected = projected[~projected.geometry.is_empty & projected.geometry.notna()]
    return projected


def compute_parcel_impact(
    parcels: gpd.GeoDataFrame,
    changed_districts: gpd.GeoDataFrame,
    cfg: OverlayConfig = OverlayConfig(),
) -> gpd.GeoDataFrame:
    """Return impacted parcels with the fraction of each parcel's area affected.

    parcels           authoritative parcel fabric (any CRS with a defined CRS)
    changed_districts new / modified zoning district polygons (the delta only)
    """
    parcels = _prepare(parcels, cfg.working_epsg, "parcels")
    districts = _prepare(changed_districts, cfg.working_epsg, "districts")

    if parcels.empty or districts.empty:
        logger.info("Nothing to overlay (parcels=%d, districts=%d)", len(parcels), len(districts))
        return parcels.iloc[0:0].assign(affected_frac=pd.Series(dtype="float64"))

    # Record each parcel's TRUE area before it is cut, keyed on parcel id.
    # The fraction must divide by the whole parcel, not the clipped fragment.
    parcels = parcels.copy()
    parcels["parcel_area"] = parcels.geometry.area
    degenerate = parcels["parcel_area"] <= 0
    if degenerate.any():
        logger.warning("Dropping %d parcels with zero/negative area", int(degenerate.sum()))
        parcels = parcels[~degenerate]

    # 1) Cheap, index-backed candidate filter. sjoin uses the STRtree on
    #    parcels.sindex; 'intersects' includes edge-touchers, which the area
    #    threshold will remove later — that is intentional.
    candidates = gpd.sjoin(
        parcels[[cfg.parcel_id_col, "parcel_area", parcels.geometry.name]],
        districts[[cfg.district_id_col, districts.geometry.name]],
        how="inner",
        predicate="intersects",
    )
    if candidates.empty:
        logger.info("No parcel intersects any changed district.")
        return parcels.iloc[0:0].assign(affected_frac=pd.Series(dtype="float64"))

    # 2) True geometric intersection ONLY on the candidate parcels/districts.
    #    Restricting the overlay to matched ids keeps it O(candidates), not O(n*m).
    hit_parcels = parcels[parcels[cfg.parcel_id_col].isin(candidates[cfg.parcel_id_col])]
    hit_districts = districts[districts[cfg.district_id_col].isin(candidates[cfg.district_id_col])]
    pieces = gpd.overlay(hit_parcels, hit_districts, how="intersection", keep_geom_type=True)
    if pieces.empty:
        logger.info("Overlay produced no polygonal overlap (all edge-touch).")
        return parcels.iloc[0:0].assign(affected_frac=pd.Series(dtype="float64"))

    # 3) Area of each overlap fragment, then the fraction of the whole parcel.
    pieces["overlap_area"] = pieces.geometry.area
    # A parcel may split across several districts; sum the fragments per parcel.
    per_parcel = (
        pieces.groupby(cfg.parcel_id_col, as_index=False)
        .agg(overlap_area=("overlap_area", "sum"), parcel_area=("parcel_area", "first"))
    )
    per_parcel["affected_frac"] = per_parcel["overlap_area"] / per_parcel["parcel_area"]
    # Clamp floating-point overshoot from repaired geometry (e.g. 1.0000003).
    per_parcel["affected_frac"] = per_parcel["affected_frac"].clip(upper=1.0)

    # 4) Drop slivers: both a relative fraction AND an absolute-area floor, so a
    #    huge parcel touched by a 2 ft^2 shard is not "0.6% affected => keep".
    keep = (
        (per_parcel["affected_frac"] >= cfg.sliver_frac)
        & (per_parcel["overlap_area"] >= cfg.sliver_area_min)
    )
    dropped = int((~keep).sum())
    if dropped:
        logger.info("Dropped %d sliver overlaps below threshold", dropped)
    impacted = per_parcel[keep].copy()

    # 5) Re-attach parcel attributes and geometry for the output layer.
    out = parcels.merge(
        impacted[[cfg.parcel_id_col, "overlap_area", "affected_frac"]],
        on=cfg.parcel_id_col,
        how="inner",
    )
    out = gpd.GeoDataFrame(out, geometry=parcels.geometry.name, crs=parcels.crs)
    out["affected_pct"] = (out["affected_frac"] * 100).round(2)
    logger.info("Impact run: %d parcels affected across %d districts",
                len(out), districts[cfg.district_id_col].nunique())
    return out.sort_values("affected_frac", ascending=False)

Several choices in this code are load-bearing. The candidate filter and the overlay are deliberately split so the expensive geometric cut runs only on parcels that could plausibly be affected — the STRtree behind sjoin turns a quadratic scan into a near-linear one. The affected fraction divides by parcel_area, captured before the overlay cuts the geometry, because dividing the fragment by itself would always yield 1.0. Sliver removal uses two gates at once — a relative fraction and an absolute floor — because either alone is exploitable: a relative-only threshold keeps meaningless shards of enormous parcels, while an absolute-only threshold discards real impacts on tiny urban lots. And every early return preserves the output schema (affected_frac column present, correct CRS) so downstream consumers never branch on an empty-versus-populated frame.

Edge cases and gotchas jump to heading

Overlay math breaks in specific, reproducible ways. Budget for these before they reach an underwriting report.

• Geographic CRS silently accepted. gdf.area on EPSG:4326 returns a number, so nothing raises — but it is square degrees. Guard against it explicitly (the _prepare check above) rather than trusting callers to reproject.
• keep_geom_type off by default in older GeoPandas. An intersection of two polygons can yield a GeometryCollection containing stray lines or points along shared edges. Without keep_geom_type=True, those non-polygonal parts inflate the fragment count and can crash a .area sum. Force polygonal output.
• Fractions above 1.0. Repaired or slightly overlapping parcel geometry produces overlap areas marginally larger than the parent, so an honest division yields 1.0000004. Clip to 1.0; do not let a 100.00004% affected figure escape into a compliance artifact.
• Duplicate parcel ids across the fabric. If the same parcel_id appears twice (a common annexation artifact), the groupby sums both, over-reporting impact. Deduplicate on the identifier plus geometry before the run.
• District polygons that overlap each other. When two changed overlays cover the same ground, summing per-parcel fragments double-counts the shared area. If overlays can stack, dissolve the district layer (unary_union or dissolve) first so the impact reflects unique affected land.
• Empty overlay treated as an error. A refresh where nothing intersects is a valid, common result, not a failure. Return an empty frame with the right columns so the alerting stage reads “zero impacted” cleanly instead of throwing.
• Precision-grid snapping. GEOS may drop micro-fragments below its precision grid, so overlap areas can differ by nanometres between runs on identical input. Round affected_pct and never assert bit-exact equality across runs.

Overlay predicate reference jump to heading

The predicate you pass to sjoin (and the mental model you apply to overlay results) determines which parcels survive the filter. Choosing the wrong one is the most common source of over- and under-reporting.

The rule of thumb: filter with intersects, measure with an overlay, and decide with a threshold. Never use intersects alone as the definition of “impacted”, and never use touches to mean anything but adjacency.

Integration points jump to heading

Overlay analysis is a midstream stage with a clear contract on both sides. Its input is the isolated set of new and modified district polygons produced by change detection and geometry diffing — feeding it only the delta, rather than the full zoning map, is what keeps each run bounded and cheap. Its output is the impacted-parcel table, and the primary consumer of that table is zoning change alerting, which turns each row — a parcel id, its percent affected, and the district that hit it — into a notice, a webhook, or a watchlist update. The stability of that output schema matters more than any internal detail: as long as overlay emits the same columns with the same CRS and the same sliver semantics, the diffing and alerting stages can evolve independently. For the exact end-to-end walkthrough of one rezoning event, see computing parcel-level zoning impact with GeoPandas overlay.

Compliance and audit artifacts jump to heading

For PropTech underwriting and regulatory review, an impact number is only as good as its provenance. A parcel reported as “37% inside the new inclusionary-housing overlay” must be reconstructable months later, from the same inputs, to the same figure. Every overlay run should emit:

• A run manifest recording the working EPSG, the sliver thresholds, the source vintages of both the parcel fabric and the changed-district set, and the GEOS/GeoPandas versions — the geometry library version matters because precision handling changes between releases.
• Per-parcel provenance linking each impacted parcel back to the specific district id that hit it and the overlap area in CRS units, not just the percentage, so a reviewer can recompute the fraction independently.
• A reconciliation total — the sum of overlap areas against the total area of the changed districts — which surfaces double-counting from overlapping overlays or duplicate parcels as a mismatch.
• A dropped-sliver ledger capturing every fragment removed by the threshold with its parcel id and area, so a challenged “not affected” determination can be defended with the exact number that fell below the line.

These artifacts pair naturally with the schema and lineage discipline that CRS alignment strategies establish upstream, extending the provable chain from raw municipal publication all the way to an underwriting-grade impact determination.

FAQ jump to heading

Related jump to heading

• Computing parcel-level zoning impact with GeoPandas overlay — the end-to-end walkthrough of one rezoning event
• Change Detection & Geometry Diffing — isolates the changed districts that feed this stage
• Zoning Change Alerting — consumes the impacted-parcel list this stage emits
• CRS Alignment Strategies — the projected-CRS contract every overlay depends on
• Spatial Database Indexing & Performance — the indexing that keeps candidate filtering near-linear

Up: Spatial Impact Analysis & Zoning Change Detection