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

On this page

Problem statement

Registered land titles should not overlap: a piece of ground is in one freehold title or none. HM Land Registry publishes the index polygon of every freehold title in England and Wales, and here are the 7,304 of them that lie in one Leeds ward, Headingley & Hyde Park (the setup is under The data used on this page):

from shapely import union_all

print(f"{len(titles):,} titles, all valid: {titles.is_valid.all()}")
print(f"sum of title areas: {titles.area.sum():>12,.1f} m²")
print(f"area of their union: {union_all(titles.geometry.values, grid_size=0.001).area:>11,.1f} m²")
7,304 titles, all valid: True
sum of title areas:  2,170,942.6 m²
area of their union: 2,156,007.2 m²

Fifteen thousand square metres are counted twice. Somewhere in 7,304 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 coarse grid moves real boundaries, and on this data — as step 2 shows — barely closes anything anyway.

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, fix the small ones, report the big ones:

import numpy as np
import shapely
from shapely import Polygon

GRID = 0.001        # m — the data's own precision: every coordinate is a whole millimetre
SLIVER = 1.0        # m² — below this, treat it as a digitising artefact


def coverage_errors(gdf):
    """Every overlap with positive area, and every hole in the union."""
    a, b = gdf.sindex.query(gdf.geometry.values, predicate="intersects")
    a, b = a[a < b], b[a < b]                               # each pair once
    shapes = shapely.intersection(gdf.geometry.values[a], gdf.geometry.values[b], grid_size=GRID)
    overlaps = gpd.GeoDataFrame({"a": a, "b": b}, geometry=shapes, crs=gdf.crs)
    overlaps["area"] = overlaps.area
    overlaps = overlaps[overlaps["area"] > 0].reset_index(drop=True)

    union = shapely.union_all(gdf.geometry.values, grid_size=GRID)
    rings = [r for p in shapely.get_parts(union) for r in p.interiors]
    holes = gpd.GeoDataFrame(geometry=[Polygon(r) for r in rings], crs=gdf.crs)
    holes["area"] = holes.area
    return overlaps, holes


overlaps, holes = coverage_errors(titles)
print(f"{len(overlaps)} overlaps: {(overlaps['area'] < SLIVER).sum()} under 1 m², "
      f"{overlaps['area'].sum():,.1f} m² in all")
print(f"{len(holes)} holes: {(holes['area'] < SLIVER).sum()} under 1 m²")
print(f"overlap sizes (m²), median {overlaps['area'].median():.3f}, "

      f"90th percentile {overlaps['area'].quantile(0.9):.2f}, largest {overlaps['area'].max():,.1f}")

# a sliver's mean width is about twice its area over its perimeter
small = overlaps[overlaps["area"] < SLIVER]
width_mm = 2 * small.area / small.length * 1000
print(f"small overlaps: median width {width_mm.median():.1f} mm, 90th percentile {width_mm.quantile(0.9):.1f} mm")
572 overlaps: 417 under 1 m², 21,840.5 m² in all
474 holes: 327 under 1 m²
overlap sizes (m²), median 0.025, 90th percentile 156.35, largest 1,154.4
small overlaps: median width 0.9 mm, 90th percentile 11.8 mm

Most of the errors are tiny — the median overlap is a few hundredths of a square metre, a strip about a millimetre wide where two neighbouring titles were drawn a hair apart. A handful are over a thousand square metres. They are not the same problem and must not get the same fix.

Symptom Size What it is What to do
overlap < 1 m² two neighbours digitised the same line slightly differently clip it out of one side
overlap > 1 m² two records claiming the same ground report, do not resolve
hole < 1 m² sliver between mismatched edges give it to a neighbour
hole > 1 m² genuinely unassigned — here, roads and unregistered land investigate
missing at the edge any a polygon that was never supplied compare against a boundary

Artefacts versus disagreements

Coverage errors in 7,304 Leeds land titles split by size into the ones a script may fix and the ones it may not. one threshold decides whether a script may actMeasuredAction417 overlaps under 1 m²clip out of the higher-numberedtitle155 overlaps over 1 m²write to a review file, decide byhand327 holes under 1 m²give to the longest-shared-edgeneighbour147 holes over 1 m²roads and unregistered land —investigateThree passes fixed every small one and moved 0.87 m² in all spatialworkflow.io
One threshold decides whether a script may act. Everything above it is a person's decision.

