Snapping and Simplifying Polygons Without Creating Gaps

Run GeoSeries.simplify() over a layer of administrative boundaries and the file gets smaller, the map gets faster, and hairline white cracks appear along every border you thought was shared. This guide is for anyone thinning a polygon layer for a web map, a tile build, or a storage budget without corrupting the areas and adjacencies downstream analysis depends on. It sits under Topology Validation & Repair in Geospatial Data Ingestion & Processing Workflows, and picks up where Fixing Self-Intersecting Polygons Programmatically stops: that guide repairs defects inside one polygon, this one repairs defects between polygons.

Why This Approach / What Goes Wrong

A layer of parcels, census tracts, or flood zones is usually a coverage: the polygons tile a region with no gaps and no overlaps, and any two neighbours store the border between them as the same list of coordinates, once per feature. Douglas-Peucker simplification does not know that. GeoSeries.simplify() hands each geometry to GEOS on its own, so the shared border is fed into the algorithm twice — once as part of the west polygon's ring, once as part of the east polygon's ring — with different neighbouring vertices around it each time. Douglas-Peucker anchors on the first and last point of a segment and recursively keeps the vertex furthest from the chord, so a different surrounding ring produces a different set of survivors. The two copies of the border come out as two different lines.

Where the west copy bows east of the east copy, both features claim the same ground and you get an overlap; where it bows the other way, nothing covers the strip and you get a gap. Neither is an invalid geometry — every individual polygon passes is_valid, which is why this survives the validity gate described on the parent topic and only surfaces later as white seams in a vector tile, double-counted population in a dissolve, or points that fall into no zone at all in a spatial join.

preserve_topology=True, the default, does not help. It preserves the topology of each geometry in isolation: no new self-intersection, no collapsed ring, no hole escaping its shell. It has no view of the other 4,000 rows in the table. The same blind spot applies to buffer(0), to per-feature rounding, and to writing rounded WKT — anything that moves vertices one feature at a time will desynchronise a shared edge.

The fix has three parts. First, make the shared edges exactly identical by snapping every coordinate onto a grid — otherwise the topology engine cannot tell "shared" from "0.4 mm apart". Second, simplify the coverage rather than the features, so each shared edge is thinned once and the result is written back to both neighbours. Third, audit the output with a union-area check that would have caught the problem on the naive path. GEOS 3.12 added a coverage simplifier for exactly this, exposed as shapely.coverage_simplify and, since GeoPandas 1.1, as GeoSeries.simplify_coverage(); the topojson package solves the same problem by construction, storing shared edges as arcs referenced by both polygons.

Prerequisites

• geopandas>=1.1 — adds GeoSeries.simplify_coverage(), is_valid_coverage() and invalid_coverage_edges()
• shapely>=2.1 — exposes coverage_simplify, coverage_is_valid and coverage_invalid_edges; requires GEOS 3.12 or newer
• pyproj>=3.7 — CRS objects and the UTM zone estimator
• topojson>=1.9 — the arc-based fallback when your GEOS is older, or when you want the arcs themselves

python -m pip install "geopandas>=1.1" "shapely>=2.1" "pyproj>=3.7" "topojson>=1.9"

import shapely
import geopandas as gpd

print(shapely.__version__, shapely.geos_version)      # 2.1.2 (3, 13, 1)
print(hasattr(gpd.GeoSeries, "simplify_coverage"))    # True

If shapely.geos_version is below (3, 12, 0) the coverage functions raise GEOSException; skip to the topojson path in step 5. Conda-forge builds keep GEOS consistent across Shapely, GeoPandas and pyproj — mixing pip wheels with a system GEOS is the usual reason a machine reports 3.11 when the changelog promised 3.13.

Step-by-Step Implementation

1. Move into a metric CRS before you name any tolerance. Every number in this guide — grid size, snap distance, simplification tolerance, minimum sliver area — is in the units of the layer's CRS. In EPSG:4326 those units are degrees, where a 0.0001 tolerance means about 11 m of latitude but only 7 m of longitude in Berlin and 11 m at the equator, so the same call thins the north-south edges harder than the east-west ones. Project to the local UTM zone, never to Web Mercator, whose scale factor grows with latitude and would make a 12 m tolerance mean 20 m on the ground in Scandinavia. The estimator is covered in Choosing a UTM Zone Automatically in Python.

import geopandas as gpd

flood_zones = gpd.read_file("flood_zones.gpkg", layer="zones", engine="pyogrio")
if flood_zones.crs is None:
    raise ValueError("Layer has no CRS — a tolerance in metres is meaningless without one.")

