How to Find and Fix Gaps and Overlaps in a Polygon Coverage

Problem statement

The parcels layer should tile the district exactly. It does not.

parcels.geometry.area.sum()              # 4_834_112.6
unary_union(parcels.geometry).area       # 4_813_880.2

Twenty thousand square metres are counted twice. Somewhere in 12,400 polygons there are pairs that overlap, and holes where nothing is claimed — and is_valid reports everything is fine, because validity is a per-polygon property and this is a property of the set.

Fixing it by hand is not an option at this scale. Fixing it blindly is worse: a script that snaps everything to a 1 m grid will close the 214 slivers and also move real boundaries by up to half a metre.

The workable approach is to separate the artefacts from the disagreements, fix the artefacts mechanically, and report the rest.

Quick answer

Measure, split by size, snap the small ones, report the big ones:

import geopandas as gpd
from shapely.ops import unary_union
from shapely import set_precision

SLIVER = 1.0        # m² — below this it is a digitising artefact
GRID    = 0.001     # m  — snap coordinates to 1 mm

def fix_coverage(gdf, *, sliver_area=SLIVER, grid=GRID):
    report = {"in": len(gdf)}

    # 1. every geometry valid, or the rest of the analysis is unreliable
    gdf = gdf.assign(geometry=gdf.geometry.make_valid())

    # 2. snap coordinates to a common grid — this is what closes the slivers
    gdf = gdf.assign(geometry=set_precision(gdf.geometry.values, grid))
    gdf = gdf[~gdf.geometry.is_empty & gdf.geometry.notna()]

    # 3. measure what is left
    parts = gdf.geometry.area.sum()
    union = unary_union(gdf.geometry)
    report["overlap_area"] = round(parts - union.area, 3)
    report["holes"] = [
        Polygon(r) for p in getattr(union, "geoms", [union]) for r in p.interiors
    ]
    report["sliver_holes"] = sum(1 for h in report["holes"] if h.area < sliver_area)
    report["real_holes"] = sum(1 for h in report["holes"] if h.area >= sliver_area)
    return gdf, report

before:  214 overlapping pairs, 20,232.4 m² double-counted, 87 holes
after :    3 overlapping pairs,  3,211.7 m² double-counted,  4 holes

211 of 214 overlaps were sub-millimetre and closed by snapping. The three that remain are between 400 and 3,200 m² — two parcels genuinely claiming the same field, which is a records question, not a geometry one.

Symptom	Size	What it is	What to do
overlap	< 1 m²	two people digitised the same line	snap to a grid
overlap	> 1 m²	a real dispute about ownership	report, do not resolve
hole	< 1 m²	sliver between mismatched edges	snap, then fill remainder
hole	> 1 m²	genuinely unassigned land	investigate
missing at the edge	any	a polygon that was never supplied	compare against a boundary

Artefacts versus disagreements

Triage rows separating sliver-sized coverage errors from real ones and the action each takes. — One threshold decides whether a script may act. Everything above it is a person's decision.

Step-by-step solution

Vertical steps from validity through snapping, measuring, filling slivers and reporting the rest. — Snap before you measure. Measuring first counts artefacts you are about to remove.

1. Make everything valid first

Every measurement below is unreliable on invalid geometry, because a self-intersecting polygon has an ambiguous interior — so overlaps and area both lie.

invalid = (~gdf.is_valid).sum()
if invalid:
    print(f"repairing {invalid} invalid geometries first")
    gdf = gdf.assign(geometry=gdf.geometry.make_valid())

# make_valid can change the type — keep only polygonal parts
gdf = gdf[gdf.geom_type.isin(["Polygon", "MultiPolygon"])]
assert gdf.is_valid.all()

2. Snap to a precision grid — the single most effective fix

shapely.set_precision rounds every coordinate to a grid and, crucially, does it consistently across the whole set. Two boundaries that differed by 10⁻⁹ now round to the same number, so the edge becomes genuinely shared.

from shapely import set_precision

gdf = gdf.assign(geometry=set_precision(gdf.geometry.values, 0.001))   # 1 mm

Choosing the grid size is the one judgement call:

# what precision does the data actually carry?
import numpy as np
coords = np.concatenate([np.asarray(g.exterior.coords) for g in gdf.geometry
                         if g.geom_type == "Polygon"])
decimals = [len(str(c).split(".")[-1]) for c in coords[:1000, 0]]
print(f"typical decimal places: {np.median(decimals):.0f}")

1 mm (0.001) — safe for survey-grade metric data; closes floating-point noise and nothing else.
1 cm (0.01) — closes noise from reprojection round trips.
1 m (1.0) — too coarse for parcels; will move real boundaries. Only for coarse regional data.

