Topological Coverage: What "Clean" Means for Adjacent Polygons

On this page

Problem statement

Every polygon in the layer is valid. is_valid.all() returns True. The layer is still broken. Here is HM Land Registry's layer of 21,253 freehold parcels for Rutland:

print(parcels.is_valid.all(), parcels.geometry.is_empty.any())

total = parcels.geometry.area.sum()
union = parcels.geometry.union_all().area
print(f"{total / 1e6:.1f} km² summed, {union / 1e6:.1f} km² covered, ratio {total / union:.4f}")

True False
429.0 km² summed, 398.0 km² covered, ratio 1.0779

The parcels add up to 31 km² more land than they cover — 7.8% more. Somewhere, the same ground is claimed by more than one polygon.

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.

GeoPandas has no is_valid for it. Shapely 2.1 has one on GEOS 3.12 or later: shapely.coverage_is_valid() flags overlaps and shared edges that do not match exactly, shapely.coverage_invalid_edges() says where, and a gap_width argument makes both flag gaps narrower than that width too. Most code in the wild predates it, so it helps to know what it checks.

Quick answer

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

Rule	Test	Typical violation
No overlaps — interiors are disjoint	`sjoin` with `overlaps` and `contains` returns nothing	double-counted area
No gaps — the union has no interior holes	the 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
import shapely

def coverage_report(gdf):
    frame = gpd.GeoDataFrame(geometry=gdf.geometry.reset_index(drop=True))
    overlapping = gpd.sjoin(frame, frame, predicate="overlaps")
    nested = gpd.sjoin(frame, frame, predicate="contains")
    nested = nested[nested.index != nested["index_right"]]

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

    return {
        "features": len(frame),
        "coverage_is_valid": bool(shapely.coverage_is_valid(frame.geometry.values)),
        "overlapping_pairs": len(overlapping) // 2,      # each pair appears twice
        "nested_pairs": len(nested),                     # one polygon wholly inside another
        "interior_holes": holes,
        "sum_of_parts_km2": round(float(frame.geometry.area.sum()) / 1e6, 2),
        "area_of_union_km2": round(merged.area / 1e6, 2),
    }

print(coverage_report(parcels))
print(coverage_report(wards))

{'features': 21253, 'coverage_is_valid': False, 'overlapping_pairs': 11453, 'nested_pairs': 4106, 'interior_holes': 9256, 'sum_of_parts_km2': 428.99, 'area_of_union_km2': 397.99}
{'features': 15, 'coverage_is_valid': True, 'overlapping_pairs': 0, 'nested_pairs': 0, 'interior_holes': 0, 'sum_of_parts_km2': 393.73, 'area_of_union_km2': 393.73}

The parcels fail every rule: 11,453 pairs overlap, 4,106 parcels sit wholly inside another, and the union has 9,256 holes. Rutland's 15 wards from the Office for National Statistics pass all of them — sum of parts and union agree to the square metre.

Note nested_pairs. The overlaps predicate is False when one polygon lies entirely inside the other, so a check that only asks about overlaps misses every nested parcel. The two errors — double claims and holes — are independent, and a layer can have either without the other.

The three coverage rules

Each polygon is valid in all three. The set is not.

Per-feature validity against per-set coverage

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 cross?
gpd.sjoin(gdf, gdf, predicate="contains")   # does any polygon swallow another?
gdf.geometry.union_all()                    # is the merged shape hole-free?
gdf.geometry.area.sum() - union.area        # how much is claimed twice?
shapely.coverage_is_valid(gdf.geometry.values)   # overlaps and mismatched edges, GEOS 3.12+

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 — and on the Rutland parcels there was nothing for it to fix, since every one was already valid.

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.

import pandas as pd

pairs = gpd.sjoin(parcels, parcels, predicate="overlaps")
pairs = pairs[pairs.index < pairs["index_right"]]          # each pair once

# how big are they? that tells you the cause
sizes = pd.Series([
    parcels.geometry.iloc[left].intersection(parcels.geometry.iloc[right]).area
    for left, right in zip(pairs.index, pairs["index_right"])
])
print(sizes.describe(percentiles=[0.5]).map("{:,.4f}".format).to_string())
print((sizes < 1).sum(), "under 1 m²,", (sizes >= 1).sum(), "of 1 m² or more")

