How to Snap and Align Geometries to Fix Slivers and Gaps in Python

On this page

Adjacent polygons that should share an edge often don't. After an overlay, a union, or a reprojection you end up with a layer that looks fine at full extent but, when you zoom in, is riddled with hairline overlaps and gaps. These artefacts break dissolves, inflate feature counts, and make area statistics wrong. This guide finds and fixes them in real data — Rutland's Land Registry parcels and the ONS ward boundaries of Oakham, two layers drawn independently by two organisations — with shapely and geopandas.

Problem statement

You are dealing with this article if you see any of these symptoms:

Thin sliver polygons appear after an overlay() or union — long, skinny shapes with almost no area that sit along the boundary between two real features.
Gaps open up between neighbours that are supposed to be coincident — small slivers of empty space where two polygons fail to touch.
Vertices that almost-but-not-quite match — two boundaries that should be identical differ by a fraction of a millimetre, or by a metre, so the geometries neither snap together nor dissolve cleanly.
A dissolve() that should merge two adjacent polygons into one leaves a visible internal boundary or a tiny hole.
Feature counts jump after an intersection because every shared edge spawns one or two extra micro-polygons.

Quick answer

Work in a projected, metric CRS so tolerances and areas are in metres. Find slivers by area and shape, decide which near-misses are noise and which are real, align the layers where one is authoritative, and drop or absorb what is left:

Aligning vertices to a shared grid turns near-miss boundaries into a clean shared edge.

pieces = gpd.overlay(parcels, wards, how="intersection").explode(index_parts=False)
area, thinness = pieces.area, pieces.area / pieces.length ** 2
is_sliver = area < 5                                          # m², in British National Grid
print(f"{len(pieces):,} pieces from {len(parcels):,} parcels; {int(is_sliver.sum())} under 5 m², "
      f"{area[is_sliver].sum():,.0f} m² between them")

cleaned = pieces[~is_sliver]
print(f"{len(cleaned):,} pieces kept, {cleaned.area.sum() / pieces.area.sum():.4%} of the area")

5,147 pieces from 5,006 parcels; 164 under 5 m², 113 m² between them
4,983 pieces kept, 99.9993% of the area

Detect slivers, snap, then close the gaps

A gap of a few centimetres is invisible on screen and fatal to a dissolve.

Identify slivers by area and thinness

A sliver is small and thin. Filtering on area alone catches some, but a long thin sliver can have a deceptively large area. The thinness ratio $\text{area}/\text{perimeter}^2$ helps: compact shapes approach 1 / (4*pi) (~0.0796), slivers approach zero.

small, thin = area < 5, thinness < 0.01
print(f"small: {int(small.sum())}, thin: {int(thin.sum())}, both: {int((small & thin).sum())}, either: {int((small | thin).sum())}")
thin_only = pieces[thin & ~small]
print(f"thin but not small: {len(thin_only)} pieces, the largest {thin_only.area.max():,.0f} m², "
      f"median {thin_only.area.median():,.0f} m²")

small: 164, thin: 209, both: 122, either: 251
thin but not small: 87 pieces, the largest 75,084 m², median 1,152 m²

Thinness catches long artefacts that an area filter misses — and it also catches long real parcels. The largest "thin" piece here is several hectares: a road or track parcel, long and narrow by nature. Use thinness to find candidates, and an area threshold (or a look) to decide.

Align vertices with `set_precision` — for noise only

shapely.set_precision snaps every coordinate to a fixed grid. When two boundaries differ only by floating-point noise, rounding both to the same grid makes them identical:

def sliver_count(p, w):
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        result = gpd.overlay(p, w, how="intersection").explode(index_parts=False)
    return len(result), int((result.area < 5).sum())