Valid first, precision, measure, clip overlaps, fill slivers

The order of operations for repairing a polygon coverage. check, measure, fix the small, report the large1make_valid everything firstor measurements lie2check the coordinates’ precisionsnap only if they are noisy3measure with intersects + areaoverlaps misses contained titles4clip small overlaps, fill sliversrepeat until a pass finds nothing5export the rest for reviewdo not resolveOn these titles every coordinate was already a whole millimetre spatialworkflow.io
Snapping is a check, not a step: it only helps when the coordinates carry floating-point noise.

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.

gdf = titles.copy()
invalid = (~gdf.is_valid).sum()
print(f"{invalid} invalid geometries")
if invalid:
    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()
0 invalid geometries

These titles are clean polygon by polygon. Everything wrong with them is between polygons.

2. Check what precision the data carries — and whether snapping can help

The usual advice is to snap every coordinate to a grid with shapely.set_precision, so that two copies of a shared boundary that differ in the ninth decimal place round to the same number. First look at what precision the data has:

xy = shapely.get_coordinates(titles.geometry.values)
on_mm = np.isclose(xy * 1000, np.round(xy * 1000), atol=1e-6).mean()
print(f"{len(xy):,} coordinates, {on_mm:.1%} already on a whole millimetre")
100,580 coordinates, 100.0% already on a whole millimetre

Every coordinate is already on a millimetre grid, so there is no floating-point noise for a 1 mm snap to remove. The slivers are not rounding error; they are two people putting the same fence a millimetre or so apart, and in one case in ten more than a centimetre. Snapping to a coarser grid does not fix that either, because two vertices 7 mm apart round to different grid points unless the grid is coarse enough to swallow real geometry:

for grid in (0.001, 0.01, 0.1, 0.5):
    snapped = titles.assign(geometry=shapely.set_precision(titles.geometry.values, grid))
    snapped = snapped[~snapped.is_empty].reset_index(drop=True)

    ov, ho = coverage_errors(snapped)
    print(f"grid {grid:>5} m: {len(ov):>3} overlaps, {len(ho):>3} holes, "
          f"{len(titles) - len(snapped)} titles collapsed to nothing")
grid 0.001 m: 572 overlaps, 475 holes, 0 titles collapsed to nothing
grid  0.01 m: 882 overlaps, 768 holes, 0 titles collapsed to nothing
grid   0.1 m: 865 overlaps, 701 holes, 1 titles collapsed to nothing
grid   0.5 m: 669 overlaps, 607 holes, 3 titles collapsed to nothing

A 1 cm grid makes things worse — rounding two nearly-coincident edges independently pushes some of them across each other. Even at half a metre, a grid that visibly moves every boundary, hundreds of overlaps survive and titles start to vanish. Snapping to a grid only closes gaps between vertices that were meant to be identical and nearly are. It is the right fix for coordinates that went through a reprojection or a text export; it is the wrong tool for digitising disagreements, which need the neighbour-aware fixes in steps 4 and 5.

Overlapping pairs among 7,304 Leeds land titles after snapping to each grid size. overlaps left after snapping 7,304 titles to a gridas supplied5721 mm grid572 — no change1 cm grid882 — worse10 cm grid865, 1 title gone50 cm grid669, 3 titles goneclip + fill155, all over 1 m²0882No grid gets near what three passes of clip-and-fill do spatialworkflow.io
Snapping to a coarser grid first makes the overlaps worse, then starts deleting titles. The clip-and-fill bar is from example 1.

3. Measure what is actually there

coverage_errors from the quick answer does the measuring, and its first line is the one most versions of it get wrong. The obvious predicate for "which polygons overlap" is overlaps:

g = titles.geometry.values
for predicate in ("overlaps", "intersects"):
    a, b = titles.sindex.query(g, predicate=predicate)
    a, b = a[a < b], b[a < b]
    area = shapely.area(shapely.intersection(g[a], g[b], grid_size=GRID))
    print(f"{predicate:10s} {len(a):>6,} pairs, {(area > 0).sum():>4} sharing area, {area.sum():>9,.1f} m²")


