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

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.

Three symptoms, and the check that identifies each cause A three-row by four-column matrix of overlay symptoms. Rows are: the result has more rows than the input, the result has fewer rows than expected, and areas do not sum to the parcel area. Columns are the likely cause, the check that confirms it, the fix, and whether the symptom would be visible in a row-count sanity check. Row inflation comes from a parcel spanning several districts, which is correct behaviour rather than a bug. Row loss comes from an inner join dropping parcels with no overlap. Area mismatch comes from a CRS in degrees, where area is meaningless. Three symptoms, and the check that identifies each cause likely cause check fix caught by a row-count check? More rows out than in parcel spans districts group by parcel_id, count expected — keep the split yes, misleadingly Fewer rows than expected inner join dropped no-overlap anti-join the inputs use how="left" yes Areas do not sum to the parcel area computed in degrees check crs.is_geographic project to an equal-area CRS no correct or actionable ambiguous silently wrong
The first row is not a bug at all: a parcel spanning three districts should produce three rows, and "fixing" it by de-duplicating destroys the split-zoning information the analysis exists to produce. The third row is the dangerous one — computing area in a geographic CRS returns a number, in square degrees, that means nothing.
  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.

Five-step flow for parcel-level zoning impact with GeoPandas overlay A linear five-step flow. Step one reprojects both the parcel fabric and the overlay district to an equal-area or State Plane CRS and captures each parcel's true area. Step two runs a GeoPandas overlay intersection to produce overlap fragments. Step three computes the area fraction of each parcel affected. Step four drops slivers below a threshold. Step five attaches the impact percentage to the impacted parcels. A verification note beneath checks that summed overlap area reconciles against the district area. 1 · Reproject equal-area / ft save true area 2 · Overlay how=intersection keep_geom_type 3 · Fraction overlap / parcel clip to 1.0 4 · Drop slivers frac & area gates log dropped 5 · Attach % area impacted Verify: Σ overlap area ≈ district area mismatch => overlapping districts or duplicate parcels
Reproject and capture true area first; overlay, weight, and threshold in sequence; then reconcile summed overlap against district area before trusting the output.

Step-by-step implementation

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

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

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

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

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

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

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

Area error introduced by computing in the wrong CRS Lollipop chart of the error in a computed parcel area under five CRS choices, expressed as a percentage of the true area for a parcel at 40 degrees north. Computing in a local equal-area projection is the reference at zero. A state plane projection introduces 0.02 per cent. UTM at the zone centre introduces 0.04 per cent, and at the zone edge 0.16 per cent. Web Mercator introduces about 70 per cent at this latitude. Computing in geographic degrees does not produce a percentage error at all because the result is in square degrees and is not an area. Area error introduced by computing in the wrong CRS 0 20 40 60 80 1% — beyond here density figures are wrong local equal-area projection 0% state plane (correct zone) 0.02% UTM at the zone centre 0.04% UTM at the zone edge 0.16% Web Mercator at 40°N 70% Area error (% of true area) Geographic degrees are omitted: gdf.area returns square degrees, which is not an area at all.
Web Mercator is the one that surprises people: it is the default of every web map, and at 40°N it inflates area by roughly 70% — a density figure computed there is not slightly off, it is wrong. Degrees are worse than wrong: gdf.area returns square degrees, a number with no physical meaning that will happily be summed and reported.
  • 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

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

Why divide by parcel_area captured before the overlay?

The overlay cuts each parcel down to the overlap fragment, so dividing the fragment by its own post-overlay area always yields 1.0. Capturing each parcel’s whole-polygon area in Step 1, before the geometry is clipped, gives you the correct denominator — the fraction of the entire parcel that now falls inside the district.

What equal-area or projected CRS should I use for one corridor?

For a corridor inside a single State Plane zone, that zone in US feet (for example EPSG:2263 for Long Island) is exact and convenient. For an extent spanning two zones, use an equal-area projection such as a locally centred Albers so areas stay comparable across the seam. The only hard rule is that the CRS must be projected with linear units, never geographic degrees.

How do I stop overlay from returning hundreds of tiny slivers?

The slivers are real but meaningless overlaps caused by boundary disagreement between the parcel and district sources. Apply a two-part threshold after computing fractions: a relative gate (below roughly 0.5% of the parcel) and an absolute area floor in CRS units. Requiring the fragment to clear a small noise floor before either gate applies removes precision artifacts without discarding genuine partial impacts.

My overlap areas sum to more than the district area — what happened?

Either the overlay district layer self-overlaps, so shared ground is counted once per piece, or the parcel fabric contains duplicate ids that the group-by sums twice. Dissolve the district layer before the overlay and deduplicate parcels on their identifier plus geometry, then the reconciliation total should fall back in line with the district’s own area.