count       11,453.0000
mean         2,376.8146
std         30,182.8579
min              0.0000
50%              0.0029
max      1,761,624.9490
10108 under 1 m², 1345 of 1 m² or more

The size distribution is the diagnosis. Half of the overlapping pairs share less than 0.003 m² — thin slivers where two boundaries were drawn a fraction of a millimetre apart. At the other end, pairs sharing thousands of square metres are two titles genuinely claiming the same ground — a records question, not a geometry problem, and not something a script should silently resolve.

Gaps: the unassigned-land error

from shapely.geometry import Polygon

merged = parcels.geometry.union_all()
parts = getattr(merged, "geoms", [merged])
holes = [Polygon(ring) for part in parts for ring in part.interiors]
hole_areas = pd.Series([h.area for h in holes])

print(f"{len(holes)} holes, total {hole_areas.sum() / 1e6:.2f} km²")
print((hole_areas < 1).sum(), "slivers under 1 m² to close,", (hole_areas >= 1).sum(), "gaps to investigate")

9256 holes, total 1.96 km²
7974 slivers under 1 m² to close, 1282 gaps to investigate

Same diagnostic logic. A hole of a few square centimetres between two parcels is a sliver — the two boundaries were drawn from slightly different vertices. A hole of hundreds of square metres is land nobody's title covers — in a freehold register that is roads, verges and unregistered land, which is information rather than an error.

Note that "holes in the union" only finds gaps enclosed by the coverage. Land missing at the edge produces no hole — the union is simply smaller. Comparing the union against a known outer boundary catches those, and the reverse comparison finds polygons that spill out of the area:

rutland = wards.geometry.union_all()
print(f"{rutland.difference(merged).area / 1e6:.2f} km² of Rutland's {rutland.area / 1e6:.2f} km² not covered")
print(f"{merged.difference(rutland).area / 1e6:.2f} km² of parcel area outside Rutland")

32.92 km² of Rutland's 393.73 km² not covered
37.17 km² of parcel area outside Rutland

Shared edges: the cause of both

The durable fix is snapping the shared boundary onto common vertices, not patching each sliver.

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 fraction of a millimetre of each other but never coincide. This is one of the median-sized overlaps above:

by_id = parcels.set_index("INSPIREID").geometry
a, b = by_id[55832784], by_id[55835445]

print(a.touches(b), a.overlaps(b))
sliver = a.intersection(b)
print(f"{sliver.area:.4f} m² over {sliver.length / 2:.1f} m of boundary")

False True
0.0030 m² over 49.9 m of boundary

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. Spread 0.003 m² along a boundary about 50 m long and the two lines are less than a tenth of a millimetre apart on average: invisible on any map, and still enough to put a point on that line in both parcels.

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.

Two ways to store a boundary

GIS has spent forty years alternating between two storage models. In a topological model — ArcInfo coverages, GRASS vectors, PostGIS topology, TopoJSON, OpenStreetMap's nodes and ways — the boundary between two features is stored once and both reference it, so overlaps and gaps cannot arise. In a simple features model — shapefile, GeoPackage, GeoJSON, Shapely — each polygon carries its own complete ring and the shared boundary exists twice. Simple features won because every feature is independent and files are easy to write; the price is that two copies of "the same" boundary are equal only if every coordinate matches, and coordinates from different surveys, digitisers or edits essentially never do. Every coverage error on this page exists because the format allows it.

The same split applies to lines. A road network is topological when segments that meet share a node; OpenStreetMap stores it that way, so its networks are connected exactly, while a layer digitised as independent lines has undershoots and overshoots instead (cleaning a line network).

When a layer must be a coverage

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 a dissolve as holes and specks

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.

A coverage audit of Rutland parcels and wards

The data used on this page

HM Land Registry's INSPIRE Index Polygons for Rutland (21,253 freehold parcels, Open Government Licence, © Crown copyright and database rights) and Rutland's 15 wards from the ONS Wards (December 2024) full-resolution boundaries, both in British National Grid. The Land Registry download service sets a cookie before it serves the file, so it is fetched with a requests.Session.

