Diffing parcel geometries to detect zoning boundary changes

You pull this month’s parcel export from a county GIS portal, compare it against last month’s copy, and your change report claims 4,812 parcels were modified overnight. The planning department confirms three rezones and one lot-line adjustment happened that month. So where did the other 4,808 “changes” come from? Nothing was rezoned — the county’s export server simply re-serialized the layer: a new coordinate-precision setting nudged the seventh decimal place, the ring winding flipped, and the vertex that starts each polygon moved. Every one of those parcels is now byte-different and geometrically identical, and a naive comparison cannot tell the four real changes apart from the thousands of cosmetic ones. This guide is a focused, runnable companion to the broader Change Detection & Geometry Diffing topic within the Spatial Impact Analysis & Zoning Change Detection framework, and it answers one narrow question: how do you diff two parcel exports so that only genuine boundary and zoning changes survive?

Diagnosis: why identical parcels report as changed jump to heading

Before writing any comparison, reproduce the false positives and name their cause. Load both exports and run the comparison a naive pipeline would run — a direct equality test on the raw geometry — then inspect a parcel you know did not change. A representative trace from a Wisconsin county export looks like this:

import geopandas as gpd

t0 = gpd.read_file("parcels_2026_05.gpkg")
t1 = gpd.read_file("parcels_2026_06.gpkg")

a = t0.set_index("parcel_id").loc["0708-241-0913-4"]
b = t1.set_index("parcel_id").loc["0708-241-0913-4"]

print(a.geometry.equals_exact(b.geometry, tolerance=0))  # False  -> flagged as changed
print(a.geometry.equals(b.geometry))                     # True   -> topologically identical
print(list(a.geometry.exterior.coords)[:2])
# [(571043.2100001, 289114.8700000), (571051.4400000, 289120.1300000)]
print(list(b.geometry.exterior.coords)[:2])
# [(571051.44, 289120.13), (571061.07, 289126.55)]   # different start vertex + rounded

The two rows describe the same lot. equals_exact with zero tolerance returns False because the coordinate strings differ; the topological equals() returns True because the shapes are the same. Three distinct, recognizable noise sources produce this:

1. Vertex reordering and ring flips. The exporter starts the polygon at a different vertex and winds it the other way. The coordinate sequence changes while the boundary does not.
2. Precision jitter. A change in the source’s output precision, or a reprojection performed on ingest, perturbs coordinates in the sixth or seventh decimal. Sub-centimeter noise, invisible on a map, breaks any exact-string comparison.
3. CRS drift. If the two exports arrive in different projections — one in the state’s meters-based grid, one in latitude/longitude — the same parcel occupies completely different coordinate values, and a diff computed across them reports every parcel as having moved kilometers.

A real change looks different from all three: a boundary adjustment changes the area of the symmetric difference by a meaningful amount, and a rezone changes an attribute while leaving the geometry untouched. The fix is a pipeline that neutralizes the three noise sources, then measures what remains against a threshold.

Step-by-step implementation jump to heading

Each step neutralizes one failure mode from the diagnosis. Compose them in order — the sequence is load-bearing, because snapping a geometry that has not been reprojected snaps in the wrong units.

Step 1 — Reproject both exports into one CRS jump to heading

CRS drift is the noise source that produces the largest phantom movements, so eliminate it first. Force both exports into a single projected reference system measured in linear units, because every later threshold is expressed in meters. This is the parcel-diff application of the broader CRS alignment strategies discipline; here the rule is simply that the diff refuses to proceed until both sides share one metric CRS.

import geopandas as gpd

WORKING_CRS = "EPSG:3071"   # Wisconsin Transverse Mercator, meters

