SpatiaLite & GeoPackage Automation with Python

Snapping Tolerances for Topology Cleanup

Measure two things before picking a number: the largest gap you intend to close, and the smallest genuine separation you must preserve. The tolerance goes…

Measure two things before picking a number: the largest gap you intend to close, and the smallest genuine separation you must preserve. The tolerance goes between them. If those two measurements overlap — if the slivers are as wide as the real gaps — no single tolerance can be correct, and the cleanup needs a different approach entirely.

This page belongs to the Geometry Validity & Topology Repair guide. Where Detecting and Removing Duplicate Vertices is about redundancy within one geometry, this is about gaps between geometries.

Why This Matters

Adjacent parcels that were digitised independently almost never share vertices exactly. The result is a thin sliver of empty space along every shared boundary, a few centimetres wide, invisible at map scale and consequential everywhere else: an area total that omits the slivers, a point-in-polygon test that finds nothing for a point that fell in one, a dissolve that produces holes along every internal boundary.

Snapping closes them by moving nearby vertices onto common positions. The risk is symmetrical: a tolerance large enough to close a five-centimetre sliver is also large enough to merge two features genuinely five centimetres apart, and there is nothing in the geometry to distinguish the two cases. Choosing the number is therefore a data question rather than a technical one, and the technical part — the actual snap — is the easy half.

Prerequisites

  • Python 3.9+ with shapely 2.0+
  • A projected reference system, so the tolerance is in metres rather than degrees
  • A backup: snapping moves coordinates and is not reversible
  • Understanding of what validity does and does not cover, from the parent guide

Primary Method

python
# Measure the two bounds that determine a usable tolerance
from shapely import from_wkb
from shapely.strtree import STRtree


def gap_distribution(geoms, max_distance: float = 1.0) -> list[float]:
    """
    Distances between each geometry and its nearest neighbour, up to
    max_distance. Zero means they already touch; small positive values
    are candidate slivers; larger ones are genuine separations.
    """
    tree = STRtree(geoms)
    gaps = []
    for i, g in enumerate(geoms):
        # everything within max_distance, excluding the geometry itself
        for j in tree.query(g.buffer(max_distance)):
            if j == i:
                continue
            d = g.distance(geoms[j])
            if 0 < d <= max_distance:
                gaps.append(d)
    return sorted(gaps)

The distribution this produces is the whole decision. A layer with a clean answer shows two clusters: one near zero (the slivers) and one at whatever the real minimum separation is, with an empty band between them. The tolerance belongs in that band, and its width is a measure of how safe the operation is.

The usable tolerance band between slivers and real separationsA one-dimensional scale of nearest-neighbour distances. A cluster of sliver gaps sits near zero, up to about four centimetres. A gap follows with no observed distances. A cluster of genuine separations begins at about forty centimetres. Any tolerance inside the empty band closes every sliver and merges nothing real; a tolerance above the band starts merging genuine neighbours.Nearest-neighbour distance, one layerslivers0 – 4 cmthe usable band4 – 40 cm · nothing observed heregenuine separations40 cm and abovepick from the middle of the banda wide band means the operation is safea band of zero width means no tolerance is correct
The band's width, not the tolerance itself, is what tells you whether to run the cleanup at all.

Step-by-Step Walkthrough

1. Measure the gap distribution

python
import sqlite3
from shapely import from_wkb

conn = sqlite3.connect("file:parcels.gpkg?mode=ro", uri=True)
geoms = [from_wkb(gpb_to_wkb(b))
         for (b,) in conn.execute("SELECT geom FROM parcels WHERE geom IS NOT NULL")]

gaps = gap_distribution(geoms, max_distance=2.0)
print(f"{len(gaps)} neighbour pairs within 2 m")
for pct in (50, 75, 90, 95, 99):
    print(f"p{pct}: {gaps[len(gaps) * pct // 100]:.4f} m")

2. Find the band explicitly

Rather than eyeballing percentiles, look for the widest run of distances with nothing in it. That run is the band, and its bounds are the numbers to reason about.

python
def widest_empty_band(gaps: list[float]) -> tuple[float, float]:
    """The largest interval between consecutive observed distances."""
    best = (0.0, 0.0)
    for a, b in zip(gaps, gaps[1:]):
        if (b - a) > (best[1] - best[0]):
            best = (a, b)
    return best


low, high = widest_empty_band(gaps)
print(f"usable band: {low:.4f} m .. {high:.4f} m (width {high - low:.4f} m)")
print(f"suggested tolerance: {(low + high) / 2:.4f} m")

3. Sanity-check the band against the domain

A statistical band is a starting point, not an answer. Two questions decide whether to trust it. Does the lower bound match the capture accuracy — if the slivers are wider than the equipment’s error, they may be real features rather than artefacts. Does the upper bound match something meaningful on the ground — a forty-centimetre minimum separation between parcels is plausible for a boundary wall; a four-metre one probably means the layer contains two different feature classes.

4. Snap, then verify

python
# Snap a layer to the chosen tolerance, refusing anything that breaks
from shapely import set_precision, to_wkb

TOLERANCE = 0.10   # metres, from the middle of the band

conn = sqlite3.connect("parcels.gpkg", isolation_level=None)

updates, skipped = [], []
for fid, blob in conn.execute("SELECT fid, geom FROM parcels WHERE geom IS NOT NULL"):
    geom = from_wkb(gpb_to_wkb(blob))
    snapped = set_precision(geom, TOLERANCE)
    if snapped.is_empty or not snapped.is_valid:
        skipped.append(fid)
        continue
    if abs(snapped.area - geom.area) / max(geom.area, 1e-12) > 0.01:
        skipped.append(fid)          # more than 1% of this feature's area moved
        continue
    updates.append((to_wkb(snapped), fid))

