Geospatial Schema Standardization & Mapping Automation

Detecting Gaps and Overlaps in Parcel Coverages

A parcel fabric claims to partition the ground: every square metre inside the administrative boundary belongs to exactly one parcel. Two defects break that claim in opposite directions — a gap is land nobody owns, an overlap is land two people own — and both are invisible on a map at any sensible zoom because the offenders are typically a few centimetres wide. This procedure finds all of them, classifies the ones that are digitizing noise separately from the ones that are conflicting claims, and reports each exactly once even when it straddles a tile edge. It implements the detection half of topology rule enforcement, the gate inside Spatial Data Quality Validation & Geometry Integrity that owns relationships between features.

The steps map to configure (Steps 1–2, boundary and tiling), execute (Steps 3–4, the two detections), validate (Step 5, classification and dedup), and log (Step 5’s report rows). Everything here assumes the geometries have already cleared the validity gate — a union over invalid input produces a result no one can defend.

Prerequisites checklist

Step 1: Establish the outer boundary before anything else

Gap detection is subtraction, and subtraction needs a minuend. Without an explicit boundary layer, the “gap” outside the outermost parcel is unbounded, and the procedure either crashes or reports the rest of the plane as missing.

python
# geopandas >= 0.14 — Python 3.10+
import geopandas as gpd

CRS = "EPSG:25832"

boundary = gpd.read_file("reference/municipality_boundary.gpkg").to_crs(CRS)
assert len(boundary) == 1, "the coverage boundary must be a single authoritative polygon"
assert boundary.geometry.iloc[0].is_valid, "boundary must be valid before it is subtracted from"

parcels = gpd.read_file("input/parcels.gpkg").to_crs(CRS)
assert parcels.crs.is_projected, "tolerances are in linear units; a geographic CRS is meaningless here"

Use the boundary the authority publishes, not the dissolved parcel layer. Dissolving the parcels to make a boundary guarantees zero gaps at the edge, which is the one place gaps most often are — a strip of unregistered land along a river or a road reserve.

Step 2: Tile the extent with a halo

python
# shapely >= 2.0, geopandas >= 0.14 — Python 3.10+
from shapely.geometry import box

TILE = 2000.0      # CRS units
HALO = 25.0        # must exceed the largest tolerance in the manifest


def tiles(bounds, size: float):
    minx, miny, maxx, maxy = bounds
    y = miny
    while y < maxy:
        x = minx
        while x < maxx:
            yield box(x, y, min(x + size, maxx), min(y + size, maxy))
            x += size
        y += size


def tile_slice(gdf: gpd.GeoDataFrame, tile) -> gpd.GeoDataFrame:
    """Features intersecting the tile plus its halo; the halo prevents false edge gaps."""
    return gdf[gdf.intersects(tile.buffer(HALO))]

The halo is a correctness device, not a performance one. Without it, every tile boundary cuts through parcels and the subtraction reports a false gap along each cut. With a halo narrower than the largest tolerance, a genuine breach straddling the edge can be measured differently in the two tiles that see it, and the classification becomes unstable.

Step 3: Find gaps by union and subtract

python
# shapely >= 2.0 — Python 3.10+
from shapely import get_parts, set_precision, union_all

GRID = 0.001       # 1 mm — kills arithmetic-noise slivers before they are counted


def gaps_in_tile(parcel_geoms, tile, boundary_geom):
    """Every piece of the boundary inside this tile that no parcel covers."""
    area_of_interest = boundary_geom.intersection(tile)
    if area_of_interest.is_empty:
        return []

    covered = union_all([set_precision(g, GRID) for g in parcel_geoms])
    remainder = area_of_interest.difference(covered)
    if remainder.is_empty:
        return []

    return [part for part in get_parts(remainder) if part.area > 0.0]
Union, Subtract, Explode — How a Gap Is Found Three panels run left to right. The first shows the authoritative boundary clipped to the tile, drawn as one outline. The second shows the parcels within that tile unioned into a single covered area, which fills most but not all of the boundary. The third shows the difference between the two: two small remainder polygons, one a long thin strip along a shared edge and one a compact square, which are exploded into individual gap features and then classified by area and thinness. A note records that snapping every geometry to a one millimetre grid before the union is what stops the shared boundaries from producing thousands of hair-thin false gaps. boundary ∩ tile the area required to be covered union of the parcels snapped to a 1 mm grid first difference → gaps exploded into parts, then classified compact = gap · hair-thin = sliver minus equals

Two details carry the weight. set_precision before the union removes the sub-millimetre discrepancies that would otherwise become thousands of hair-thin “gaps” along every shared boundary — the noise floor of any coverage assembled from more than one source. And clipping the boundary to the tile before subtracting keeps the operation bounded: unioning a whole county’s parcels in one call is where this procedure usually runs out of memory.

Step 4: Find overlaps with an indexed pair query

python
# shapely >= 2.0 — Python 3.10+
from shapely import STRtree