import io
import zipfile
from pathlib import Path

import geopandas as gpd
import requests

Path("data").mkdir(exist_ok=True)
if not Path("data/parcels.gpkg").exists():
    reply = requests.Session().get(
        "https://use-land-property-data.service.gov.uk/datasets/inspire/download/Rutland.zip", timeout=300
    )
    reply.raise_for_status()
    zipfile.ZipFile(io.BytesIO(reply.content)).extract("Land_Registry_Cadastral_Parcels.gml", "data")
    gpd.read_file("data/Land_Registry_Cadastral_Parcels.gml")[["INSPIREID", "geometry"]].to_file("data/parcels.gpkg")
if not Path("data/wards.gpkg").exists():
    gpd.read_file(
        "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
        "Wards_December_2024_Boundaries_UK_BFC/FeatureServer/0/query"
        "?where=LAD24CD%3D%27E06000017%27&outFields=WD24CD,WD24NM&outSR=27700&f=geojson"
    ).to_file("data/wards.gpkg")

parcels = gpd.read_file("data/parcels.gpkg")
wards = gpd.read_file("data/wards.gpkg")
print(len(parcels), "parcels,", len(wards), "wards,", parcels.crs)

21253 parcels, 15 wards, EPSG:27700

Example 1: a full coverage audit

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

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)

    # double claims — crossing and nested — with sizes so you can tell noise from error
    claims = []
    for predicate in ("overlaps", "contains"):
        pairs = gpd.sjoin(gdf, gdf, predicate=predicate)
        pairs = pairs[pairs.index != pairs["index_right"]]
        if predicate == "overlaps":
            pairs = pairs[pairs.index < pairs["index_right"]]
        for left, right in zip(pairs.index, pairs["index_right"]):
            area = gdf.geometry.iloc[left].intersection(gdf.geometry.iloc[right]).area
            claims.append({"a": int(left), "b": int(right), "kind": predicate, "area": area})
    claims = pd.DataFrame(claims)

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

    report = {
        "features": len(gdf),
        "overlapping": int((claims["kind"] == "overlaps").sum()),
        "nested": int((claims["kind"] == "contains").sum()),
        "claim_slivers": int((claims["area"] < sliver_area).sum()),
        "claim_real": int((claims["area"] >= sliver_area).sum()),
        "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),
        "sum_of_parts_km2": round(float(gdf.geometry.area.sum()) / 1e6, 2),
        "union_km2": round(merged.area / 1e6, 2),
    }
    if boundary is not None:
        report["uncovered_at_edge_km2"] = round(boundary.difference(merged).area / 1e6, 2)
    return report, claims, holes

report, claims, holes = audit_coverage(parcels, boundary=wards.geometry.union_all())
for key, value in report.items():
    print(f"{key:>22}: {value}")

              features: 21253
           overlapping: 11453
                nested: 4106
         claim_slivers: 10114
            claim_real: 5445
                 holes: 9256
          hole_slivers: 7974
             hole_real: 1282
      sum_of_parts_km2: 428.99
             union_km2: 397.99
 uncovered_at_edge_km2: 32.92

The split between _slivers and _real is the point of the whole function. Slivers are mechanically fixable; each real double claim needs somebody to decide who owns the land. The Rutland parcels have plenty of both.

Measured in Example 1. The slivers can be closed by a script; the rest are disagreements about land that a person has to resolve.

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 = gdf.geometry.union_all().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)"
    )

assert_coverage(wards)
print("wards: coverage holds")
try:
    assert_coverage(parcels)
except AssertionError as error:
    print("parcels:", error)

wards: coverage holds
parcels: coverage broken: parts sum to 428,985,424.7 but union is 397,985,259.4 (7.789% double-counted)

One assertion, and it catches every double claim in the layer, crossing or nested, without enumerating them. It will not see gaps — a hole makes both sides smaller by the same amount — so pair it with the hole count when gaps matter. It belongs in the validation step of any pipeline whose output is a coverage.