aligned_parcels, aligned_wards = parcels.copy(), wards.copy()
aligned_parcels["geometry"] = shapely.set_precision(parcels.geometry.values, grid_size=0.001)
aligned_wards["geometry"] = shapely.set_precision(wards.geometry.values, grid_size=0.001)
print("original:       %s pieces, %s under 5 m²" % sliver_count(parcels, wards))
print("1 mm grid both: %s pieces, %s under 5 m²" % sliver_count(aligned_parcels, aligned_wards))

original:       5147 pieces, 164 under 5 m²
1 mm grid both: 5057 pieces, 74 under 5 m²

A millimetre grid removed 90 of the 164 small pieces: slivers less than a millimetre wide, which collapsed to lines on the grid and were dropped by the overlay. Those were floating-point noise. The 74 that remain are not: the Land Registry draws title extents, the ONS draws electoral boundaries, and where one follows a hedge and the other a road centre line they are metres apart. No grid that small reconciles them, and a grid large enough to would distort everything else.

Snap one layer to another with `shapely.snap`

When one layer is authoritative, move the other's vertices onto it. shapely.snap(geom, reference, tolerance) moves vertices of geom that fall within tolerance of the reference:

ward_edges = wards.boundary.union_all()
for tolerance in (1, 3):
    snapped = parcels.copy()
    snapped["geometry"] = snapped.geometry.apply(lambda g: shapely.snap(g, ward_edges, tolerance)).make_valid()
    snapped = snapped[snapped.geom_type.isin(["Polygon", "MultiPolygon"])]
    print(f"parcels snapped to ward edges within {tolerance} m: %s pieces, %s under 5 m²" % sliver_count(snapped, wards))

parcel_edges = parcels.boundary.union_all()
moved = wards.copy()
moved["geometry"] = moved.geometry.apply(lambda g: shapely.snap(g, parcel_edges, 3)).make_valid()
print("wards snapped to parcel edges within 3 m: %s pieces, %s under 5 m²" % sliver_count(parcels, moved))

parcels snapped to ward edges within 1 m: 5104 pieces, 114 under 5 m²
parcels snapped to ward edges within 3 m: 5122 pieces, 117 under 5 m²
wards snapped to parcel edges within 3 m: 5413 pieces, 330 under 5 m²

Snapping helped, modestly. It only moves vertices: where a parcel edge runs alongside the ward edge with no vertex near it, nothing moves, and snapping the other way — the wards to the parcels — pulled ward boundaries onto whichever parcel edge was nearest and made things worse, as the last line shows. Snap toward the layer you trust, with a tolerance below the size of anything real.

The alignment tests above, counted.

Remove sliver polygons after overlay

Overlays are the classic sliver factory. After an intersection, filter the result on area — having exploded multi-part features so each part is judged on its own:

before = len(pieces)
kept = pieces[pieces.area >= 5].reset_index(drop=True)
print(f"{before - len(kept)} pieces under 5 m² dropped; {len(kept):,} kept")

164 pieces under 5 m² dropped; 4,983 kept

Close gaps between adjacent polygons

Gaps are the inverse of overlaps. The parcels themselves have them — holes in their union where no title covers the ground:

merged = parcels.union_all()
holes = lambda g: [shapely.Polygon(r) for part in shapely.get_parts(g) for r in part.interiors]
print(f"{len(holes(merged))} holes in the parcel coverage, {sum(h.area for h in holes(merged)):,.0f} m² in all")

per_feature = parcels.geometry.buffer(0.05, join_style="mitre").buffer(-0.05, join_style="mitre").union_all()
print(f"closing each parcel, then merging: {len(holes(per_feature))} holes")
for tolerance in (0.01, 0.05, 0.5):
    closed = merged.buffer(tolerance, join_style="mitre").buffer(-tolerance, join_style="mitre")
    print(f"closing the merged coverage at {tolerance:>4} m: {len(holes(closed)):>3} holes, "
          f"{sum(h.area for h in holes(closed)):,.0f} m², area added {closed.area - merged.area:,.1f} m²")