Snapping is not free: it moves every vertex by up to half the grid. At 1 mm on a parcel layer that is invisible; at 1 m it is vandalism. Never pick a grid larger than the smallest feature you care about.

# check the damage before accepting it
moved = (gdf.geometry.area - original.geometry.area).abs()
print(f"largest area change from snapping: {moved.max():.4f} m²")

3. Measure what snapping did not fix

def coverage_errors(gdf, sliver_area=1.0):
    pairs = gpd.sjoin(gdf, gdf, predicate="overlaps")
    pairs = pairs[pairs.index != pairs.index_right]

    seen, overlaps = set(), []
    for a, b in zip(pairs.index, pairs["index_right"]):
        key = tuple(sorted((a, b)))
        if key in seen:
            continue
        seen.add(key)
        shape = gdf.geometry.iloc[a].intersection(gdf.geometry.iloc[b])
        overlaps.append({"a": int(a), "b": int(b), "area": shape.area, "geometry": shape})

    union = unary_union(gdf.geometry)
    holes = [Polygon(r) for p in getattr(union, "geoms", [union]) for r in p.interiors]

    return (
        gpd.GeoDataFrame(overlaps, crs=gdf.crs) if overlaps else gpd.GeoDataFrame(),
        gpd.GeoDataFrame(geometry=holes, crs=gdf.crs) if holes else gpd.GeoDataFrame(),
    )

overlaps, holes = coverage_errors(gdf)
print(f"{len(overlaps)} overlaps ({(overlaps.area > 1).sum()} real), "
      f"{len(holes)} holes ({(holes.area > 1).sum()} real)")

4. Fill sliver holes by assigning them to a neighbour

A sliver between two parcels belongs to one of them. Which one is arbitrary — so pick a rule and apply it consistently. The usual rule: give it to the neighbour it shares the longest boundary with.

def fill_slivers(gdf, holes, max_area=1.0):
    """Assign each sliver hole to its longest-shared-edge neighbour."""
    slivers = holes[holes.geometry.area < max_area]
    filled = 0

    for sliver in slivers.geometry:
        touching = gdf[gdf.geometry.intersects(sliver)]
        if touching.empty:
            continue
        # longest shared boundary wins; ties break on the lower index
        shares = touching.geometry.apply(
            lambda g: g.intersection(sliver.buffer(1e-9)).length
        )
        winner = shares.idxmax()
        gdf.loc[winner, "geometry"] = gdf.loc[winner, "geometry"].union(sliver)
        filled += 1

    return gdf, filled

Two properties make this defensible: the rule is deterministic (same input, same assignment every run), and it conserves area — every square metre of sliver ends up in exactly one parcel rather than being deleted.

5. Resolve small overlaps by clipping the later feature

For sub-threshold overlaps, one polygon must give up the disputed strip. Again the rule matters more than which side wins:

def resolve_small_overlaps(gdf, overlaps, max_area=1.0):
    """Clip the overlap out of the higher-index feature. Deterministic."""
    resolved = 0
    for row in overlaps[overlaps.geometry.area < max_area].itertuples():
        loser = max(row.a, row.b)                    # deterministic tie-break
        gdf.loc[loser, "geometry"] = (
            gdf.loc[loser, "geometry"].difference(row.geometry)
        )
        resolved += 1
    return gdf, resolved

Using max(a, b) rather than "whichever comes first in the file" matters: file order changes between runs and between machines, and a non-deterministic cleaning step makes every downstream total unreproducible. See idempotency explained.

6. Report the rest — do not resolve it

real_overlaps = overlaps[overlaps.geometry.area >= 1.0]
real_holes = holes[holes.geometry.area >= 1.0]

if len(real_overlaps):
    real_overlaps.to_file(out / "review_overlaps.gpkg", driver="GPKG")
    print(f"⚠ {len(real_overlaps)} overlaps over 1 m² need a human decision:")
    for r in real_overlaps.itertuples():
        print(f"    parcels {r.a} and {r.b} share {r.area:,.1f} m²")

Writing them to a file that somebody opens in QGIS is the deliverable. A script that silently picks a winner for a 3,000 m² dispute has destroyed the only evidence that a dispute existed.

Code examples

Example 1: the complete repair, with a report

import geopandas as gpd
from shapely import set_precision
from shapely.geometry import Polygon
from shapely.ops import unary_union

