SpatiaLite & GeoPackage Automation with Python

Detecting and Removing Duplicate Vertices

Compare consecutive coordinate pairs against a tolerance, not for equality — exactly coincident vertices are rare and near-coincident ones are everywhere.…

Compare consecutive coordinate pairs against a tolerance, not for equality — exactly coincident vertices are rare and near-coincident ones are everywhere. Remove them with shapely.set_precision(geom, grid_size), which snaps coordinates to a grid and collapses the resulting duplicates, then re-check validity because collapsing a vertex can change ring topology.

This page belongs to the Geometry Validity & Topology Repair guide. Duplicate vertices are the clearest example of a defect that validity checking cannot see.

Why This Matters

Two vertices a nanometre apart are perfectly valid. The ring does not cross itself, the holes are properly contained, and ST_IsValid returns 1 without hesitation. What those vertices do instead is generate a segment of essentially zero length, and zero-length segments are where geometry engines produce results nobody can explain: a spatial join that matches on one machine and not another, an overlay that produces slivers, a simplification that removes the wrong vertex because the one it should have removed had no length to measure.

They also cost real storage. A polygon captured by a GPS logger at one-second intervals while the surveyor stood still can carry a hundred vertices covering half a metre. Across a layer that is megabytes of coordinates describing nothing, and every predicate that touches those geometries pays for all of them.

Prerequisites

  • Python 3.9+ with shapely 2.0+ — set_precision arrived in 2.0
  • GEOS 3.9+ for the precision model this relies on
  • A clear idea of the data’s real precision, in the layer’s own units
  • A backup: collapsing vertices is not reversible

Primary Method

python
# Find near-duplicate consecutive vertices and report the worst offenders
from shapely import from_wkb, set_precision
from shapely.geometry import Polygon, MultiPolygon
import math


def duplicate_runs(geom, tolerance: float) -> int:
    """Count consecutive vertex pairs closer together than the tolerance."""
    parts = geom.geoms if isinstance(geom, MultiPolygon) else [geom]
    count = 0
    for part in parts:
        for ring in [part.exterior, *part.interiors]:
            coords = list(ring.coords)
            for (x1, y1), (x2, y2) in zip(coords, coords[1:]):
                if math.hypot(x2 - x1, y2 - y1) < tolerance:
                    count += 1
    return count


def deduplicate(geom, grid_size: float):
    """
    Snap coordinates to a grid, which collapses vertices closer than grid_size.
    Returns (new_geom, vertices_removed). new_geom is None if the result is
    empty or invalid — a sign the tolerance was too coarse for this feature.
    """
    before = sum(len(r.coords) for p in (geom.geoms if isinstance(geom, MultiPolygon) else [geom])
                 for r in [p.exterior, *p.interiors])
    snapped = set_precision(geom, grid_size)
    if snapped.is_empty or not snapped.is_valid:
        return None, 0
    after = sum(len(r.coords) for p in (snapped.geoms if isinstance(snapped, MultiPolygon) else [snapped])
                for r in [p.exterior, *p.interiors])
    return snapped, before - after

set_precision is the right tool because it is topology-aware: it snaps coordinates and then repairs the consequences, rather than deleting vertices and leaving whatever results. Removing coordinates directly from the array is faster and occasionally produces a self-intersection where two segments crossed after one of them lost its middle.

What near-duplicate vertices cost, and what notices themFour consequences of near-duplicate vertices. Validity checking does not notice them at all. Storage grows in proportion to the redundant coordinates. Predicate cost grows with vertex count, so every spatial test pays. Overlay operations can produce slivers where a zero-length segment meets another geometry. Only the last is usually noticed, and it is noticed as a mysterious result rather than as a vertex problem.One redundant vertex, four effectsvalidity checkingno signal whatsoevera duplicate vertex is legalstorage16 bytes per redundant vertexmegabytes across a layerpredicate costscales with vertex countevery spatial test paysoverlay resultsslivers from zero-length segmentsreported as an unexplainable result
The pink box is the only one anyone reports, and it is reported as something else entirely.

Step-by-Step Walkthrough