390 holes in the parcel coverage, 55,389 m² in all
closing each parcel, then merging: 386 holes
closing the merged coverage at 0.01 m: 106 holes, 55,385 m², area added 4.8 m²
closing the merged coverage at 0.05 m:  68 holes, 55,376 m², area added 15.1 m²
closing the merged coverage at  0.5 m:  51 holes, 55,647 m², area added 205.8 m²

Two lessons. A buffer-out-buffer-in "closing" applied to each feature separately does nothing to the gaps between features — each parcel is restored to itself before they meet. Applied to the merged coverage, it removed most holes at a centimetre, adding a few square metres. What it did not remove are the large holes, and it should not: they are land with no registered title in this extract — roads and verges, among others — and they are real.

Verify the result

Always confirm the cleanup worked: no remaining tiny parts, every geometry valid, and total area roughly unchanged:

print("all valid:", bool(kept.is_valid.all()), "| smallest piece:", round(kept.area.min(), 1), "m² | empty:",
      int(kept.is_empty.sum()), f"| area kept: {kept.area.sum() / pieces.area.sum():.4%}")

all valid: True | smallest piece: 5.1 m² | empty: 0 | area kept: 99.9993%

Sliver repair on Oakham's parcels and wards

The data used on this page

HM Land Registry's INSPIRE parcels for Rutland (© Crown copyright and database rights; OGL) — the 5,006 that touch the three Oakham wards — and those wards from the ONS December 2024 full-resolution boundaries (OGL), both in British National Grid.

import warnings
from pathlib import Path

import geopandas as gpd
import requests
import shapely

Path("data").mkdir(exist_ok=True)
if not Path("data/Rutland.zip").exists():
    Path("data/Rutland.zip").write_bytes(requests.get(
        "https://use-land-property-data.service.gov.uk/datasets/inspire/download/Rutland.zip", timeout=300).content)
if not Path("data/wards.geojson").exists():
    Path("data/wards.geojson").write_bytes(requests.get(
        "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
        "Wards_December_2024_Boundaries_UK_BFC/FeatureServer/0/query",
        params={"where": "LAD24CD='E06000017'", "outFields": "WD24NM", "outSR": 27700, "f": "geojson"},
        timeout=120).content)
wards = gpd.read_file("data/wards.geojson")
wards = wards[wards["WD24NM"].str.startswith("Oakham")].reset_index(drop=True)
all_parcels = gpd.read_file("/vsizip/data/Rutland.zip/Land_Registry_Cadastral_Parcels.gml",
                            columns=["INSPIREID"], DOWNLOAD_SCHEMA="NO")
parcels = all_parcels[all_parcels.intersects(wards.union_all())].reset_index(drop=True)
warnings.filterwarnings("ignore", message=".*keep_geom_type.*")
print(len(parcels), "parcels touching", wards["WD24NM"].tolist())

5006 parcels touching ['Oakham North East', 'Oakham North West', 'Oakham South']

Example 1: A reusable sliver detector

def flag_slivers(gdf, min_area=5.0, min_thinness=0.01, max_thin_area=100.0):
    """True where a feature is small, or thin and not large enough to be a real long parcel."""
    area = gdf.geometry.area
    thinness = area / gdf.geometry.length ** 2
    return (area < min_area) | ((thinness < min_thinness) & (area < max_thin_area))

flags = flag_slivers(pieces)
print(flags.value_counts().to_dict())

{False: 4960, True: 187}

The max_thin_area cap is what keeps road parcels out of the sliver list; set it from your own data after looking at what the thinness test catches.

work = pieces.reset_index(drop=True)
is_sliver = work.area < 5
keep = work[~is_sliver].copy()
absorbed = orphaned = 0
for _, sliver in work[is_sliver].iterrows():
    candidates = keep[keep.geometry.intersects(sliver.geometry.buffer(0.01))]
    if candidates.empty:
        orphaned += 1
        continue
    shared = candidates.geometry.intersection(sliver.geometry.buffer(0.01)).length
    target = shared.idxmax()
    keep.loc[target, "geometry"] = keep.loc[target, "geometry"].union(sliver.geometry)
    absorbed += 1
