🧭 Concept Data cleaning Geometry & topology Updated

Topological Coverage: What "Clean" Means for Adjacent Polygons

Problem statement

Every polygon in the layer is valid. is_valid.all() returns True. The layer is still broken.

parcels.is_valid.all()          # True
parcels.geometry.is_empty.any() # False

# and yet
total = parcels.geometry.area.sum()
boundary_area = wards.geometry.area.sum()
total / boundary_area           # 1.0043  — 0.4% more land than exists

Validity is a property of one polygon: does its ring self-intersect, are its holes inside it. A set of polygons that are each individually valid can still overlap each other, leave slivers of unclaimed space between them, and share boundaries that do not quite match.

That set-level property is called a coverage, and it is what "clean" means for administrative boundaries, parcels, land use, soil types, catchments — any layer where the polygons are supposed to tile a region without gaps or overlaps.

There is no is_valid for it. You have to check it deliberately.

Quick answer

A polygon layer is a valid coverage when three things hold:

Rule Test Typical violation
No overlaps — interiors are disjoint sjoin(predicate="overlaps") returns nothing double-counted area
No gaps — the union has no interior holes unary_union has no interior rings sliver of unassigned land
Shared edges match exactly — adjacent polygons use identical vertices touches rather than overlaps between neighbours boundary drawn twice, slightly differently 
import geopandas as gpd
from shapely.ops import unary_union

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

    merged = unary_union(gdf.geometry)
    holes = sum(len(g.interiors) for g in getattr(merged, "geoms", [merged]))

    return {
        "features": len(gdf),
        "overlapping_pairs": len(pairs) // 2,      # each pair appears twice
        "interior_holes": holes,
        "sum_of_parts": round(gdf.geometry.area.sum(), 1),
        "area_of_union": round(merged.area, 1),
        "overlap_area": round(gdf.geometry.area.sum() - merged.area, 1),
    }

print(coverage_report(parcels))
# {'features': 12400, 'overlapping_pairs': 214, 'interior_holes': 87,
#  'sum_of_parts': 4834112.6, 'area_of_union': 4813880.2, 'overlap_area': 20232.4}

That last line is the number that matters: 20,232 m² of land is counted twice, because 214 pairs of parcels overlap. And 87 holes are land that belongs to nobody.

The two errors are independent, and a layer can have either without the other.

The three coverage rules

Three panels showing an overlap, a gap and a mismatched shared edge between adjacent polygons.
Each polygon is valid in all three. The set is not.

Step-by-step solution

Grid contrasting what validity checks against what coverage checks, per feature and per set.
is_valid answers a question about one shape. A coverage question needs the whole set.

Validity is per-feature; coverage is per-set 

# these all ask about ONE geometry
gdf.is_valid          # does this ring self-intersect?
gdf.is_empty          # does this shape enclose anything?
gdf.geom_type         # what kind of shape is this?

# these ask about the SET
gpd.sjoin(gdf, gdf, predicate="overlaps")   # do any two interiors intersect?
unary_union(gdf.geometry)                    # is the merged shape hole-free?
gdf.geometry.area.sum() - union.area         # how much is double-counted?

This distinction is why a cleaning script that runs make_valid on every row and stops feels complete but is not. make_valid fixes each polygon in isolation. It has no idea the polygon next door claims some of the same ground.

Overlaps: the double-counting error

An overlap means two polygons claim the same land. Every sum over the layer is then wrong by the overlapping area, and every point falling in the overlap joins to two parcels instead of one.

pairs = gpd.sjoin(parcels, parcels, predicate="overlaps")
pairs = pairs[pairs.index != pairs.index_right]

# how big are they? that tells you the cause
import pandas as pd
sizes = []
for left, right in zip(pairs.index, pairs["index_right"]):
    if left < right:                                  # each pair once
        a = parcels.geometry.iloc[left]
        b = parcels.geometry.iloc[right]
        sizes.append(a.intersection(b).area)

s = pd.Series(sizes)
print(s.describe())
# count    214
# 50%        0.004        ← sub-centimetre: digitising noise
# max     3211.700        ← square metres: a genuine data error

The size distribution is the diagnosis. Overlaps of a few square centimetres are floating-point artefacts from two people digitising the same line. Overlaps of thousands of square metres are two parcels genuinely claiming the same field — a records problem, not a geometry problem, and not something a script should silently resolve.

Gaps: the unassigned-land error 

merged = unary_union(parcels.geometry)
parts = getattr(merged, "geoms", [merged])
holes = [Polygon(ring) for part in parts for ring in part.interiors]

print(f"{len(holes)} holes, total {sum(h.area for h in holes):.1f} m²")
# 87 holes, total 412.6 m²

Same diagnostic logic. A hole of 0.02 m² between two parcels is a sliver — the two boundaries were drawn from slightly different vertices. A hole of 300 m² is a plot that nobody recorded, which is information rather than an error.

slivers = [h for h in holes if h.area < 1.0]
real_gaps = [h for h in holes if h.area >= 1.0]
print(f"{len(slivers)} slivers to close, {len(real_gaps)} genuine gaps to investigate")