Example 3: visualising where the problems are

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

oakham = (484000, 307500, 487000, 310000)          # a 3 x 2.5 km window on Oakham
view = parcels.cx[oakham[0]:oakham[2], oakham[1]:oakham[3]]

fig, ax = plt.subplots(figsize=(11, 9))
view.plot(ax=ax, facecolor="none", edgecolor="#94a3b8", linewidth=0.4)

real = claims[claims["area"] >= 1.0].copy()
real["geometry"] = [parcels.geometry.iloc[r.a].intersection(parcels.geometry.iloc[r.b]) for r in real.itertuples()]
real = gpd.GeoDataFrame(real, crs=parcels.crs).cx[oakham[0]:oakham[2], oakham[1]:oakham[3]]
polygons = real.explode(index_parts=False)
polygons = polygons[polygons.geom_type == "Polygon"]            # drop the stray lines and points
polygons.plot(ax=ax, facecolor="#ef4444", edgecolor="#991b1b", linewidth=0.3)

big_holes = gpd.GeoSeries([h for h in holes if h.area >= 1.0], crs=parcels.crs).cx[oakham[0]:oakham[2], oakham[1]:oakham[3]]
big_holes.plot(ax=ax, facecolor="#f59e0b", edgecolor="#b45309", linewidth=0.3)

ax.set_xlim(oakham[0], oakham[2]); ax.set_ylim(oakham[1], oakham[3])
ax.set_title(f"{len(real)} double claims (red) and {len(big_holes)} gaps (amber) of 1 m² or more")
ax.set_axis_off()
fig.savefig("coverage_errors_oakham.png", dpi=150, bbox_inches="tight")
print(ax.get_title())
print(real["kind"].value_counts().to_dict())

4229 double claims (red) and 85 gaps (amber) of 1 m² or more
{'contains': 3736, 'overlaps': 493}

A map of the errors is where the pattern shows. Most of the red is nesting: large titles whose polygons enclose whole streets of smaller ones. The amber gaps are enclosed plots that no title covers. The roads, although no freehold title covers them either, are not holes at all — they connect to the land outside the parcels, so they only show up in the comparison against the county boundary.

Measured in Example 4 below with Shapely 2.1 on GEOS 3.13.

Example 4: simplification breaks a coverage, coverage simplification does not

Every boundary in a shapefile is stored twice, once in each neighbour. Simplify each polygon on its own and the two copies of a shared edge are thinned differently:

simple = gpd.GeoDataFrame(wards[["WD24NM"]], geometry=wards.geometry.simplify(50), crs=wards.crs)
shared = gpd.GeoDataFrame(
    wards[["WD24NM"]], geometry=shapely.coverage_simplify(wards.geometry.values, 50), crs=wards.crs
)

for name, layer in [("original", wards), ("simplify(50)", simple), ("coverage_simplify(50)", shared)]:
    r = coverage_report(layer)
    vertices = shapely.get_num_coordinates(layer.geometry.values).sum()
    print(f"{name:>22}: valid coverage {r['coverage_is_valid']}, {r['overlapping_pairs']} overlapping pairs, "
          f"{r['interior_holes']} holes, {vertices} vertices")

              original: valid coverage True, 0 overlapping pairs, 0 holes, 26086 vertices
          simplify(50): valid coverage False, 29 overlapping pairs, 72 holes, 1208 vertices
 coverage_simplify(50): valid coverage True, 0 overlapping pairs, 0 holes, 2083 vertices

coverage_simplify (Shapely 2.1, GEOS 3.12+) simplifies each shared edge once and gives both neighbours the same result.

Reprojection, by contrast, does not break a coverage on its own: each vertex is transformed by the same formula, so two copies of the same vertex land on the same new coordinate.

for crs in ["EPSG:4326", "EPSG:3857"]:
    there = wards.to_crs(crs)
    print(crs, shapely.coverage_is_valid(there.geometry.values),
          shapely.coverage_is_valid(there.to_crs(wards.crs).geometry.values))

EPSG:4326 True True
EPSG:3857 True True

Coverage errors have a size, and the size says who fixes them