def repair_coverage(gdf, *, grid=0.001, sliver=1.0):
    original_area = gdf.geometry.area.sum()
    report = {"features_in": len(gdf)}

    # validity
    bad = int((~gdf.is_valid).sum())
    if bad:
        gdf = gdf.assign(geometry=gdf.geometry.make_valid())
        report["repaired_invalid"] = bad
    gdf = gdf[gdf.geom_type.isin(["Polygon", "MultiPolygon"])].reset_index(drop=True)

    # snap
    gdf = gdf.assign(geometry=set_precision(gdf.geometry.values, grid))
    gdf = gdf[~gdf.geometry.is_empty].reset_index(drop=True)
    report["snapped_to_grid_m"] = grid

    # measure
    overlaps, holes = coverage_errors(gdf, sliver)
    report["overlaps_found"] = len(overlaps)
    report["holes_found"] = len(holes)

    # mechanical fixes
    if len(holes):
        gdf, n = fill_slivers(gdf, holes, sliver)
        report["slivers_filled"] = n
    if len(overlaps):
        gdf, n = resolve_small_overlaps(gdf, overlaps, sliver)
        report["small_overlaps_clipped"] = n

    # what is left for a person
    overlaps_after, holes_after = coverage_errors(gdf, sliver)
    report["overlaps_remaining"] = len(overlaps_after)
    report["holes_remaining"] = len(holes_after)
    report["area_change_m2"] = round(gdf.geometry.area.sum() - original_area, 3)
    report["features_out"] = len(gdf)

    return gdf, report, overlaps_after, holes_after

clean, report, review_overlaps, review_holes = repair_coverage(parcels)
print(report)
# {'features_in': 12400, 'repaired_invalid': 31, 'snapped_to_grid_m': 0.001,
#  'overlaps_found': 214, 'holes_found': 87, 'slivers_filled': 83,
#  'small_overlaps_clipped': 211, 'overlaps_remaining': 3, 'holes_remaining': 4,
#  'area_change_m2': -18.442, 'features_out': 12400}

features_in == features_out and a small area_change_m2 are the two numbers to check. If the feature count changed, something dropped a polygon; if the area moved by thousands, the grid is too coarse.

Example 2: rebuilding a coverage from its boundaries

When a layer is badly broken — hundreds of real overlaps, systematic mismatches — repairing feature by feature is hopeless. Rebuilding from the lines is often faster and always gives a true coverage:

from shapely.ops import unary_union, polygonize

def rebuild_from_edges(gdf, grid=0.001):
    """Explode to boundaries, node them, re-polygonise. Guarantees a coverage."""
    snapped = set_precision(gdf.geometry.values, grid)
    edges = unary_union([g.boundary for g in snapped])      # nodes every crossing
    faces = list(polygonize(edges))                          # gap- and overlap-free by construction

    out = gpd.GeoDataFrame(geometry=faces, crs=gdf.crs)
    # re-attach attributes by point-in-polygon against the original
    reps = out.copy()
    reps["geometry"] = out.geometry.representative_point()
    out = out.join(
        gpd.sjoin(reps, gdf, predicate="within", how="left")
           .drop(columns=["geometry", "index_right"])
    )
    return out

polygonize builds faces from a noded line network, so gaps and overlaps are structurally impossible in the output. The cost is that attributes have to be re-attached, and any face that falls in a former overlap gets its attributes from whichever original polygon its representative point lands in. Use it when the alternative is worse.

Example 3: guarding against regression

def test_output_is_a_coverage(cleaned_parcels):
    parts = cleaned_parcels.geometry.area.sum()
    union = unary_union(cleaned_parcels.geometry).area
    assert abs(parts - union) / union < 1e-6, "output overlaps itself"

    holes = [r for p in getattr(unary_union(cleaned_parcels.geometry), "geoms", [])
             for r in p.interiors]
    big = [h for h in holes if Polygon(h).area > 1.0]
    assert not big, f"{len(big)} unexplained gaps over 1 m²"

Cheap, and it catches the regression that matters: a later pipeline step — a reprojection, a simplify, a dissolve — quietly breaking the coverage you just fixed.

Explanation

Two panels showing two near-identical boundaries before and after snapping to a common grid. — Snapping does not move the boundary. It makes two copies of it agree.

The reason snapping works so well is that it attacks the cause rather than the symptom.

In a shapefile or GeoPackage, the boundary between two parcels is stored twice — once in each feature — as two independent coordinate lists. Those lists start identical if the layer was built by splitting, and diverge the moment anything touches one side: an edit, a reprojection, a simplify, an import through a format with different precision. A divergence of 10⁻⁹ metres is enough to turn touches into overlaps and produce a sliver you can measure but not see.

set_precision rounds every coordinate in the layer to the same grid. Two vertices that differed in the twelfth decimal place now round to the identical value, so the two copies of the boundary become byte-identical and the shapes genuinely share the edge. The overlap is not clipped away — it ceases to exist, because the two polygons now agree where the line is.

