Computing parcel-level zoning impact with GeoPandas overlay

A transit-oriented-development overlay is adopted over a six-block corridor, and by Monday morning your underwriting team needs one thing: the exact list of parcels it touches, each with the percentage of its area now inside the new district. You load the parcel fabric and the overlay polygon, run gpd.overlay(parcels, overlay, how="intersection"), and the result looks reasonable — until someone notices the report lists 900 parcels for a corridor that visibly contains about 140, half the “affected” fractions are 0.03%, and one lot reads 340% affected. Every one of those symptoms traces to the same three root causes: a geographic CRS under the area math, no sliver threshold, and dividing by the wrong denominator. This guide is the concrete, reproducible fix for that scenario — the linear how-to behind the broader spatial overlay analysis topic. We will reproject to an equal-area or State Plane CRS, run a clean intersection, compute the fraction against each parcel’s real area, drop the slivers, and verify the totals reconcile.

Diagnosis: three symptoms, three root causes jump to heading

Reproduce the broken run before changing anything, so each fix can be measured against a known-bad baseline.

import geopandas as gpd

parcels = gpd.read_file("parcels.gpkg")          # EPSG:4326, degrees
overlay = gpd.read_file("tod_overlay.geojson")   # EPSG:4326, degrees

hits = gpd.overlay(parcels, overlay, how="intersection")
hits["frac"] = hits.geometry.area / hits.geometry.area   # <-- denominator bug
print(len(hits), "fragments")
print(hits["frac"].describe())
print("CRS:", parcels.crs.is_geographic)

Three things in that snippet are wrong, and each maps to a visible symptom in the report.

1. CRS: True — the area is in square degrees. With both layers in EPSG:4326, geometry.area returns a fraction of a square degree. It runs without error, which is exactly why it is dangerous: the numbers are internally consistent and completely meaningless as land area. Any per-parcel fraction built on them is noise.
2. No sliver filter — the 900-versus-140 gap. overlay(..., how="intersection") returns every polygonal overlap, including the paper-thin fragments where a parcel’s boundary runs a few centimetres inside the district edge because the two survey sources disagree. Those coordinate-precision slivers are real polygons with a tiny positive area, so they show up as “affected” parcels at 0.02–0.05%.
3. The 340% parcel — the denominator bug. Dividing the fragment area by itself always yields 1.0; the real intent is to divide the overlap by the whole parcel. When the code instead joins fragments back and divides by an already-clipped or duplicated geometry, a parcel split across overlapping district pieces sums past its own area and reports an impossible percentage.

The fix addresses all three in order: project first, measure against the true parent area, then threshold.

Step-by-step implementation jump to heading

Each step fixes one root cause and can be run and inspected on its own.

Step 1 — Reproject to an equal-area CRS and capture true area jump to heading

Pick a projected CRS with linear units and reproject both layers into it before touching .area. For a single county, the local State Plane zone in US feet is exact enough; for a corridor spanning two zones, an equal-area projection keeps areas comparable. Critically, record each parcel’s whole-polygon area now, while the geometry is still intact — this is the denominator every later fraction divides by. The choice of projection is the same discipline covered under CRS alignment strategies; here it is the difference between a real percentage and square-degree noise.

import geopandas as gpd

WORKING_EPSG = 2263  # NY State Plane Long Island, US survey feet

def prepare(layer: gpd.GeoDataFrame, epsg: int) -> gpd.GeoDataFrame:
    if layer.crs is None:
        raise ValueError("layer has no CRS; refusing to assume one")
    proj = layer.to_crs(epsg=epsg)
    if proj.crs.is_geographic:
        raise ValueError(f"EPSG:{epsg} is geographic; area needs a projected CRS")
    # Repair before measuring so a bowtie parcel does not report a wrong area.
    proj["geometry"] = proj.geometry.make_valid()
    return proj[~proj.geometry.is_empty & proj.geometry.notna()]

parcels = prepare(gpd.read_file("parcels.gpkg"), WORKING_EPSG)
overlay = prepare(gpd.read_file("tod_overlay.geojson"), WORKING_EPSG)

parcels = parcels.copy()
parcels["parcel_area"] = parcels.geometry.area          # the true denominator
parcels = parcels[parcels["parcel_area"] > 0]           # drop degenerate lots

Step 2 — Run the overlay intersection jump to heading

With both layers in the same projected space, run the intersection. Pass keep_geom_type=True so an overlap that GEOS returns as a GeometryCollection (a polygon plus a stray edge line) is reduced to its polygonal parts instead of polluting the area sum. If several overlay districts can stack on the same ground, dissolve them first so shared area is not counted twice.

# If the overlay layer can self-overlap, collapse it to unique coverage first.
overlay_dissolved = overlay.dissolve().reset_index(drop=True)

pieces = gpd.overlay(
    parcels[["parcel_id", "parcel_area", "geometry"]],
    overlay_dissolved[["geometry"]],
    how="intersection",
    keep_geom_type=True,
)
if pieces.empty:
    raise SystemExit("No parcel overlaps the district — verify inputs and CRS.")