a, b = titles.sindex.query(g, predicate="contains")
a, b = a[a != b], b[a != b]
same = shapely.equals(g[a], g[b])
print(f"one title wholly inside another: {(~same).sum()} pairs; exact duplicates: {same.sum() // 2} pairs")
overlaps      980 pairs,  470 sharing area,  10,117.0 m²
intersects 11,184 pairs,  572 sharing area,  21,840.5 m²
one title wholly inside another: 55 pairs; exact duplicates: 47 pairs

In the DE-9IM sense that the predicate implements, two polygons overlap only when each has some area outside the other. A title lying wholly inside another, or an exact duplicate, is not an overlap by that definition — and in this ward those account for more than half of all the double-claimed area. intersects catches every pair that shares anything; filtering on the intersection's area then drops the thousands of neighbours that only share an edge.

Computing every overlay with grid_size keeps the results on the data's own millimetre grid, so a repaired polygon does not pick up coordinates of 429672.33000000007 and the next comparison does not see a new sliver.

4. Resolve small overlaps by clipping one side

For sub-threshold overlaps, one polygon must give up the disputed strip. Which one is arbitrary, so the rule must be deterministic:

def polygonal(geom):
    """Keep the polygon parts of an overlay result that also produced lines."""
    if geom.geom_type in ("Polygon", "MultiPolygon"):
        return geom
    parts = [p for p in shapely.get_parts(geom) if p.geom_type == "Polygon"]
    return shapely.MultiPolygon(parts) if len(parts) > 1 else (parts[0] if parts else Polygon())


def resolve_small_overlaps(gdf, overlaps, max_area=SLIVER):
    """Clip each small overlap out of the higher-numbered feature."""
    geoms = gdf.geometry.values.copy()
    small = overlaps[overlaps["area"] < max_area]
    for a, b in zip(small["a"], small["b"]):                 # a < b always
        geoms[b] = polygonal(shapely.difference(geoms[b], geoms[a], grid_size=GRID))
    return gdf.assign(geometry=geoms), len(small)


gdf, clipped = resolve_small_overlaps(gdf, overlaps)
overlaps, holes = coverage_errors(gdf)
print(f"clipped {clipped}; now {len(overlaps)} overlaps ({(overlaps['area'] < SLIVER).sum()} small), "
      f"{len(holes)} holes")
clipped 417; now 158 overlaps (3 small), 473 holes

Two details are there because the naive version failed. It subtracts the whole neighbour rather than the measured overlap shape, because the overlap shape was computed separately and its edges do not coincide exactly with either polygon's — subtracting it leaves a hairline. And polygonal() exists because a difference on a fixed grid can return a GeometryCollection holding a polygon and a stray line, and the next overlay call then refuses it with GEOSException: Overlay input is mixed-dimension.

Using b — the higher index — 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.

5. Fill sliver holes by assigning them to a neighbour

A sliver between two titles belongs to one of them. The usual rule: give it to the neighbour it shares the longest boundary with.

def fill_slivers(gdf, holes, max_area=SLIVER):
    """Union each small hole into the neighbour sharing most of its edge."""
    geoms = gdf.geometry.values.copy()
    filled = 0
    for sliver in holes.geometry[holes["area"] < max_area]:
        near = gdf.sindex.query(sliver, predicate="intersects")
        if len(near) == 0:
            continue
        shared = shapely.length(shapely.intersection(geoms[near], sliver.boundary, grid_size=GRID))
        winner = near[np.argmax(shared)]                       # ties go to the lower index
        geoms[winner] = polygonal(shapely.union(geoms[winner], sliver, grid_size=GRID))
        filled += 1
    return gdf.assign(geometry=geoms), filled


gdf, filled = fill_slivers(gdf, holes)
overlaps, holes = coverage_errors(gdf)
print(f"filled {filled}; now {len(overlaps)} overlaps ({(overlaps['area'] < SLIVER).sum()} small), "
      f"{len(holes)} holes ({(holes['area'] < SLIVER).sum()} small)")
print(f"total area moved: {gdf.area.sum() - titles.area.sum():+.2f} m² across {len(gdf):,} titles")
filled 326; now 162 overlaps (7 small), 148 holes (1 small)
total area moved: -0.87 m² across 7,304 titles

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 title rather than being deleted. Filling can re-create a few tiny overlaps where a sliver touched a neighbour's corner; a second clip pass removes them, and example 1 runs one.