source_crs = flood_zones.crs
zones = flood_zones.to_crs(flood_zones.estimate_utm_crs())   # e.g. EPSG:32633

assert zones.crs.is_projected
print(zones.crs.axis_info[0].unit_name)     # metre

2. Derive the two tolerances from the data and the output, not from a hunch. They answer different questions. The snap grid asks how precisely the source was digitised: it must be fine enough to leave real vertices alone and coarse enough to weld coordinates that differ only in their last floating-point digits. Two orders of magnitude below the finest observed vertex spacing is a safe default. The simplify tolerance asks how the layer will be seen: at the largest zoom you publish, a vertex closer than one screen pixel to its neighbour cannot be drawn distinctly, so that pixel's ground size is the natural ceiling.

import numpy as np
import shapely


def vertex_spacing(geoms) -> np.ndarray:
    """Distance between consecutive vertices on every exterior ring, in CRS units."""
    rings = shapely.get_exterior_ring(shapely.get_parts(geoms))
    coords, ring_idx = shapely.get_coordinates(rings, return_index=True)
    steps = np.hypot(*np.diff(coords, axis=0).T)
    return steps[ring_idx[1:] == ring_idx[:-1]]      # drop jumps between rings


def tile_pixel_metres(zoom: int, latitude: float, tile_px: int = 512) -> float:
    """Ground size of one screen pixel in a Web Mercator tile at this zoom."""
    metres_per_px_z0 = 156543.03392804097 * 256 / tile_px
    return metres_per_px_z0 * np.cos(np.radians(latitude)) / (2 ** zoom)


steps = vertex_spacing(zones.geometry.values)
finest = float(np.percentile(steps, 1))
snap_grid = max(10 ** np.floor(np.log10(finest) - 2), 0.001)     # never finer than 1 mm
simplify_tolerance = tile_pixel_metres(zoom=12, latitude=52.5)

print(f"finest spacing : {finest:.2f} m")            # finest spacing : 4.00 m
print(f"snap grid      : {snap_grid} m")             # snap grid      : 0.01 m
print(f"tolerance      : {simplify_tolerance:.1f} m")  # tolerance      : 11.6 m

3. Weld the shared nodes with set_precision. Rounding every coordinate onto the same grid is what turns "two borders that look identical" into "two borders that are identical". This matters more than it sounds: a layer that has been through a reprojection, a segmentize, or a round-trip via Shapefile will typically have shared edges that differ in the eleventh decimal, and GEOS will treat those as two separate lines with a nanometre-wide gap between them. mode="valid_output" (the default) guarantees each rounded geometry stays valid.

zones["geometry"] = shapely.set_precision(
    zones.geometry.values, grid_size=snap_grid, mode="valid_output"
)

When the misalignment is larger than a rounding error — a newly digitised layer that has to meet an authoritative cadastre — use shapely.snap to pull the loose vertices onto the reference linework first, with a tolerance smaller than the narrowest real feature so you do not collapse a lane or a hedgerow.

reference_edges = shapely.union_all(cadastre.geometry.values)
new_parcels["geometry"] = shapely.snap(
    new_parcels.geometry.values, reference_edges, tolerance=0.5    # 0.5 m
)

4. Prove the layer is a coverage before you simplify it. is_valid_coverage() returns a single boolean for the whole GeoSeries; invalid_coverage_edges() returns the offending edge of each row as a LineString (empty where the row is fine), which is what you write to a debug GeoPackage and open in QGIS.

if not zones.geometry.is_valid_coverage():
    bad = zones.geometry.invalid_coverage_edges()
    zones.loc[~bad.is_empty].to_file("coverage_errors.gpkg", driver="GPKG")
    raise ValueError(f"{int((~bad.is_empty).sum())} rows have mismatched shared edges")

5. Simplify the coverage, not the features. simplify_coverage() thins each shared edge once and writes the identical result back into both neighbours, so adjacency and total area survive by construction. simplify_boundary=False pins the outer boundary of the whole coverage and thins only the interior borders — the right choice when the outline is a legal extent (a national border, a study-area clip) that must not move.

zones["geometry"] = zones.geometry.simplify_coverage(
    simplify_tolerance, simplify_boundary=True
)

On GEOS older than 3.12, or when you want the shared arcs themselves to feed a TopoJSON or vector-tile build, use the topojson package. It decomposes the layer into arcs, simplifies each arc once, and reassembles the polygons from the shared result. Pass prequantize=False to keep full coordinate precision (quantisation is itself a snap, and doing it twice is how you lose a metre you never budgeted for) and shared_coords=False so junctions are detected from shared line segments rather than exact coordinate equality.