1. Measure before choosing a tolerance

A tolerance chosen without looking at the data will either do nothing or destroy detail. Measure the distribution of segment lengths first.

python
import sqlite3, math
from shapely import from_wkb

conn = sqlite3.connect("file:field.gpkg?mode=ro", uri=True)

lengths = []
for (blob,) in conn.execute("SELECT geom FROM parcels WHERE geom IS NOT NULL"):
    geom = from_wkb(gpb_to_wkb(blob))
    parts = geom.geoms if geom.geom_type == "MultiPolygon" else [geom]
    for part in parts:
        for ring in [part.exterior, *part.interiors]:
            cs = list(ring.coords)
            lengths += [math.hypot(b[0] - a[0], b[1] - a[1]) for a, b in zip(cs, cs[1:])]

lengths.sort()
for pct in (1, 5, 25, 50):
    print(f"p{pct:>2}: {lengths[len(lengths) * pct // 100]:.6f}")

A healthy layer shows a smooth distribution starting near the capture resolution. A layer with duplicate vertices shows a cluster at or near zero, well separated from the rest — and the gap between that cluster and the real segments is where the tolerance belongs.

2. Choose the tolerance from the data’s real precision

Two anchors are useful. The first is the capture accuracy: a handheld GPS gives a few metres, so vertices under a decimetre apart carry no information. The second is the coordinate unit: in a projected system the units are metres and the numbers are intuitive; in a geographic system they are degrees, and a tolerance of 1e-7 is roughly a centimetre — a mistake of three orders of magnitude here is easy and expensive.

3. Deduplicate one feature and inspect

python
from shapely import from_wkb

geom = from_wkb(gpb_to_wkb(blob))
snapped, removed = deduplicate(geom, grid_size=0.05)   # 5 cm, projected units

print(f"removed {removed} vertices")
print(f"area {geom.area:.6f} -> {snapped.area:.6f}")
print(f"max displacement <= {0.05 * 1.415:.4f} units")

The displacement bound is worth stating explicitly: snapping to a grid of size g can move any vertex by at most half a grid cell diagonally, so no coordinate moves further than g × √2 / 2. That is the guarantee to quote when someone asks whether the operation changed the data.

4. Run the pass with a validity guard

python
# Deduplicate a layer, refusing any feature the tolerance would break
conn = sqlite3.connect("field.gpkg", isolation_level=None)

updates, skipped, removed_total = [], [], 0
for fid, blob in conn.execute("SELECT fid, geom FROM parcels WHERE geom IS NOT NULL"):
    geom = from_wkb(gpb_to_wkb(blob))
    snapped, removed = deduplicate(geom, grid_size=0.05)
    if snapped is None:
        skipped.append(fid)
        continue
    if removed:
        updates.append((to_wkb(snapped), fid))
        removed_total += removed

conn.execute("BEGIN IMMEDIATE")
conn.executemany(
    "UPDATE parcels SET geom = GeomFromWKB(?, 4326) WHERE fid = ?", updates
)
conn.execute("COMMIT")
print(f"{removed_total} vertices removed across {len(updates)} features; "
      f"{len(skipped)} skipped as too small for the tolerance")

The skip list matters. A feature whose whole extent is smaller than the grid collapses to nothing, and those are exactly the features — a survey marker, a small structure — that a coarse tolerance would silently delete.

Choosing a snapping tolerance against feature sizeThree tolerance choices against the same layer. A tolerance far below the capture accuracy removes almost nothing and leaves the redundant vertices in place. A tolerance matched to the capture accuracy removes the redundant vertices while moving no coordinate further than one grid diagonal. A tolerance approaching the size of the smallest features collapses them entirely, which is why the pass must skip and report rather than write the result.too finebelow capture accuracyremoves almost nothingsafe and pointlessthe redundancy remainsmatchedat the capture accuracyremoves the redundancybounded displacementwhat you wanttoo coarsenear smallest feature sizecollapses small featuressilently, if unguardedskip and report insteadthe usable band sits between capture accuracy and the smallest feature you must keep
The band is usually wide. The failure is choosing a tolerance without measuring either edge of it.