Overlaps plus holes in 7,304 Leeds land titles through each repair step. overlaps + holes, Headingley & Hyde Parkfound1,046after clipping 417 overlaps631÷ 2after filling 326 slivers310every one over 1 m²after passes 2 and 3302The 302 left are the review file: 155 overlaps and 147 holes spatialworkflow.io
Two thirds of the errors were slivers a script may remove. The rest shrink by eight in two more passes and then stop — they were never slivers.

6. Report the rest — do not resolve it

review = overlaps[overlaps["area"] >= SLIVER].copy()
review["title_a"] = titles["INSPIREID"].values[review["a"]]
review["title_b"] = titles["INSPIREID"].values[review["b"]]
review.to_file("review_overlaps.gpkg", driver="GPKG")

print(f"⚠ {len(review)} overlaps over 1 m² need a human decision; the largest:")
for r in review.nlargest(3, "area").itertuples():
    print(f"    titles {r.title_a} and {r.title_b} share {r.area:,.1f} m²")
print(f"{(holes['area'] >= SLIVER).sum()} holes over 1 m², the largest {holes['area'].max():,.0f} m²")
⚠ 155 overlaps over 1 m² need a human decision; the largest:
    titles 43968482 and 51430302 share 1,154.4 m²
    titles 26314979 and 43968482 share 1,103.8 m²
    titles 26314979 and 26317727 share 1,101.5 m²
147 holes over 1 m², the largest 3,971 m²

The large holes are roads and land that has never been registered; the large overlaps are titles whose index polygons claim the same ground. Neither is a geometry problem. Writing them to a file that somebody opens in QGIS is the deliverable. A script that silently picks a winner for a thousand-square-metre overlap has destroyed the only evidence that there was a question to ask.

Repairing Headingley & Hyde Park titles, with a report

The data used on this page

HM Land Registry's INSPIRE Index Polygons for Leeds — the extent of every registered freehold title — cut to the titles lying wholly inside one electoral ward, Headingley & Hyde Park (ONS ward boundaries, December 2024). British National Grid, metres.

import sys

import geopandas as gpd
import pyogrio
import requests

ONS = "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services"
wards = gpd.read_file(f"{ONS}/Wards_December_2024_Boundaries_UK_BFC/FeatureServer/0/query"
                      "?where=LAD24NM%3D%27Leeds%27&outFields=WD24CD,WD24NM&outSR=27700&f=geojson")
ward = wards.set_crs(27700, allow_override=True).query("WD24NM == 'Headingley & Hyde Park'").geometry.iloc[0]

# HM Land Registry's download sets a cookie and redirects to itself, so use a session
reply = requests.Session().get("https://use-land-property-data.service.gov.uk/datasets/inspire/download/Leeds.zip",
                               timeout=600)
reply.raise_for_status()
with open("leeds_inspire.zip", "wb") as f:
    f.write(reply.content)
pyogrio.set_gdal_config_options({"GML_DOWNLOAD_SCHEMA": "NO"})   # the file's schema URL is unreachable
leeds = gpd.read_file("zip://leeds_inspire.zip!Land_Registry_Cadastral_Parcels.gml", columns=["INSPIREID"])
titles = leeds[leeds.within(ward)].sort_values("INSPIREID").reset_index(drop=True)
titles.to_parquet("titles.parquet")
print(f"{len(leeds):,} Leeds titles, {len(titles):,} inside the ward")
293,729 Leeds titles, 7,304 inside the ward

Example 1: the complete repair, with a report


def repair_coverage(gdf, *, sliver=SLIVER, max_passes=5):
    report = {"features_in": len(gdf), "area_in": round(gdf.area.sum(), 1)}
    gdf = gdf.reset_index(drop=True)

    bad = int((~gdf.is_valid).sum())
    if bad:
        gdf = gdf.assign(geometry=gdf.geometry.make_valid())
    report["repaired_invalid"] = bad

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


    # clipping can open a hairline hole and filling can touch a neighbour's corner,
    # so repeat until a pass finds nothing small left to do
    for n in range(1, max_passes + 1):
        small = (overlaps["area"] < sliver).sum() + (holes["area"] < sliver).sum()
        if not small:
            break
        gdf, clipped = resolve_small_overlaps(gdf, overlaps, sliver)
        overlaps, holes = coverage_errors(gdf)
        gdf, filled = fill_slivers(gdf, holes, sliver)
        overlaps, holes = coverage_errors(gdf)
        report[f"pass_{n}"] = f"clipped {clipped}, filled {filled}"
    report["overlaps_remaining"] = len(overlaps)
    report["small_overlaps_remaining"] = int((overlaps["area"] < sliver).sum())
    report["holes_remaining"] = len(holes)
    report["small_holes_remaining"] = int((holes["area"] < sliver).sum())
    report["area_change_m2"] = round(gdf.area.sum() - report.pop("area_in"), 3)
    report["features_out"] = int((~gdf.is_empty).sum())
    return gdf, report, overlaps, holes