That is why the order in step 3 matters. Measuring before snapping counts 214 overlaps, 211 of which are about to evaporate; the interesting number is what survives.

The threshold is the other half of the method, and it encodes a judgement worth stating explicitly: a script may fix things that are wrong because of how data is stored, and may not fix things that are wrong because people disagree. A 0.003 m² overlap is a storage artefact — no surveyor intended it, and no information is lost by removing it. A 3,000 m² overlap is two records asserting incompatible things about the world, and any automatic resolution silently picks a winner in a dispute the script cannot understand.

Rebuilding from edges (example 2) is the structural version of the same insight. Instead of repairing polygons that each carry their own copy of a shared line, it throws away the polygons, keeps the lines, nodes them so every crossing becomes a vertex, and derives the faces. The output cannot have gaps or overlaps, because the faces are defined by the edges rather than duplicating them — which is how topological GIS models work natively.

Edge cases or notes

set_precision needs Shapely 2.0+. On 1.x, the nearest equivalents are shapely.wkt.dumps(g, rounding_precision=…) or a manual coordinate round trip.
set_precision can produce empty or invalid geometry for features smaller than the grid. Filter afterwards and count what went.
Snapping changes area slightly, always. Assert the change is small rather than zero.
polygonize drops attributes entirely — they must be re-attached by point-in-polygon, and that is lossy where the original overlapped.
representative_point() is guaranteed inside the polygon; centroid is not, for concave shapes. Use the former for attribute transfer.
Do not snap before reprojecting. Reproject to the working CRS first, then snap once in that CRS.
A layer with holes at the outer edge is not detected by interiors — the union is simply smaller. Compare against an authoritative boundary.
Very large layers make the pairwise overlap join expensive. Use the area-sum check to decide whether it is worth running at all.
QGIS "Fix geometries" and v.clean in GRASS do a comparable job interactively, and are worth reaching for on a one-off.

Internal links

Topological coverage: what "clean" means for adjacent polygons — the concept and the audit
How to snap and align geometries to fix slivers and gaps — snapping in more depth
Topology in GIS explained: slivers, gaps and shared boundaries — where slivers come from
How to fix invalid geometries in Python (GeoPandas) — step 1
My cleaned layer has more rows than the original — what make_valid can do to a count
Coordinate precision and floating point in GIS explained — why the twelfth decimal place matters
Repair, reject or flag? Choosing what cleaning should do — why the big overlaps are reported, not fixed
Idempotency explained: why a GIS job must be safe to re-run — why the tie-breaks are deterministic

FAQ

What grid size should I snap to?

The smallest that closes the artefacts: 1 mm for survey-grade metric data, 1 cm if the data has been through a reprojection round trip. Never larger than the smallest feature you care about.

Will snapping move my boundaries?

By up to half the grid size, yes. At 1 mm that is invisible; at 1 m it is destructive. Measure the area change and assert it is small.

Should I fix large overlaps automatically?

No. An overlap of thousands of square metres is two records disagreeing about the world. Report it, write it to a review file, and let somebody decide.

Which parcel should get a sliver?

Whichever shares the longest boundary with it — that is the conventional rule, it is deterministic, and it conserves area. What matters most is that the rule is fixed, not which rule you pick.

Why measure after snapping rather than before?

Because most of the errors disappear when you snap. Measuring first gives you a count dominated by artefacts you are about to remove.

When should I rebuild from edges instead of repairing?

When there are hundreds of real overlaps or the mismatches are systematic. polygonize guarantees a coverage; the cost is re-attaching attributes.

Does dissolve fix coverage problems?

No. dissolve merges by attribute and carries slivers straight through — it often makes them more visible rather than removing them.

How to Find and Fix Gaps and Overlaps in a Polygon Coverage #

Problem statement #

Quick answer #

Artefacts versus disagreements #

Step-by-step solution #

1. Make everything valid first #

2. Snap to a precision grid — the single most effective fix #

3. Measure what snapping did not fix #

4. Fill sliver holes by assigning them to a neighbour #

5. Resolve small overlaps by clipping the later feature #

6. Report the rest — do not resolve it #

Code examples #

Example 1: the complete repair, with a report #

Example 2: rebuilding a coverage from its boundaries #

Example 3: guarding against regression #

Explanation #

Edge cases or notes #

Internal links #

FAQ #

What grid size should I snap to? #

Will snapping move my boundaries? #

Should I fix large overlaps automatically? #

Which parcel should get a sliver? #

Why measure after snapping rather than before? #

When should I rebuild from edges instead of repairing? #

Does dissolve fix coverage problems? #