Verification

python
# Nothing should remain below the tolerance, and area should barely move
remaining = 0
for (blob,) in conn.execute("SELECT geom FROM parcels WHERE geom IS NOT NULL"):
    remaining += duplicate_runs(from_wkb(gpb_to_wkb(blob)), tolerance=0.05)
assert remaining == 0, f"{remaining} near-duplicate pairs remain"
sql
-- Area should move by far less than 1% of the total
SELECT round(sum(ST_Area(geom)), 6) AS total_area FROM parcels;
SELECT count(*) AS invalid FROM parcels WHERE ST_IsValid(geom) = 0;

The validity count is the one to watch. Snapping can turn a valid geometry invalid when two rings collapse onto one another, and set_precision repairs most but not all of those cases — so a validity scan belongs immediately after the pass, not in next week’s gate.

The displacement guarantee a grid snap gives youSnapping to a grid of a given size moves any coordinate by at most half a grid cell in each axis, so the largest possible displacement is half the grid diagonal. For a five centimetre grid that bound is about three and a half centimetres. Because the bound is a property of the grid rather than of the data, it can be quoted before the pass runs and verified afterwards.grid size gchosen from capture accuracymax displacementg × √2 ÷ 2g = 0.05 mno vertex moves > 3.5 cmthe bound depends on the grid, not the data — quote it before the pass runs
This is the sentence to have ready when someone asks whether the cleanup changed the survey.

Alternative Approaches or Edge Cases

Simplification instead of snapping. shapely.simplify with a small tolerance removes vertices that lie close to the line between their neighbours, which is a different goal: it thins detail rather than removing redundancy. Where the aim is to shrink a layer for rendering, simplification is the right tool; where the aim is to remove coordinates that carry no information, snapping is.

Doing it at capture time. A logger that records a vertex per second while stationary is the source of the problem, and a capture-side minimum-displacement filter removes it before it reaches the container. That is far cheaper than a repair pass and does not require choosing a tolerance retrospectively.

Layers in geographic coordinates. A grid size in degrees is not uniform on the ground — a degree of longitude shrinks toward the poles — so a single tolerance over a wide extent snaps more aggressively at high latitude. For anything beyond a small area, project first, snap, and project back, or accept that the effective tolerance varies.

Troubleshooting

set_precision returns an empty geometry

Cause: The whole feature is smaller than the grid cell. Fix: Skip and report the feature, as the guarded pass does. This is not a repair failure; it is the tolerance telling you it is too coarse for part of the data.

Vertex count barely drops

Cause: The tolerance is below the actual spacing of the redundant vertices. Fix: Re-measure the segment-length distribution and pick a value above the near-zero cluster. A GPS logger’s “stationary” vertices are typically metres apart, not millimetres, because of positional noise.

Geometries become invalid after snapping

Cause: Two rings, or two parts of one ring, snapped onto the same grid line. Fix: Run the validity repair from How to Repair Self-Intersecting Polygons with Shapely afterwards, or use a finer grid. A tolerance that produces invalid output on more than a handful of features is too coarse.

Frequently Asked Questions

How much storage does this actually recover?

It depends entirely on the capture method, and the range is wide. A layer digitised from imagery by a person rarely has meaningful redundancy. A layer walked with a GPS logger can be half redundant vertices, because a stationary receiver keeps recording. Since each coordinate pair is sixteen bytes, halving the vertex count of a large polygon layer is a substantial reduction — and the predicate speed-up that comes with it is usually the bigger win.

Should the tolerance be stored with the data?

Yes, alongside whatever else the pipeline records about a container. A consumer that knows coordinates were snapped to five centimetres can interpret a comparison accordingly; one that does not will eventually treat a five-centimetre difference as a discrepancy. The provenance table pattern in Datum Transformations & Projection Accuracy is a reasonable place to put it.

Does this help the spatial index?

Only indirectly. The R-tree stores bounding boxes, and removing interior vertices does not change a bounding box, so the index itself is unaffected. What improves is everything past the index: the exact predicate applied to the candidates has fewer vertices to consider, which is where the time in an indexed spatial query actually goes.