clean, report, review_overlaps, review_holes = repair_coverage(titles)
for k, v in report.items():
    print(f"{k:26s} {v}")
features_in                7304
repaired_invalid           0
overlaps_found             572
holes_found                474
pass_1                     clipped 417, filled 326
pass_2                     clipped 7, filled 3
pass_3                     clipped 1, filled 1
overlaps_remaining         155
small_overlaps_remaining   0
holes_remaining            147
small_holes_remaining      0
area_change_m2             -0.866
features_out               7304

features_in == features_out and a small area_change_m2 are the two numbers to check. If the feature count changed, something clipped a title away entirely; if the area moved by thousands, the threshold is too high.

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 always gives a true coverage:

from shapely.ops import polygonize


def rebuild_from_edges(gdf):
    """Node every boundary, re-polygonise, and re-attach attributes."""
    edges = shapely.union_all(gdf.geometry.boundary.values, grid_size=GRID)   # nodes every crossing
    faces = gpd.GeoDataFrame(geometry=list(polygonize(shapely.get_parts(edges))), crs=gdf.crs)

    points = faces.assign(geometry=faces.representative_point())
    hits = gpd.sjoin(points, gdf, predicate="within", how="left")
    claims = hits.groupby(level=0)["INSPIREID"].count()
    faces["INSPIREID"] = hits.groupby(level=0)["INSPIREID"].first()
    return faces, claims


faces, claims = rebuild_from_edges(titles)
print(f"{len(titles):,} titles → {len(faces):,} faces")
print(f"faces claimed by no title: {(claims == 0).sum():,}, by one: {(claims == 1).sum():,}, "
      f"by two or more: {(claims > 1).sum():,}")
7,304 titles → 8,356 faces
faces claimed by no title: 483, by one: 7,223, by two or more: 650

polygonize builds faces from a noded line network, so gaps and overlaps are structurally impossible in the output. The cost shows in the counts: every hole and every overlap becomes a face of its own, the unclaimed faces have no title to inherit from, and each face in a former overlap gets its attributes from whichever title the first() happens to name. Use it when the alternative is worse.

Example 3: guarding against regression

%%writefile test_coverage.py
import geopandas as gpd
import pytest
import shapely

GRID = 0.001


def small_holes_and_overlaps(gdf, limit=1.0):
    a, b = gdf.sindex.query(gdf.geometry.values, predicate="intersects")
    a, b = a[a < b], b[a < b]
    areas = shapely.area(shapely.intersection(gdf.geometry.values[a], gdf.geometry.values[b], grid_size=GRID))
    union = shapely.union_all(gdf.geometry.values, grid_size=GRID)
    holes = [shapely.area(shapely.Polygon(r)) for p in shapely.get_parts(union) for r in p.interiors]
    return int(((areas > 0) & (areas < limit)).sum()), sum(h < limit for h in holes)


@pytest.fixture(scope="module")
def raw():
    return gpd.read_parquet("titles.parquet")


@pytest.fixture(scope="module")
def cleaned():
    return gpd.read_parquet("titles_clean.parquet")


def test_raw_titles_are_not_a_coverage(raw):
    overlaps, holes = small_holes_and_overlaps(raw)
    assert overlaps > 0 and holes > 0          # the fixture really is broken


def test_cleaned_output_has_no_slivers(cleaned):
    overlaps, holes = small_holes_and_overlaps(cleaned)
    assert overlaps == 0, f"{overlaps} overlaps under 1 m²"
    assert holes == 0, f"{holes} holes under 1 m²"


def test_cleaning_moved_almost_no_area(raw, cleaned):
    assert abs(cleaned.area.sum() - raw.area.sum()) < 10
