This guide is part of Automated Vector & Raster Cleaning Workflows.

Geometry Repair with Shapely & GeoPandas

Invalid geometries are the silent killers of spatial ETL pipelines. A single self-intersecting polygon, a collapsed ring, or an incorrectly oriented boundary can cascade into failed spatial joins, corrupted raster extractions, and silent data loss downstream. This guide outlines a production-tested, deterministic methodology for identifying, repairing, and validating geometries at scale using Shapely 2.0’s vectorized API and GeoPandas, ensuring your spatial data pipelines remain resilient across heterogeneous source datasets.

Unlike ad-hoc fixes applied per-project, a standardised repair workflow belongs in the transformation stage of every ingestion pipeline — slotted in after CRS Normalization Across Mixed Datasets and before spatial joins or raster extractions.

Why Naive Approaches Fail

The classic buffer(0) workaround — applying a zero-distance buffer to coerce invalid polygons into valid ones — was the de facto standard before Shapely 2.0. It fails in three predictable ways:

1. Topology destruction: buffer(0) discards bowtie (figure-eight) polygons entirely rather than splitting them into two valid parts. The geometry silently vanishes from your dataset.
2. GEOS version sensitivity: The output of buffer(0) varies across GEOS minor versions, making pipelines non-reproducible when containers or environments are updated.
3. Row-wise .apply() overhead: Pre-2.0 workflows typically called .apply(lambda g: g.buffer(0)) on large GeoDataFrames, iterating in Python rather than using C-accelerated array operations. On datasets above 50k features, this adds minutes of unnecessary latency.

shapely.make_valid() in Shapely 2.0 replaces all of these patterns. It accepts NumPy arrays of geometry objects, delegates to the GEOS topology engine for deterministic reconstruction, and preserves spatial extent rather than discarding problematic rings.

Prerequisites & Environment

Ensure your environment meets these minimum requirements before implementing any repair routine. Legacy Shapely 1.x APIs are incompatible with the vectorized patterns shown here.

pip install "geopandas>=1.0" "shapely>=2.0" pyogrio numpy pandas

Verify GEOS availability — Shapely’s validation and repair functions depend entirely on the underlying GEOS topology engine:

import shapely
print(shapely.geos_version)       # e.g. (3, 12, 1)
print(shapely.geos_version_string) # e.g. '3.12.1-CAPI-1.18.1'

Pin the GEOS version in containerised deployments. Silent behavioural drift between GEOS 3.10 and 3.12 affects make_valid() output geometry types in edge cases involving nested rings.

Version / Compatibility Matrix

For pyogrio chunked I/O (used in the performance section below), version 0.7+ is required for the skip_features / max_features parameters.

Step-by-Step Implementation

Step 1 — Ingest with Schema Preservation

Load the dataset using pyogrio for optimal I/O throughput. Capture the geometry column name before any transformations, as some operations silently rename it.

import geopandas as gpd

gdf = gpd.read_file("source_data.gpkg", engine="pyogrio")
original_geom_col = gdf.geometry.name
print(f"Loaded {len(gdf)} features, geometry column: '{original_geom_col}'")

Avoid implicit type coercion during reads. If your source has mixed-type numeric columns, use read_info() first to inspect field types and apply explicit dtypes overrides where needed.

Step 2 — Diagnose Invalid Topologies

Use gdf.geometry.is_valid to generate a boolean mask. Log counts, indices, and geometry types before touching any data.

import shapely
import numpy as np

invalid_mask = ~gdf.geometry.is_valid
invalid_count = int(invalid_mask.sum())

if invalid_count > 0:
    print(f"Found {invalid_count} invalid geometries ({invalid_count/len(gdf)*100:.1f}%)")
    invalid_types = gdf.loc[invalid_mask, "geometry"].geom_type.value_counts()
    print(invalid_types)

Do not assume all invalid geometries are self-intersections. Common topology violations include:

• Self-intersection (bowtie/figure-eight): polygon rings cross themselves
• Unclosed rings: first and last coordinate differ beyond floating-point tolerance
• Duplicate vertices: consecutive identical coordinates violate simple geometry rules
• Nested shells: inner ring lies outside the outer ring
• Collapsed rings: a polygon degenerates to a line or point after coordinate rounding

For deeper diagnostics, shapely.is_valid_reason() returns human-readable GEOS error strings. This is especially valuable before deciding whether Spatial Deduplication & Topology Simplification should precede the repair step.

# Vectorized — returns an array of strings, one per geometry
reasons = shapely.is_valid_reason(gdf.geometry.values)

for idx, reason in zip(gdf.index[invalid_mask], reasons[invalid_mask.values]):
    print(f"  Feature {idx}: {reason}")

Step 3 — Apply Deterministic Repair

Route all geometries through shapely.make_valid(). Pass the underlying NumPy array directly to avoid GeoSeries overhead and trigger the C-accelerated ufunc path.

gdf = gdf.copy()  # avoid modifying the original in-place

# Vectorized repair — operates on entire geometry array at once
repaired = shapely.make_valid(gdf.geometry.values)