The size of the error decides whether a script may fix it.

Rutland's parcels sum to 7.8% more land than they cover, and the overlap sizes split cleanly in two. Half of the overlapping pairs share less than 0.003 m²: one fence line digitised twice, a fraction of a millimetre apart. Those are artefacts of storing the boundary twice, and they can be closed mechanically by snapping one copy onto the other (finding and fixing gaps and overlaps). Pairs sharing thousands of square metres are two titles claiming the same ground, and 4,106 parcels sit wholly inside another. Those are disagreements about the world. Resolving them by giving the land to whichever polygon comes first in the file is worse than leaving them, because it destroys the evidence that a disagreement existed.

Measuring and cleaning a coverage once is straightforward in Python. Maintaining one — so that editing a parcel moves its neighbour's edge too — needs a topological store: PostGIS topology, GRASS, or QGIS with topological editing. If boundaries are edited repeatedly, that is the right tool.

Gaps carry the same split with one more twist: double claims are always errors in a coverage; gaps sometimes are not. A hole where a road reserve or an unregistered plot lies is real, and the county's roads do not even appear as holes, because they connect to land outside the parcels. That is why the audit reports slivers, real overlaps and gaps separately rather than as one "problems" count.

Non-coverages, containment and coverage_union_all

Buffers, catchments and ranges are not coverages. Overlap is their normal state; do not audit them for it.
overlaps misses containment. A polygon wholly inside another is contains/within, not overlaps; 4,106 Rutland parcels are in that position. The sum-of-parts check sees both.
union_all is not expensive at this size — about 2.5 s for the 21,253 parcels here — but it grows with the layer. The area-sum assertion needs the union too; shapely.coverage_is_valid does not, and took under a second.
overlaps can be True with a zero-area intersection. Nine Rutland pairs of valid parcels test as overlapping, yet intersection returns only lines: the predicate is computed exactly, while the overlay snaps a crossing too thin to survive. Filter pairs by intersection area, not by the predicate alone.
coverage_union_all assumes a valid coverage. On the Rutland parcels it stops with TopologyException: side location conflict; use union_all until the layer passes coverage_is_valid.
simplify() reliably breaks a coverage. Each boundary is simplified separately, so shared edges diverge — Example 4. Use shapely.coverage_simplify, topojson, or PostGIS ST_CoverageSimplify.
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; the Rutland download reaches 37 km² beyond the county.
Dissolve does not repair coverage. It merges by attribute and carries slivers through into the result as holes.

Internal links

How to find and fix gaps and overlaps in a polygon coverage — the repair procedure
What makes a geometry valid? The OGC rules explained — the per-feature check this page contrasts with
Spatial predicates explained — why touches flipping to overlaps is the signature
How to snap and align geometries to fix slivers and gaps — the snapping approach
How to dissolve polygons in Python (GeoPandas) — the operation that exposes coverage errors
Spatial data quality: the six dimensions that matter — where consistency sits
Coordinate precision and floating point in GIS explained — why identical-looking edges are not identical

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?

shapely.coverage_is_valid(gdf.geometry.values) on Shapely 2.1 with GEOS 3.12 or later: under a second for Rutland's 21,253 parcels; pass gap_width to have it flag narrow gaps as well. Without it, compare gdf.geometry.area.sum() against gdf.geometry.union_all().area; any difference is area claimed twice. Neither tells you where; shapely.coverage_invalid_edges does.

Are gaps always errors?

No. A hole where a river, road reserve or unregistered plot sits is real; in Rutland's freehold parcels the roads are uncovered land of this kind, although they run out to the edge and so count as missing coverage rather than as holes. 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. All 21,253 Rutland parcels are valid, and they add up to 31 km² more land than they cover.

Does reprojecting break a coverage?

Not by itself. Every vertex goes through the same transformation, so two copies of a shared vertex land on the same new coordinate: Rutland's wards stayed a valid coverage through EPSG:4326 and EPSG:3857 and back. If a coverage breaks around a reprojection, look for the step that treats each feature on its own — simplify, densifying, or snapping each feature to a grid separately — and check coverage_is_valid after each step.

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.