Note that "holes in the union" only finds gaps enclosed by the coverage. A missing polygon at the edge produces no hole — the union is simply smaller. Comparing the union against a known outer boundary catches those:

missing = ward_boundary.difference(merged)
print(f"{missing.area:.1f} m² inside the ward and not covered")

Shared edges: the cause of both

The deepest version of the problem is that adjacent polygons do not share vertices. Two parcels that "obviously" abut may have boundaries digitised independently, so their edges run within a few millimetres of each other but never coincide.

a, b = parcels.geometry.iloc[0], parcels.geometry.iloc[1]

a.touches(b)      # False — they do not share a boundary exactly
a.overlaps(b)     # True  — one crosses the other by 3 mm
a.intersection(b).area   # 0.0031

touches flipping to overlaps is the signature. In a properly built coverage every neighbouring pair returns touches=True and overlaps=False — see spatial predicates explained.

This is why fixing overlaps and gaps one at a time is unsatisfying: they are symptoms of the same underlying condition, which is that the layer was never built as a coverage in the first place. The durable fix is to snap the shared boundaries onto common vertices — see how to find and fix gaps and overlaps in a polygon coverage.

What "planar" means, and when you need it

A planar coverage is the formal version: the polygons partition the plane, every point belongs to exactly one polygon, and shared edges are shared exactly. GIS tools that enforce it — PostGIS topology, ArcGIS topology rules, GRASS vector topology — store the edges once and derive the polygons from them, which makes gaps and overlaps structurally impossible.

Shapefile, GeoPackage and GeoJSON store each polygon independently, with its boundary repeated in both neighbours. That is simpler, portable, and offers no protection at all. Every coverage error in this article exists because the format allows it.

You need coverage to hold when:

  • you sum an attribute over the layer — overlaps inflate it, gaps deflate it
  • you join points to polygons — an overlap produces duplicate rows
  • you compute proportions — "40% of the district is residential" needs a denominator that is not double-counted
  • you dissolve or aggregate — slivers survive dissolve and multiply

You do not need it when the polygons are genuinely allowed to overlap: buffers, catchment areas, viewsheds, habitat ranges. Checking coverage on those is a false alarm.

Code examples

Example 1: a full coverage audit 

import geopandas as gpd
import pandas as pd
from shapely.geometry import Polygon
from shapely.ops import unary_union

def audit_coverage(gdf, *, sliver_area=1.0, boundary=None):
    """Everything wrong with a layer that is supposed to tile a region."""
    gdf = gdf.reset_index(drop=True)

    # overlaps, with sizes so you can tell noise from error
    pairs = gpd.sjoin(gdf, gdf, predicate="overlaps")
    pairs = pairs[pairs.index != pairs.index_right]
    seen, overlaps = set(), []
    for left, right in zip(pairs.index, pairs["index_right"]):
        key = tuple(sorted((left, right)))
        if key in seen:
            continue
        seen.add(key)
        area = gdf.geometry.iloc[left].intersection(gdf.geometry.iloc[right]).area
        overlaps.append({"a": int(left), "b": int(right), "area": area})

    # gaps enclosed by the coverage
    merged = unary_union(gdf.geometry)
    parts = getattr(merged, "geoms", [merged])
    holes = [Polygon(r) for p in parts for r in p.interiors]

    report = {
        "features": len(gdf),
        "overlaps": len(overlaps),
        "overlap_slivers": sum(1 for o in overlaps if o["area"] < sliver_area),
        "overlap_real": sum(1 for o in overlaps if o["area"] >= sliver_area),
        "overlap_area": round(sum(o["area"] for o in overlaps), 3),
        "holes": len(holes),
        "hole_slivers": sum(1 for h in holes if h.area < sliver_area),
        "hole_real": sum(1 for h in holes if h.area >= sliver_area),
        "hole_area": round(sum(h.area for h in holes), 3),
        "sum_of_parts": round(gdf.geometry.area.sum(), 1),
        "union_area": round(merged.area, 1),
    }
    if boundary is not None:
        report["uncovered_at_edge"] = round(boundary.difference(merged).area, 1)
    return report, pd.DataFrame(overlaps), holes

The split between _slivers and _real is the point of the whole function. A layer with 200 slivers and 0 real errors is mechanically fixable; a layer with 3 real overlaps needs somebody to decide who owns the land.

Example 2: the invariant worth asserting 

def assert_coverage(gdf, *, tolerance=0.001):
    """Sum of parts must equal the union, within a tolerance."""
    parts = gdf.geometry.area.sum()
    union = unary_union(gdf.geometry).area
    drift = abs(parts - union) / union

    assert drift < tolerance, (
        f"coverage broken: parts sum to {parts:,.1f} but union is {union:,.1f} "
        f"({drift:.3%} double-counted)"
    )

One assertion, and it catches every overlap in the layer without enumerating them. It belongs in the validation step of any pipeline whose output is a coverage — and it is far cheaper than the pairwise join, so it is the right check to run every time.