# Write only repaired geometries back into the GeoDataFrame
gdf.loc[invalid_mask, original_geom_col] = repaired[invalid_mask.values]

When dealing with complex bowtie polygons or overlapping rings, make_valid() may split a single polygon into multiple valid components, converting a Polygon to a MultiPolygon. For a focused walkthrough of this specific case, see Fixing Self-Intersecting Polygons in GeoPandas.

Step 4 — Post-Repair Validation

Never assume make_valid() succeeds unconditionally. Pathological inputs — zero-area polygons, extreme coordinate values near floating-point limits, geometries with thousands of self-intersections — occasionally produce empty geometries or remain technically invalid under stricter GEOS versions.

post_invalid_mask = ~gdf.geometry.is_valid
remaining_invalid = int(post_invalid_mask.sum())

empty_mask = gdf.geometry.is_empty
remaining_empty = int(empty_mask.sum())

if remaining_invalid > 0:
    print(f"WARNING: {remaining_invalid} geometries remain invalid after repair.")
    print("Indices:", gdf.index[post_invalid_mask].tolist())

if remaining_empty > 0:
    print(f"WARNING: {remaining_empty} empty geometries produced by repair.")
    print("Indices:", gdf.index[empty_mask].tolist())

Step 5 — Log and Route Fallbacks

Maintain a structured audit trail for every repair run. This supports compliance requirements, pipeline observability, and manual topology correction workflows.

import json
from datetime import datetime

audit = {
    "timestamp": datetime.utcnow().isoformat(),
    "source": "source_data.gpkg",
    "total_features": len(gdf),
    "invalid_pre_repair": invalid_count,
    "invalid_post_repair": remaining_invalid,
    "empty_post_repair": remaining_empty,
    "dead_letter_indices": gdf.index[post_invalid_mask | empty_mask].tolist(),
}

with open("geometry_repair_audit.json", "w") as f:
    json.dump(audit, f, indent=2)

Route dead-letter features — those remaining invalid or empty after repair — to a separate file for manual review rather than silently dropping them:

dead_letter = gdf[post_invalid_mask | empty_mask].copy()
if len(dead_letter) > 0:
    dead_letter.to_file("dead_letter_geometries.gpkg", engine="pyogrio")
    gdf = gdf[~(post_invalid_mask | empty_mask)].copy()

Advanced Patterns & Edge Cases

Precision Snapping to Eliminate False Invalidity

Floating-point precision errors are the most common cause of invalid geometry flags when merging datasets from different coordinate systems, digitising software, or jurisdictions. Two vertices that should coincide differ by 1e-10 degrees, creating a micro-gap that GEOS flags as an unclosed ring.

Apply shapely.set_precision() before the main repair pass to snap vertices to a uniform grid:

# Grid size in CRS units — for geographic CRS (degrees), 0.0001° ≈ 11 m at equator
# For projected CRS (metres), use 0.001 (1 mm) or smaller
gdf = gdf.copy()
gdf.geometry = shapely.set_precision(gdf.geometry.values, grid_size=0.0001)

Choose grid_size relative to the legal precision of your source data. Cadastral parcels warrant finer grids (0.00001°) than environmental monitoring polygons. For guidance on matching precision to CRS, see Handling Precision & Coordinate Rounding.

Collapsed Geometries and Empty Geometry Handling

Repair routines occasionally produce GEOMETRYCOLLECTION EMPTY or POINT EMPTY when topology violations are irreconcilable — for example, a polygon whose only non-self-intersecting component degenerates to a line at the chosen precision. Filter these explicitly before any spatial index construction:

empty_mask = gdf.geometry.is_empty
if empty_mask.any():
    empty_count = int(empty_mask.sum())
    dead_letter = gdf[empty_mask].copy()
    gdf = gdf[~empty_mask].copy()
    print(f"Removed {empty_count} empty geometries; routed to dead-letter store.")

Never pass empty geometries into spatial joins, rasterisation calls, or to_file(). GeoPandas and GDAL handle empty geometries inconsistently across backends, producing silent drops or file-write errors.

Multipart Geometry Splits and Schema Alignment

When make_valid() converts a Polygon to a MultiPolygon, row-level aggregations — area calculations, attribute joins keyed on a unique parcel ID — break silently because the count of rows no longer matches the count of source features.

Use explode() to normalise geometry types, but preserve original feature identifiers for traceability:

# Preserve the original row index as a column before exploding
gdf["original_fid"] = gdf.index

gdf = gdf.explode(index_parts=True).reset_index(drop=True)

# Count how many parts each original feature was split into
split_counts = gdf["original_fid"].value_counts()
multi_splits = split_counts[split_counts > 1]
if not multi_splits.empty:
    print(f"{len(multi_splits)} features were split into multiple parts:")
    print(multi_splits.head(10))

If your downstream schema requires single-part geometries, select the largest part per original FID:

gdf["area"] = gdf.geometry.area
gdf = (
    gdf.sort_values("area", ascending=False)
       .drop_duplicates(subset="original_fid", keep="first")
       .drop(columns=["area"])
       .reset_index(drop=True)
)