conn.execute("BEGIN IMMEDIATE")
conn.executemany(
    "UPDATE parcels SET geom = GeomFromWKB(?, 27700) WHERE fid = ?", updates
)
conn.execute("COMMIT")
print(f"snapped {len(updates)}; skipped {len(skipped)} for review")

The per-feature area guard is what a global area check cannot give you. Across a layer, slivers closing in one place and opening in another largely cancel, so the total barely moves while individual features change substantially. Checking per feature catches exactly the rows a reviewer needs to see.

Why a per-feature area guard catches what a total cannotTwo views of the same snap. Measured across the whole layer, area barely moves, because slivers absorbed on one boundary are offset by slivers absorbed on another — so a total-area check reports success. Measured per feature, a small number of features changed by several per cent, which is the population a reviewer needs to see. Only the second view distinguishes a good cleanup from a damaging one.layer totalmoves by 0.02%gains and losses cancelreports successand tells you nothingper featuremost move under 0.1%a handful move several %those are the review listthis is the useful checkan aggregate check on a symmetrical operation is close to no check at all
The same reasoning applies to any cleanup whose errors cancel — aggregate agreement is not evidence of per-row correctness.
Order of operations around a snapValidity is repaired first, because snapping an invalid geometry gives the engine an ill-defined input. The snap runs second. Validity is checked again third, because snapping can create new self-intersections where two segments land on the same grid line. The spatial index is rebuilt last, since every snapped geometry has a new bounding box.1 · repair validityclean input to the snap2 · snaptolerance from the band3 · re-check validitythe snap can break it4 · rebuildevery bbox movedrunning steps 1 and 2 the other way round makes each pass undo part of the other
Step three is the one people leave out, and it is the one that catches the damage a slightly-too-coarse tolerance does.

Verification

Re-measure the gap distribution after the pass. Every distance below the tolerance should be gone, and nothing above it should have moved.

python
after = gap_distribution(reloaded_geoms, max_distance=2.0)
below = [g for g in after if g < TOLERANCE]
assert not below, f"{len(below)} gaps remain below the tolerance"

# and the upper cluster should be intact
assert min(g for g in after) >= TOLERANCE
sql
-- No feature may have become invalid, and the count must be unchanged
SELECT count(*) AS invalid FROM parcels WHERE ST_IsValid(geom) = 0;
SELECT count(*) AS features FROM parcels;

A dropped feature count is the failure to watch for: a feature that collapsed entirely under the snap disappears, and unlike an invalid geometry, nothing else reports it.

Alternative Approaches or Edge Cases

When the band has zero width. If sliver gaps and genuine separations overlap in size, no tolerance separates them and snapping will damage the layer whatever value you choose. The honest options are to fix the source — re-digitise shared boundaries against a common topology — or to snap only the pairs you can identify as artefacts by some other signal, such as sharing an attribute value that says they are neighbours.

Snapping to a reference layer. Where an authoritative boundary layer exists, snapping the working layer to it rather than to itself is both safer and more meaningful: the tolerance then answers “how far may a vertex have drifted from the authority”, which is a question with a defensible answer. shapely.snap(geom, reference, tolerance) does this per geometry.

Geographic coordinates. A tolerance in degrees is not a distance, and it varies with latitude. Project to a suitable planar system before measuring or snapping — the alternative is an effective tolerance that shrinks toward the poles, which produces a cleanup that works in the south of an extent and not in the north.

Troubleshooting

Features disappear after snapping

Cause: A feature smaller than the tolerance collapsed to nothing. Fix: Guard on is_empty before writing, as the pass does, and compare feature counts afterwards. Small features — survey markers, structures, isolated fragments — are exactly the ones a tolerance chosen from parcel boundaries will destroy.

Gaps remain after the pass

Cause: The two geometries were further apart than the tolerance, or set_precision snapped them to different grid lines. Grid snapping aligns coordinates to a lattice rather than to each other, so two vertices either side of a grid boundary can stay apart. Fix: Use shapely.snap between the specific pair, or increase the grid size — but re-check the band first, because increasing it is what merges real neighbours.

New self-intersections appear

Cause: Two segments of one ring snapped onto the same line. Fix: Repair with the pass in How to Repair Self-Intersecting Polygons with Shapely, and treat a high rate of new invalidity as evidence the tolerance is too coarse.

Frequently Asked Questions

Should snapping run before or after a validity repair?

Validity repair first, snapping second, validity check again. Snapping an already-invalid geometry gives GEOS an ill-defined input and the result is unpredictable; snapping a valid one can create new invalidity, which is why the check repeats. Running the two in the other order tends to produce a layer where each pass undoes some of the other’s work.

Does this fix gaps between layers as well as within one?

Not as written — the measurement and the snap both operate within a single layer. For cross-layer alignment, snap the working layer to the authoritative one with shapely.snap, which moves the working geometry onto the reference’s vertices and leaves the reference untouched. That asymmetry is the point: one layer is the authority and the other is being corrected toward it.

How does this interact with the spatial index?

Every snapped geometry has a new bounding box, so the index must be rebuilt afterwards exactly as with any geometry-changing update. The measurement phase, by contrast, benefits enormously from an index — the STRtree in the primary method is a Shapely-side equivalent, built once in memory, because doing the nearest-neighbour search without one is quadratic in the feature count.