Example 3: visualising where the problems are 

import matplotlib.pyplot as plt

report, overlaps, holes = audit_coverage(parcels)

fig, ax = plt.subplots(figsize=(11, 11))
parcels.plot(ax=ax, facecolor="none", edgecolor="#cbd5e1", linewidth=0.4)

if holes:
    gpd.GeoSeries(holes, crs=parcels.crs).plot(
        ax=ax, facecolor="#f59e0b", edgecolor="#b45309", label="gaps")

if len(overlaps):
    shapes = [
        parcels.geometry.iloc[r.a].intersection(parcels.geometry.iloc[r.b])
        for r in overlaps.itertuples()
    ]
    gpd.GeoSeries(shapes, crs=parcels.crs).plot(
        ax=ax, facecolor="#ef4444", edgecolor="#991b1b", label="overlaps")

ax.set_title(f"{report['overlaps']} overlaps · {report['holes']} gaps")
ax.set_axis_off()

Coverage errors cluster. A map of them almost always shows a single boundary line where two datasets were joined, or one surveyor's block — which turns "fix 214 slivers" into "re-request one file".

Explanation

Triage rows separating sliver-sized coverage errors from genuine data errors and what each needs.
The size of the error decides whether a script may fix it.

The reason coverage errors are so common comes down to how the data is stored. In a shapefile or GeoPackage, each polygon carries its own complete boundary. Two adjacent parcels store the line between them twice, once in each feature, as two independent lists of coordinates.

Nothing keeps those two lists identical. They start identical if the layer was built by splitting one polygon; they diverge the moment anyone edits one side, imports from a different source, reprojects at a different precision, or simplifies. And a divergence of one micrometre is enough to turn touches into overlaps.

Topological data models solve this by storing the shared edge once — the parcels are defined by reference to the edge rather than by copying it. Move the edge and both parcels move together, so overlaps and gaps cannot arise. GRASS, PostGIS topology and ArcGIS geodatabase topology all work this way. The cost is complexity and portability, which is why the simple formats won for interchange and why this article exists.

The practical consequence is the distinction that runs through everything above: coverage errors have a size, and the size tells you who should fix them. Sub-millimetre errors are artefacts of the storage model and can be closed mechanically — snapping to a grid, or a tiny buffer-and-dissolve round trip. Errors measured in metres are disagreements about the world, and resolving them by picking whichever polygon happens to come first in the file is worse than leaving them, because it destroys the evidence that a disagreement existed.

One last asymmetry worth knowing: overlaps are always errors in a coverage; gaps sometimes are not. A hole where a river runs, a road reserve, or an unregistered plot is real. That is why the audit reports gaps and slivers separately rather than presenting a single "problems" count.

Edge cases or notes

  • Buffers, catchments and ranges are not coverages. Overlap is their normal state; do not audit them for it.
  • unary_union on a large layer is expensive — minutes for hundreds of thousands of polygons. The area-sum assertion is the cheap proxy.
  • Overlaps at exactly zero area are shapes that touch along an edge; overlaps should return False for them. If it returns True with area 0, the geometry is probably invalid.
  • Reprojecting can create coverage errors where none existed, because coordinates are recomputed independently per feature.
  • simplify() reliably breaks a coverage. Each boundary is simplified separately, so shared edges diverge. Use topology-preserving simplification (topojson, or PostGIS ST_SimplifyPreserveTopology across the set).
  • MultiPolygons complicate hole counting — iterate .geoms before reading .interiors.
  • A layer can be a valid coverage of the wrong region — covering land that is not in scope. Compare the union against an authoritative boundary.
  • Dissolve does not repair coverage. It merges by attribute and happily carries slivers through into the result.

FAQ

Why does is_valid return True for a layer with overlaps?

is_valid is a per-geometry check — it asks whether one polygon breaks the OGC rules. Overlaps are a relationship between two polygons, so no per-geometry check can see them.

How do I know if my layer should be a coverage?

Ask whether two features are allowed to claim the same ground. Parcels, wards and land use: no, so it is a coverage. Buffers, catchments and habitat ranges: yes, so it is not.

What is the fastest coverage check?

Compare gdf.geometry.area.sum() against unary_union(gdf.geometry).area. Any difference is double-counted area. It will not tell you where, but it answers "is there a problem" in one line.

Are gaps always errors?

No. A hole where a river, road reserve or unregistered plot sits is real. Only sliver-sized gaps — well under a square metre — are reliably artefacts.

Does make_valid fix coverage errors?

No. It repairs each polygon in isolation and knows nothing about its neighbours. A layer can be 100% valid and still overlap itself everywhere.

Why did reprojecting break my coverage?

Each feature's coordinates are transformed independently, so a shared edge stored twice can round to slightly different values. Reproject the coverage once, early, and check afterwards.

Can I prevent this rather than repair it?

Yes — build the layer in a topological model (PostGIS topology, GRASS) where shared edges are stored once, or derive all polygons from a single set of lines by polygonising.