Performance Optimization for Large Datasets

Vectorized operations are mandatory for datasets exceeding 50k features. Row-wise .apply() calls introduce Python interpreter overhead and bypass Shapely 2.0’s C-accelerated ufunc paths.

Key rules:

1. Always pass gdf.geometry.values (a NumPy array of geometry objects) to Shapely 2.0 functions — not the GeoSeries itself.
2. Use boolean array indexing rather than iterating over rows.
3. For files exceeding available RAM, chunk reads via pyogrio’s skip_features / max_features parameters.

import geopandas as gpd
import shapely
import pyogrio

CHUNK_SIZE = 50_000
OUTPUT_PATH = "repaired_output.gpkg"

info = pyogrio.read_info("large_dataset.gpkg")
total_features = info["features"]
first_chunk = True

for offset in range(0, total_features, CHUNK_SIZE):
    chunk = gpd.read_file(
        "large_dataset.gpkg",
        engine="pyogrio",
        skip_features=offset,
        max_features=CHUNK_SIZE,
    )

    # Vectorized validity check on the underlying geometry array
    valid_flags = shapely.is_valid(chunk.geometry.values)
    invalid_idx = ~valid_flags

    if invalid_idx.any():
        chunk = chunk.copy()
        chunk.geometry.values[invalid_idx] = shapely.make_valid(
            chunk.geometry.values[invalid_idx]
        )

    # Remove empties
    empty_flags = shapely.is_empty(chunk.geometry.values)
    chunk = chunk[~empty_flags].copy()

    write_mode = "w" if first_chunk else "a"
    chunk.to_file(OUTPUT_PATH, engine="pyogrio", mode=write_mode)
    first_chunk = False
    print(f"Processed offset {offset}–{offset + CHUNK_SIZE}")

On a 2M-feature municipal parcel dataset, this chunked pattern sustains constant ~400 MB RAM usage regardless of total file size, versus loading the entire GeoDataFrame into memory which would require 8+ GB.

Integration into ETL Pipelines

Geometry repair must be idempotent, observable, and tightly coupled to schema validation gates. The following patterns make it safe to embed in automated ingestion triggers and CI/CD pipelines.

Schema enforcement hooks: Use pydantic or great_expectations to assert geometry validity before loading into analytical databases. Reject or quarantine batches that exceed a configurable invalidity threshold:

INVALIDITY_THRESHOLD = 0.05  # 5% of features

pre_repair_invalid_rate = invalid_count / len(gdf)
if pre_repair_invalid_rate > INVALIDITY_THRESHOLD:
    raise ValueError(
        f"Source data quality check failed: {pre_repair_invalid_rate:.1%} invalid "
        f"geometries exceeds threshold of {INVALIDITY_THRESHOLD:.1%}. "
        "Quarantine batch for manual review."
    )

Dead-letter queue pattern: Rather than dropping problematic features silently, route them to a quarantine store with attached GEOS error reasons. This prevents pipeline halts while preserving data for human review:

reasons = shapely.is_valid_reason(gdf.geometry.values)
gdf["_geos_error"] = reasons
dead_letter = gdf[post_invalid_mask].copy()
# Write to a quarantine layer for manual correction
if len(dead_letter) > 0:
    dead_letter.to_file("quarantine.gpkg", layer="invalid_geometries", engine="pyogrio")
gdf = gdf[~post_invalid_mask].drop(columns=["_geos_error"]).copy()

Cross-module consistency: Geometry repair is rarely isolated. It typically precedes Attribute Mapping & Schema Harmonization and raster alignment operations. Coordinate your repair logic with broader data quality gates documented in Automated Vector & Raster Cleaning Workflows to ensure consistent invalidity thresholds and audit formats across all ingestion stages.

CI/CD embedding: Add a geometry validation step to your data pipeline tests. A pytest fixture that loads a known-bad fixture and asserts zero post-repair invalids catches regression when GEOS or Shapely versions change:

import pytest
import geopandas as gpd
import shapely

def test_repair_leaves_no_invalid_geometries(bad_geom_fixture_path):
    gdf = gpd.read_file(bad_geom_fixture_path, engine="pyogrio")
    gdf = gdf.copy()
    gdf.geometry = shapely.make_valid(gdf.geometry.values)
    gdf = gdf[~gdf.geometry.is_empty].copy()
    assert gdf.geometry.is_valid.all(), "Post-repair validation failed"

Related

• Fixing Self-Intersecting Polygons in GeoPandas — targeted walkthrough for bowtie and figure-eight polygon repair with make_valid() and MultiPolygon splitting
• CRS Normalization Across Mixed Datasets — deterministic projection alignment to run before topology operations
• Spatial Deduplication & Topology Simplification — remove duplicate features and simplify over-detailed rings after geometry repair
• Handling Precision & Coordinate Rounding — choose grid sizes and rounding strategies that avoid introducing new invalidity
• Attribute Mapping & Schema Harmonization — align column names and types after geometry-type changes caused by explode()