clean.to_parquet("titles_clean.parquet")
!{sys.executable} -m pytest -q --no-header -p no:cacheprovider --color=no --tb=short test_coverage.py
...                                                                      [100%]
3 passed in 4.40s

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.

Why snapping did almost nothing here

The two ways two stored copies of a shared boundary can disagree, and which one snapping fixes. two ways the two copies of an edge disagreeNoise — snapping fixes it✓copies started identical✓differ at 1e-9 m✓same vertices, nudged✓both round to one grid point✓the sliver ceases to existDisagreement — it cannot✕copies drawn separately✕differ by ~1 mm to over 1 cm✕different vertices on the fence✕round to different grid points✕clip or fill insteadLand Registry title polygons are the right-hand kind spatialworkflow.io
Snapping makes two copies of an edge agree only when they were meant to be identical and nearly are.

Snapping is the fix most often recommended for coverage errors, and step 2 showed it doing almost nothing here. The reason is worth understanding, because it tells you in advance which fix a layer needs.

In a shapefile or GeoPackage, the boundary between two parcels is stored twice — once in each feature — as two independent coordinate lists. There are two ways those copies can disagree.

The first is noise: the copies started identical, and something nudged one of them — a reprojection round trip, a simplify, an export through a format with less precision. The vertices are in the same places to within 10⁻⁹ metres. set_precision rounds every coordinate in the layer to the same grid, both copies land on the same grid points, and the overlap does not get clipped away so much as cease to exist. For this kind of error, snapping is exactly right.

The second is disagreement: the copies were drawn separately, by different people or at different times, against different base maps. The Land Registry's index polygons are this kind — every coordinate is a whole millimetre, and neighbouring titles put the same fence about a millimetre apart at the median and more than a centimetre apart in one sliver in ten, with different vertices along it. There is no pair of nearly-identical vertices for a grid to merge. Rounding each copy independently can even push two nearly-coincident edges across each other, which is why the 1 cm grid made things worse. This kind of error needs a neighbour-aware fix: decide which feature owns the strip, and clip or fill accordingly.

The quickest way to tell them apart is the precision check in step 2. Coordinates with long, ragged decimal tails have been computed; snapping will probably help. Coordinates already on a round grid were captured that way; snapping will not.

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.025 m² overlap between neighbouring titles is a drawing artefact — no surveyor intended it, and no information is lost by removing it. A 1,154 m² overlap is two records asserting incompatible things about the ground, and any automatic resolution silently picks a winner in a question 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.

Choosing between snapping and a neighbour-aware repair from the coordinates themselves. what do the coordinates look like?Are the coordinates already on a roundgrid?ragged decimalsNoise — snap1 mm grid, then re-measureround gridDisagreement — clip and filldeterministic, repeat to a fixpointlarge errorsReport, do not fixover the sliver thresholdEvery Land Registry coordinate here was a whole millimetre — the middle branch spatialworkflow.io
Look at the coordinates before choosing the fix: they tell you whether the copies were computed or drawn.

set_precision, polygonize and representative points

  • 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. On this page, clipping and filling moved 0.87 m² in total across 7,304 titles.
  • 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 check expensive. Use the area-sum check from the problem statement to decide whether it is worth running at all; intersects on 7,304 titles found 11,184 candidate pairs, almost all of them neighbours that only share an edge.
  • shapely.coverage_is_valid (Shapely 2.1, GEOS 3.12+) answers the yes/no question in one call — it returns False for these titles — and shapely.coverage_invalid_edges returns the offending edges to map.
  • GRASS v.clean (with tool=snap,break,rmdupl,rmsa) does a comparable job interactively and is worth reaching for on a one-off. QGIS "Fix geometries" is not the equivalent: it repairs each polygon on its own and will not touch an overlap between two valid ones.

FAQ

What grid size should I snap to?

First check whether snapping will help at all: if every coordinate is already on a round grid, as the Land Registry's titles are, the errors are digitising disagreements and no grid will fix them. If the coordinates have ragged decimal tails, snap to the smallest grid that closes the noise — 1 mm for metric data is usually enough — and never to one 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 half a metre it collapsed three of these titles to nothing. 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 did overlaps miss some of the double claims?

Because in DE-9IM terms two polygons overlap only when each has area outside the other. A title wholly inside another, or an exact duplicate, is not an overlap — in this ward those were more than half of the double-claimed area. Query with intersects and keep the pairs whose intersection has area.

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.