print(f"{absorbed} slivers absorbed, {orphaned} with no neighbour; area kept {keep.area.sum() / work.area.sum():.4%}")

159 slivers absorbed, 5 with no neighbour; area kept 99.9997%

Absorbing rather than deleting keeps the total area. The small buffer matters: slivers usually meet their neighbour at a shared edge, and touches on floating-point geometry is unreliable, so this looks for neighbours within a centimetre.

Example 3: Dissolve a clean coverage with `coverage_union_all`

print("wards form a valid coverage:", bool(shapely.coverage_is_valid(wards.geometry.values)))
dissolved = shapely.coverage_union_all(wards.geometry.values)
print(dissolved.geom_type, f"{dissolved.area / 1e6:.2f} km², same as union_all: "
      f"{abs(dissolved.area - wards.union_all().area) < 1}")
try:
    shapely.coverage_union(shapely.GeometryCollection(wards.geometry.tolist()))
except TypeError as exc:
    print("coverage_union with one collection:", exc)

wards form a valid coverage: True
Polygon 9.10 km², same as union_all: True
coverage_union with one collection: coverage_union() missing 1 required positional argument: 'b'

coverage_union_all takes an array of polygons and assumes they already form a clean coverage — edge-matched, non-overlapping. coverage_union is the two-argument form, for a pair of geometries.

Why boundaries that should meet do not

Why near-misses happen. Some are floating-point noise: reproject, round-trip through a file format, or run a constructive operation, and a vertex that should be identical lands at 529384.1200000001 instead of 529384.12. Those are what a millimetre grid fixes. Many more, in real work, are two organisations drawing the same line differently: the parcels and wards on this page disagree by metres, because they record different things. Noise is removed by alignment; disagreement is resolved by deciding which layer is right.

The thinness ratio is what separates artefacts from genuinely small features.

Snapping vs precision reduction. These solve overlapping but different problems:

shapely.snap(geom, reference, tolerance) moves vertices of geom toward a specific reference geometry when they fall within tolerance. It is directional and relational — use it to make one layer conform to another. It does not enforce a global grid, and it can leave self-intersections, so re-validate.
shapely.set_precision(geom, grid_size) rounds every coordinate to a fixed grid with no reference geometry. It is global and absolute — use it to make coincident edges across a whole dataset bit-for-bit identical. It can collapse parts smaller than the grid, which is why it returns valid output by default.

The thinness ratio. For any polygon, $\text{area}/\text{perimeter}^2$ is a scale-invariant shape descriptor. A circle maximises it at 1 / (4*pi) ≈ 0.0796; a square sits at 0.0625. A sliver drives it toward zero — and so does a road parcel, which is why thinness finds candidates rather than deciding.

Tolerances that merge real features

Snapping is a blunt instrument; these four keep it from welding real features together.

You must be in a projected CRS for metric tolerances

grid_size, tolerance, area thresholds, and buffer distances are all in the units of your CRS. In EPSG:4326 those units are degrees, so a tolerance=0.5 means half a degree (tens of kilometres). Always to_crs() to a projected, metric system (a local UTM zone or a national grid; EPSG:3857 metres are stretched by 1/cos(latitude)) before applying any distance- or area-based fix.

A tolerance that is too large merges real features

Snapping and gap-closing are blunt instruments. If your tolerance or grid_size exceeds the size of genuine small features, you will weld neighbours together, delete real narrow polygons, or pull boundaries onto the wrong edge. Start small (millimetres to centimetres), inspect the result, and increase only if needed.

Snapping can create invalid geometry

Both shapely.snap and set_precision can produce self-intersections or zero-area rings as vertices move. Always follow them with make_valid() (GeoPandas: gdf.geometry.make_valid()) and check is_valid.all() before continuing. Skipping this step causes failures further down the pipeline that are hard to trace back.

set_precision can collapse tiny parts