def overlaps_in_tile(gdf: gpd.GeoDataFrame):
    geoms = gdf.geometry.values
    ids = gdf["feature_id"].tolist()
    tree = STRtree(geoms)
    left_idx, right_idx = tree.query(geoms, predicate="intersects")

    seen: set[tuple[str, str]] = set()
    for i, j in zip(left_idx.tolist(), right_idx.tolist()):
        if i == j:
            continue
        key = tuple(sorted((ids[i], ids[j])))
        if key in seen:
            continue
        seen.add(key)

        shared = geoms[i].intersection(geoms[j])
        if shared.is_empty or shared.area <= 0.0:
            continue                     # a shared edge is contact, not overlap
        yield key, shared

The shared.area <= 0.0 test is what separates this from a naive intersects report. In a healthy coverage, every pair of neighbouring parcels intersects — along their shared boundary, in a line of zero area. Reporting those as overlaps produces a breach list the length of the layer and teaches everyone to ignore it.

Step 5: Classify, deduplicate across tiles, and report once

python
# shapely >= 2.0, pyarrow >= 14 — Python 3.10+
import hashlib
import math

AREA_TOL = 0.05          # m2
SHAPE_INDEX_MAX = 0.02


def shape_index(geom) -> float:
    p = geom.length
    return 0.0 if p <= 0 else (4.0 * math.pi * geom.area) / (p * p)


def classify(geom, kind: str) -> str:
    thin = shape_index(geom) < SHAPE_INDEX_MAX
    small = geom.area <= AREA_TOL
    if thin and small:
        return "sliver"
    return kind          # "gap" or "overlap"


def breach_key(geom) -> str:
    """Stable hash of the breach geometry so the same one seen in two tiles collapses."""
    return hashlib.sha256(geom.wkb).hexdigest()[:32]


emitted: set[str] = set()

def emit(geom, kind: str, left_id: str | None, right_id: str | None, rows: list):
    key = breach_key(set_precision(geom, GRID))
    if key in emitted:
        return                                  # already reported from a neighbouring tile
    emitted.add(key)
    rows.append({
        "breach_id":      key,
        "kind":           kind,
        "classification": classify(geom, kind),
        "area":           geom.area,
        "shape_index":    shape_index(geom),
        "left_id":        left_id,
        "right_id":       right_id,
        "wkb":            geom.wkb,             # the evidence a steward opens
    })

Hashing the precision-snapped WKB is what makes exactly-once reporting work without coordinating between tiles: the same gap polygon computed in two neighbouring tiles produces byte-identical WKB after snapping, so the second sighting collapses. Hashing the unsnapped geometry does not work — the two computations differ in the last bits and produce two hashes, which is exactly the double-count the halo was supposed to prevent.

Finally, the dataset-level verdict:

python
total_gap_area = sum(r["area"] for r in rows if r["kind"] == "gap")
coverage_area = boundary.geometry.iloc[0].area
closure_ratio = total_gap_area / coverage_area

if closure_ratio > 1e-6:                 # published budget: 1 part per million
    raise CoverageNotClosed(
        f"gap area {total_gap_area:.3f} m² = {closure_ratio:.2e} of coverage — above budget"
    )

A coverage that fails closure is not publishable as a partition, however few individual gaps it has. That is a different judgement from the per-breach one: a hundred 1 cm² slivers is a healthy fabric with digitizing noise; one 400 m² gap is a missing parcel.

Verification

python
# pytest >= 7, geopandas >= 0.14
def test_known_gap_is_found_exactly_once(fixture_coverage):
    rows = run_coverage_check(fixture_coverage)
    gaps = [r for r in rows if r["kind"] == "gap"]
    assert len(gaps) == 1
    assert abs(gaps[0]["area"] - 12.5) < 1e-6        # the fixture's deliberate gap

def test_edge_straddling_gap_is_not_double_counted(fixture_coverage):
    # The fixture's gap sits exactly on a tile boundary at x = 2000.
    rows = run_coverage_check(fixture_coverage, tile_size=2000.0)
    assert len({r["breach_id"] for r in rows}) == len(rows)

def test_shared_boundaries_are_not_reported_as_overlaps(fixture_coverage):
    assert [r for r in rows if r["kind"] == "overlap" and r["area"] == 0.0] == []

A healthy county run logs something close to this — the ratio of slivers to genuine breaches is the number to watch over time:

text
INFO quality.topology tiles=196 parcels=214,338 pairs=612,004
INFO quality.topology gaps=41 (slivers 38) overlaps=17 (slivers 14)
INFO quality.topology closure_ratio=3.10e-08 budget=1.00e-06 result=pass

Troubleshooting

Symptom Likely cause Fix
Tens of thousands of hair-thin gaps along every boundary No precision snapping before the union Apply set_precision(g, 0.001) to each geometry first; the noise floor disappears
A long gap reported twice with slightly different areas Halo narrower than the breach, so each tile saw a fragment Widen tiling.overlap beyond the largest expected breach, then re-run
Memory exhausted during union_all The union is being computed over the whole layer rather than per tile Clip the boundary to the tile and union only that tile’s parcels plus halo
Every parcel pair reported as an overlap The area test is missing, so shared edges count Keep only intersections with area > 0.0
Closure fails but no individual gap looks wrong Many small gaps summing above the budget, or a boundary layer from a different vintage than the parcels Compare the boundary’s date with the parcel delivery; a stale boundary shifts the whole edge
Overlaps appear only along one municipal edge Two sources digitized the same boundary in different projections Check the harmonization step — see merging UTM zone boundary datasets without slivers