def to_working_crs(gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    if gdf.crs is None:
        raise ValueError("Export has no CRS; cannot diff safely.")
    out = gdf.to_crs(WORKING_CRS)
    if out.crs.is_geographic:
        raise ValueError("Working CRS must be projected (meters), not degrees.")
    return out

t0 = to_working_crs(gpd.read_file("parcels_2026_05.gpkg"))
t1 = to_working_crs(gpd.read_file("parcels_2026_06.gpkg"))

Step 2 — Normalize geometry to canonical form jump to heading

Vertex reordering and ring flips vanish once both geometries are rewritten into a canonical form. Shapely’s normalize() fixes ring orientation, the interior-ring order, and the starting vertex so that two descriptions of the same boundary become directly comparable. Repair invalid geometry first, because normalize() can raise on a self-intersecting polygon.

from shapely.geometry.base import BaseGeometry

def canonicalize(geom: BaseGeometry) -> BaseGeometry | None:
    """Return a normalized geometry, or None if it cannot be repaired."""
    if geom is None or geom.is_empty:
        return None
    if not geom.is_valid:
        geom = geom.buffer(0)                 # close bowties / self-intersections
        if geom.is_empty or not geom.is_valid:
            return None
    return geom.normalize()

Step 3 — Snap coordinates to a tolerance grid jump to heading

Precision jitter survives normalization because a re-rounded coordinate is still a different coordinate. Snapping every vertex to a fixed grid collapses that jitter: pick a grid size coarser than your reprojection’s residual error but far finer than any real survey adjustment. At 5 cm, a seventh-decimal wobble disappears while a genuine two-meter lot-line move is untouched.

from shapely import set_precision

SNAP_GRID = 0.05   # meters; coarser than reprojection noise, finer than real moves

def snap(geom: BaseGeometry) -> BaseGeometry | None:
    snapped = set_precision(geom, SNAP_GRID)
    return None if snapped.is_empty else snapped

After snapping, hash the canonical well-known binary. If both parcels hash to the same value, the difference was pure noise and you can stop — no expensive geometry math required.

import hashlib

def geom_key(geom: BaseGeometry) -> str:
    return hashlib.sha256(geom.wkb).hexdigest()

Step 4 — Threshold the symmetric-difference area jump to heading

Only pairs whose hashes differ reach this stage, and even here a difference is not yet a change. Measure the symmetric difference — the area in one polygon or the other but not both — and the Hausdorff distance, the widest gap between the two boundaries. A real boundary adjustment clears an absolute area floor; a single displaced survey vertex that barely changes area is caught by the Hausdorff floor. Requiring one or the other keeps hairline slivers from firing while still catching narrow moves.

AREA_FLOOR = 1.0        # m^2; below this a footprint difference is noise
HAUSDORFF_FLOOR = 0.5   # m; a vertex displaced beyond this is real

def geometry_moved(a: BaseGeometry, b: BaseGeometry) -> tuple[bool, float, float]:
    area = a.symmetric_difference(b).area
    hd = a.hausdorff_distance(b)
    return (area >= AREA_FLOOR or hd >= HAUSDORFF_FLOOR), round(area, 3), round(hd, 3)

Step 5 — Compute the attribute delta and classify jump to heading

A rezone leaves the footprint untouched and changes the zoning code, so geometry thresholds will never catch it — you must compare attributes independently. Combine the two signals into one verdict per parcel: geometry moved, attributes changed, both, or neither.

WATCHED = ("zoning_code", "overlay", "land_use")

def classify(a_row, b_row) -> dict | None:
    attr_delta = {c: (a_row.get(c), b_row.get(c))
                  for c in WATCHED if a_row.get(c) != b_row.get(c)}

    ga, gb = canonicalize(a_row.geometry), canonicalize(b_row.geometry)
    if ga is None or gb is None:
        return {"class": "quarantine", "attr_delta": attr_delta}
    ga, gb = snap(ga), snap(gb)

    moved, area, hd = (False, 0.0, 0.0)
    if geom_key(ga) != geom_key(gb):
        moved, area, hd = geometry_moved(ga, gb)

    changed = bool(attr_delta)
    if moved and changed:
        klass = "both"
    elif moved:
        klass = "boundary_change"
    elif changed:
        klass = "rezone"
    else:
        return None                     # survived only sub-tolerance noise
    return {"class": klass, "symdiff_area": area, "hausdorff": hd, "attr_delta": attr_delta}

Run classify over the inner join of the two exports on parcel_id and you get the four real changes back — three rezone records and one boundary_change — while the 4,808 cosmetic differences return None and never reach a reviewer.

Verification and testing jump to heading

Confirm the pipeline suppressed noise rather than hid a real change:

• The known-unchanged parcel returns None. Feed the 0708-241-0913-4 pair from the diagnosis through classify and assert the result is None. If it still classifies as a change, your snap grid is finer than the reprojection residual — coarsen it.
• The known rezone is caught. Take a parcel the planning department confirmed was rezoned and assert class == "rezone" with the correct old→new zoning code in attr_delta. A miss means the changed column is not in WATCHED or the codes were normalized differently between exports.
• Change count matches ground truth. The total number of non-None classifications should equal the number of actions the county actually recorded that month. A count far above it means the noise floor is too low; far below means it is masking real moves.
• Round-trip stability. Diff an export against an untouched copy of itself. The changeset must be empty. Any output here is a bug in the canonicalization, not a real change.
• Threshold sensitivity sweep. Re-run the diff at a few AREA_FLOOR values around your chosen one and watch where the change count stabilizes; a plateau means you are safely above the noise band.

Failure recovery jump to heading

When a parcel cannot be diffed cleanly, the run must continue and stay reproducible:

• Quarantine, never guess. A parcel whose geometry is empty or invalid in either export after a buffer(0) repair is classified quarantine, written to a holding table with both raw geometries and the exception, and skipped. Comparing a repaired polygon against a raw one manufactures a fake symmetric-difference area, so it must not enter the changeset.
• Log the suppression count. Record how many pairs returned None after passing the hash check — that is your noise metric. A sudden spike usually means the upstream export changed its serialization, which is worth investigating before it starts hiding real changes behind a numbed-up floor.
• Fail loud on CRS or key problems. A missing CRS, a geographic working CRS, or duplicate parcel IDs raise immediately rather than producing a plausible-but-wrong changeset. A hard error at diff time is cheaper than a false alert that reaches underwriting.
• Persist the thresholds with the result. The same export pair diffed at a different SNAP_GRID or AREA_FLOOR yields a different changeset, so store the parameters alongside the output; without them the result is not reproducible and cannot be audited later against the parent Change Detection & Geometry Diffing run manifest.

Frequently asked questions jump to heading

Related jump to heading

• Parent topic: Change Detection & Geometry Diffing
• Section overview: Spatial Impact Analysis & Zoning Change Detection
• CRS Alignment Strategies — forcing both exports into one projected CRS before the diff

Up: Change Detection & Geometry Diffing · Spatial Impact Analysis & Zoning Change Detection