If a polygon part is smaller than grid_size, precision reduction may shrink it to nothing. With the default mode="valid_output" these collapsed parts disappear; with mode="keep_collapsed" they survive as lower-dimensional geometry (a line or point), which can pollute a polygon layer. Choose deliberately and filter out non-polygon types afterwards if needed.

Order of operations matters

Reproject first, then reduce precision / snap, then make_valid, then filter slivers, and only then dissolve or run coverage_union. Filtering slivers before aligning vertices removes artefacts you could have fixed; dissolving before aligning leaves internal boundaries that won't merge. Get the sequence right and each step does its job cleanly.

coverage_union_all needs a valid coverage

shapely.coverage_union_all is fast because it assumes a valid polygonal coverage: features that are edge-matched, non-overlapping, and gap-free. Feed it overlapping or misaligned polygons and you get wrong or undefined results. Align with set_precision first; if you are unsure the input is a clean coverage, validate it (shapely.coverage_is_valid) or fall back to union_all().

Always verify topology after

After any cleanup, re-check: no empty geometries, all valid, plausible total area, and the expected feature count. A silent precision change that deletes 200 real parcels looks identical to a successful run until someone audits the numbers. Make verification a fixed last step, not an afterthought.

Internal links

FAQ

What grid_size should I use for set_precision?

Pick a value smaller than the real-world precision of your data but large enough to absorb floating-point noise. For cadastral or survey data in metres, 0.001 (1 mm) is a safe default. It will not reconcile two layers that disagree by metres — Land Registry parcels and ONS wards kept most of their slivers at 1 mm — and going coarser to force it distorts genuine detail.

What is the difference between snap and set_precision?

shapely.snap aligns one geometry's vertices to a reference geometry within a tolerance — it is relational and directional. shapely.set_precision rounds all coordinates to a fixed grid with no reference — it is global and absolute. Use snap to conform one layer to another; use set_precision to make a whole dataset's coincident edges identical.

Why do I get slivers after gpd.overlay?

Overlay computes geometric intersections, and wherever two boundaries that should coincide differ by floating-point noise, the operation produces a thin polygon along that edge. Reduce precision on both inputs with set_precision before the overlay to make shared edges identical, then filter any remaining slivers by area.

How do I remove slivers without deleting real small features?

Combine an area threshold with a thinness ratio ( $\text{area}/\text{perimeter}^2$ ). Real small features tend to be compact (higher ratio), while slivers are long and thin (ratio near zero). Filtering on both keeps genuine small parcels while dropping artefacts. To preserve area, absorb slivers into their largest neighbour instead of deleting them.

Do I have to reproject before snapping?

Yes, if you want tolerances in metres. In a geographic CRS like EPSG:4326, all distances are in degrees, so a tolerance meant to be a few centimetres becomes kilometres. Always to_crs() to a projected metric CRS first, do the cleanup, then reproject back if needed.

My geometries became invalid after snapping — what now?

That is expected: moving vertices can create self-intersections. Run gdf.geometry.make_valid() immediately after snap or set_precision, then confirm with gdf.is_valid.all(). If make_valid returns mixed geometry types, explode and keep only the polygon parts.

When should I use coverage_union instead of union_all?

Use shapely.coverage_union_all when your polygons already form a clean coverage — edge-matched, non-overlapping, gap-free — because it is much faster and preserves shared internal edges. Use union_all() (or dissolve()) for arbitrary, possibly overlapping geometries. Align with set_precision before either to ensure shared edges actually match.

Can I close gaps with a buffer trick?

Yes — buffering outward by a small distance and then inward by the same distance (a morphological closing) bridges sub-tolerance gaps, but only when applied to the merged geometry: closing each feature separately restores each one to itself and leaves the gaps between them. Use join_style="mitre" to keep corners sharp, and keep the distance small so you do not distort real geometry. For coincident boundaries, set_precision is usually cleaner because it does not change the shape of the polygons themselves.