import topojson as tp

topo = tp.Topology(zones, prequantize=False, shared_coords=False)
zones = topo.toposimplify(simplify_tolerance).to_gdf()   # CRS and columns preserved

6. Repair what is left, and absorb slivers rather than deleting them. A layer you inherited already torn — someone else's naive simplification, or an overlay of two independently-thinned layers — needs repair instead of re-simplification. Run make_valid first (the mechanics are in Fixing Self-Intersecting Polygons Programmatically), then deal with the fragments. Deleting a sliver re-opens the hole it was filling, so merge each one into the neighbour it shares the longest edge with; Polsby-Popper compactness, 4πA / P², separates genuine thin features from artefacts far better than area alone.

import numpy as np
import shapely
import geopandas as gpd


def absorb_slivers(gdf, min_area, max_thinness=0.05):
    """Merge sliver polygons into the neighbour they share the longest edge with."""
    parts = gdf.explode(index_parts=False, ignore_index=True)
    parts = parts[parts.geom_type == "Polygon"].copy()

    # Polsby-Popper: 1.0 is a circle, a hairline sliver approaches 0
    thinness = 4 * np.pi * parts.area / parts.length ** 2
    is_sliver = (parts.area < min_area) & (thinness < max_thinness)
    keep, slivers = parts[~is_sliver].copy(), parts[is_sliver]
    if slivers.empty:
        return keep

    pairs = slivers[["geometry"]].sjoin(keep[["geometry"]], predicate="intersects")
    best = {}
    for sliver_idx, host_idx in zip(pairs.index, pairs["index_right"]):
        shared = slivers.geometry.loc[sliver_idx].intersection(
            keep.geometry.loc[host_idx]).length
        if shared > best.get(sliver_idx, (None, -1.0))[1]:
            best[sliver_idx] = (host_idx, shared)

    for sliver_idx, (host_idx, _) in best.items():
        keep.loc[host_idx, "geometry"] = shapely.union(
            keep.geometry.loc[host_idx], slivers.geometry.loc[sliver_idx]
        )
    return keep


zones["geometry"] = shapely.make_valid(zones.geometry.values)
zones = absorb_slivers(zones, min_area=10 * simplify_tolerance ** 2)
zones = zones.to_crs(source_crs)      # back to the caller's CRS, last step only

Verification

The union-area check is the audit that catches every failure mode at once. Summing the individual areas and subtracting the area of the union isolates double-counted ground, so a non-zero result means overlaps. Counting interior rings in the union finds enclosed gaps. Comparing the union area before and after finds net drift, including gaps that open onto the outer boundary and therefore leave no ring behind. Run all three against both approaches on the same input and the difference is unambiguous.

import numpy as np
import shapely
import geopandas as gpd
from shapely.geometry import Polygon

# A three-zone coverage sharing two wavy borders, in UTM 33N (metres)
y = np.arange(0, 240.5, 4.0)
west_edge = list(zip(8 * np.sin(y / 15 + 0.3), y))
east_edge = list(zip(120 + 8 * np.sin(y / 21 + 0.7), y))
zones = gpd.GeoDataFrame(
    {"zone_id": ["W", "M", "E"]},
    geometry=[
        Polygon([(-80, 0)] + west_edge + [(-80, 240)]),
        Polygon(west_edge + east_edge[::-1]),
        Polygon(east_edge + [(220, 240), (220, 0)]),
    ],
    crs="EPSG:32633",
)
zones["geometry"] = shapely.set_precision(zones.geometry.values, grid_size=0.01)


def coverage_report(before: gpd.GeoSeries, after: gpd.GeoSeries) -> dict:
    """Union-area audit: overlaps, enclosed gaps, and net area drift."""
    union_before, union_after = before.union_all(), after.union_all()
    return {
        "gaps": int(shapely.get_num_interior_rings(shapely.get_parts(union_after)).sum()),
        "overlap_m2": round(float(after.area.sum() - union_after.area), 2),
        "drift_m2": round(float(union_after.area - union_before.area), 2),
        "is_coverage": bool(after.is_valid_coverage()),
        "vertices": int(shapely.get_num_coordinates(after.values).sum()),
    }


tolerance = 11.6      # one z12 tile pixel at latitude 52.5
print("naive   :", coverage_report(zones.geometry, zones.geometry.simplify(tolerance)))
print("coverage:", coverage_report(zones.geometry, zones.geometry.simplify_coverage(tolerance)))