Step 3 — Compute the area fraction against the parent parcel jump to heading

Now measure each fragment and divide by the parcel’s original area, not the fragment. A parcel can appear in several fragments if the dissolved district still wraps around it, so sum the fragments per parcel first, then divide once.

pieces["overlap_area"] = pieces.geometry.area

per_parcel = (
    pieces.groupby("parcel_id", 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"]
# Clip float overshoot from repaired geometry (e.g. 1.0000006 -> 1.0).
per_parcel["affected_frac"] = per_parcel["affected_frac"].clip(upper=1.0)

Step 4 — Drop slivers with two thresholds jump to heading

Remove the coordinate-precision fragments that inflated the parcel count. Use a relative fraction and an absolute area floor together: the relative gate ignores boundary-noise slivers, and the absolute gate stops a huge parcel’s tiny shard from sneaking through on a low percentage. Log what you drop so a “not affected” call is defensible.

SLIVER_FRAC = 0.005      # < 0.5% of the parcel
SLIVER_AREA = 25.0       # and < 25 sq ft absolute

keep = (per_parcel["affected_frac"] >= SLIVER_FRAC) | (per_parcel["overlap_area"] >= SLIVER_AREA)
# Require the fragment to clear the noise floor before either gate applies.
keep &= per_parcel["overlap_area"] >= SLIVER_AREA * 0.2
dropped = per_parcel[~keep]
print(f"dropped {len(dropped)} sliver overlaps")
impacted = per_parcel[keep].copy()

Step 5 — Attach the impact percentage to the parcels jump to heading

Merge the surviving fractions back onto the full parcel geometry and attributes to produce the deliverable: an impacted-parcel layer carrying affected_pct, ready for the spatial overlay analysis pipeline’s downstream stages within the Spatial Impact Analysis & Zoning Change Detection framework.

out = parcels.merge(
    impacted[["parcel_id", "overlap_area", "affected_frac"]],
    on="parcel_id", how="inner",
)
out = gpd.GeoDataFrame(out, geometry="geometry", crs=parcels.crs)
out["affected_pct"] = (out["affected_frac"] * 100).round(2)
out = out.sort_values("affected_frac", ascending=False)
out.to_file("impacted_parcels.gpkg", driver="GPKG")
print(out[["parcel_id", "affected_pct"]].head(10).to_string(index=False))

Verification and testing jump to heading

Prove the fix corrected the numbers rather than merely changing them.

• Reconcile the areas. The sum of overlap_area across all impacted parcels should be close to the dissolved district’s own area (overlay_dissolved.area.sum()), minus any district ground that falls on roads or water outside the parcel fabric. A large shortfall means missing parcels; an excess means double-counting from undissolved overlays or duplicate parcel ids.
• Sanity-check the count. For the six-block corridor, the impacted count should land near the visual estimate (~140), not in the hundreds. If it is still inflated, raise the sliver gates; if parcels you can see inside the district are missing, they are being dropped as slivers — lower the absolute floor.
• No fraction over 100%. Assert out["affected_pct"].max() <= 100.0. A failure means the denominator or the dissolve step is still wrong.
• Round-trip a known parcel. Take one lot you can measure by hand from the plat, confirm its affected_pct matches within a percent or two, and confirm a fully-inside parcel reports 100% and a merely-adjacent one is absent entirely.
• Determinism. Re-run the whole script and diff the output; identical inputs must yield identical affected_pct values. Drift signals non-deterministic ordering or an unpinned GEOS precision setting.

Failure recovery jump to heading

When a run misbehaves mid-corridor, keep it reproducible and recoverable.

• Empty overlay is not a crash. If nothing intersects, write an empty output with the same schema and exit cleanly, so the alerting stage reads “zero impacted” instead of erroring. Verify the CRS match first — an empty result on data you know overlaps almost always means the two layers were never in the same CRS.
• GeometryCollection errors on .area. If a TopologyException or a mixed-type error surfaces, confirm keep_geom_type=True is set and re-run make_valid on the offending inputs before the overlay.
• Fraction anomalies. A group of parcels above 100% points straight at the denominator or an undissolved district; recompute parcel_area in Step 1 and dissolve in Step 2 before touching thresholds.
• Quarantine, never guess. A parcel whose geometry cannot be repaired is written to a failed_geometry/ partition with its id and the exception, then skipped, so one bad polygon never voids the corridor’s impact table.

Frequently asked questions jump to heading

Related jump to heading

• Parent topic: Spatial Overlay Analysis
• Section overview: Spatial Impact Analysis & Zoning Change Detection
• CRS Alignment Strategies — choosing and pinning the projected CRS this run depends on
• Change Detection & Geometry Diffing — where the changed district polygon comes from
• Zoning Change Alerting — the consumer of the impacted-parcel table this guide produces