clean = zones.geometry.simplify_coverage(tolerance)
assert clean.is_valid_coverage(), "shared borders no longer match"
assert abs(clean.area.sum() - clean.union_all().area) < 1e-6, "zones overlap"
assert abs(clean.union_all().area - zones.geometry.union_all().area) < 1e-6, "area drifted"

Expected console output:

naive   : {'gaps': 1, 'overlap_m2': 424.79, 'drift_m2': -1261.65, 'is_coverage': False, 'vertices': 19}
coverage: {'gaps': 0, 'overlap_m2': 0.0, 'drift_m2': 0.0, 'is_coverage': True, 'vertices': 33}

The naive path is smaller — 19 vertices against 33 — and that is exactly the trap: it wins on file size by deleting the agreement between neighbours. It lost 1,261 m² of the study area outright and double-counted another 425 m². The coverage path keeps a few more vertices and every one of them is load-bearing. Sweeping the tolerance shows where the useful range ends: below the display resolution you pay displacement for vertices no one will see, and past it you buy almost no further savings.

Edge Cases & Debugging

• GEOSException from simplify_coverage. Your GEOS predates 3.12. Check shapely.geos_version, upgrade from conda-forge, or use the topojson path in step 5 — it has no GEOS version floor.
• is_valid_coverage() is False before you have touched anything. The source already has mismatched edges. Try a coarser set_precision grid first; if rows still fail, pass gap_width (for example is_valid_coverage(gap_width=0.05)) to report near-misses narrower than 5 cm and inspect them with invalid_coverage_edges().
• Overlapping polygons are not a coverage at all. Land-use layers where a wetland sits inside a forest polygon will never pass the gate. Flatten them into disjoint pieces with overlay(how="union") — see Computing Overlay Union & Difference with GeoPandas — or simplify each thematic layer separately.
• Simplification collapses a whole polygon. Features narrower than the tolerance can vanish, taking their attribute row with them. Compare len() before and after, and clamp the tolerance to a fraction of the smallest feature's width when the row count must be preserved.
• The layer is fine but the tiles still show cracks. The tile builder re-quantises coordinates into integer tile space. Simplify at or above the resolution of your highest zoom and let the builder round, rather than snapping to your own grid and letting it round again — the pattern used in Generating PMTiles from GeoParquet.
• Everything drifts after the final to_crs. Reprojection moves every vertex, and a coverage that was exact in UTM can pick up nanometre mismatches in EPSG:4326. Re-run is_valid_coverage() after the reprojection, and re-snap with a set_precision grid appropriate to degrees (1e-9 is roughly 0.1 mm) if it fails.

Frequently Asked Questions

Does preserve_topology=True prevent gaps between polygons? No. It only guarantees that the simplified version of one geometry stays valid on its own — no self-intersection, no ring collapse, no hole escaping its shell. GEOS evaluates it per geometry with no knowledge of the neighbouring rows, so a shared border is still processed twice and can still come out as two different lines. That distinction is the single most common misreading of the API, and it is why the union-area check in the verification section belongs in your pipeline rather than a trust in the flag.

Can I simplify in EPSG:4326 and skip the reprojection? You can, but the tolerance stops meaning one thing. A degree of latitude is about 111 km everywhere while a degree of longitude shrinks with the cosine of latitude, so a single numeric tolerance thins east-west edges harder than north-south ones, and the effect changes across a country-sized layer. If the layer is small and near the equator the distortion may be tolerable; otherwise project to the local UTM zone with estimate_utm_crs(), do the work, and reproject back at export. The projection mechanics are covered in Coordinate Systems with PyProj.

When should I use topojson instead of simplify_coverage? Use simplify_coverage when your GEOS is 3.12 or newer and you want the fastest path with the fewest moving parts — it is compiled C, vectorized over the whole array. Reach for topojson when GEOS is older, when you need the arc structure itself (feeding a TopoJSON file or a tile pipeline that consumes shared arcs), or when you want quantisation and simplification controlled together. Both produce a topologically consistent result; only the second keeps the shared-arc representation after the call returns.

How do I choose the minimum area for the sliver filter? Anchor it to the tolerance rather than to a round number of square metres: an artefact created by a tolerance t is rarely wider than t, so min_area around 10 * t² catches the debris without touching real features. Pair it with a compactness test — Polsby-Popper below about 0.05 — so that legitimately long, thin polygons such as rivers and rail corridors survive. Then absorb rather than delete, because deleting a sliver reopens the gap it was filling, and re-run the coverage gate afterwards to prove the merge closed cleanly.