Topology Validation & Repair in Python: A Production Pipeline Guide
Invalid geometry is the silent killer of spatial pipelines: self-intersections and bad ring orientation cause overlays and joins to fail or return wrong answers. This stage of Geospatial Data Ingestion & Processing Workflows detects and repairs those defects before analysis, sitting immediately after Shapefile & GeoJSON Parsing and before Spatial Joins & Merging. It is where every incoming layer is proven OGC-valid, snapped to a consistent precision, and gated at 100% validity — and it pairs with the deep dive on Fixing Self-Intersecting Polygons Programmatically.
Architecture & Data Structures
Topology validation ensures that vector geometries satisfy the strict constraints defined by the OGC Simple Features specification: rings must be closed and simple, exterior rings wound counter-clockwise and holes clockwise, and no two rings of the same polygon may cross. In Python these rules are enforced by the GEOS C++ engine, exposed through Shapely Geometry Operations and consumed in bulk through GeoPandas DataFrames. Every geometry carries two cheap predicates — is_valid and is_simple — and a diagnostic function, explain_validity, that returns the human-readable reason and the coordinate of the first defect.
from shapely.geometry import Polygon
from shapely.validation import explain_validity
# A "bowtie": the ring crosses itself, so it is invalid
bowtie = Polygon([(0, 0), (2, 2), (0, 2), (2, 0), (0, 0)])
print(bowtie.is_valid) # False
print(explain_validity(bowtie)) # Self-intersection[1 1]
The explain_validity output — Self-intersection[1 1] — names both the defect class and the (x, y) where GEOS detected it, which is what you log when you quarantine a feature. Common violation classes you will see in real data are self-intersections, ring self-intersections, holes lying outside the shell, nested holes, duplicate consecutive vertices, and unclosed rings.
Three properties sit next to validity and are routinely confused with it. Simplicity (is_simple) is the one that matters for lines: a LineString that crosses itself is perfectly valid — the OGC rules constrain polygon rings, not line paths — but it is not simple, and a network build, a linear-referencing step, or a route solver will produce nonsense from it. Polygons get simplicity for free, because a valid ring is by definition simple, so is_simple is only worth checking on line and multi-line layers. Emptiness is the second: is_valid returns True for an empty geometry, so a POLYGON EMPTY that arrived from a failed clip sails through the gate and then silently drops out of every join. Nullness is the third: a GeoDataFrame may hold None in the geometry column, which is not the same object as an empty polygon, and is_valid returns False for it in recent GeoPandas while older versions returned True. A production gate tests all three.
import geopandas as gpd
from shapely.geometry import LineString, Polygon
checks = gpd.GeoSeries([
LineString([(0, 0), (2, 2), (2, 0), (0, 2)]), # crosses itself
Polygon(), # POLYGON EMPTY
None, # missing geometry
])
print(checks.is_valid.tolist()) # [True, True, False]
print(checks.is_simple.tolist()) # [False, True, False]
print(checks.is_empty.tolist()) # [False, True, False]
The self-crossing line is valid and useless. That is the whole argument for treating "topology" as a set of contracts rather than a single boolean.
The fourth property is ring orientation, and it is the one that will bite you at the far end of the pipeline rather than in Python. GEOS does not care which way a ring winds: a polygon whose exterior runs clockwise and one whose exterior runs counter-clockwise are both valid and both report the same positive area. Downstream consumers disagree sharply. The Shapefile specification requires clockwise exterior rings and counter-clockwise holes; GeoJSON as standardised in RFC 7946 requires the opposite, the right-hand rule, with counter-clockwise exteriors; and the spherical-geometry engines behind BigQuery GIS and MongoDB's 2dsphere index interpret a reversed ring as the complement of the polygon — hand them a wound-backwards country and they will happily index the rest of the planet. Normalise orientation explicitly at the point of export rather than hoping the driver does it.
import shapely
# Shapely 2.1: vectorized over the array, CCW exterior + CW holes (RFC 7946)
oriented = shapely.orient_polygons(parcels.geometry.values, exterior_cw=False)
# Shapely 2.0 equivalent, per geometry
from shapely.geometry.polygon import orient
oriented = [orient(g, sign=1.0) for g in parcels.geometry]
Underneath all four properties sits GEOS's precision model. By default GEOS computes in full IEEE-754 double precision with no fixed grid, which means an intersection between two segments can land on a coordinate that neither input contained and that cannot be represented exactly — the mechanism by which a repair sometimes produces a new, tinier defect. Attaching an explicit grid with set_precision swaps that floating model for a fixed one, and is the reason the workflow below snaps after it repairs rather than before.
Environment Configuration & Dependency Resolution
Repair behaviour is tied to the GEOS version bundled with your wheels, so pin the stack explicitly. Shapely 2.0+ vectorises every predicate over NumPy arrays and requires GEOS 3.8 or newer; the make_valid method="structure" argument used later needs GEOS 3.10+.
python -m pip install "shapely>=2.0.3" "geopandas>=0.14" "pyogrio>=0.7" "pyproj>=3.6"
import shapely
import geopandas as gpd
print(shapely.__version__) # 2.0.3
print(shapely.geos_version) # (3, 12, 1) -> supports method="structure"
print(gpd.options.io_engine) # pyogrio (fast, Arrow-backed reader)
If shapely.geos_version is below (3, 10, 0), fall back to the default make_valid (the linework algorithm) rather than passing method="structure", which will otherwise raise. Keep GEOS consistent between your development machine and CI — a repair that produces a MultiPolygon under GEOS 3.12 can produce a GeometryCollection under 3.8, which changes how many rows your explode step emits.
Four GEOS boundaries govern what this stage can do, and it is worth knowing which one you are on before designing around a function that may not exist. GEOS 3.8 is the floor for Shapely 2's vectorized predicates. GEOS 3.9 introduced set_precision's fixed-precision overlay, which is what makes the snapping step deterministic rather than best-effort. GEOS 3.10 added the structure repair algorithm. GEOS 3.12 added coverage validation and coverage-aware simplification, the tools that address defects between features rather than inside them, described in Snapping and Simplifying Polygons Without Creating Gaps. On the Python side, Shapely 2.1 moved method and keep_collapsed onto the top-level shapely.make_valid, so a call written for 2.1 raises TypeError on 2.0.
The failure that wastes the most time here is a mixed installation. pip install shapely pulls a wheel with GEOS statically linked inside it; conda install shapely links against the separate geos package in the environment; installing GeoPandas from conda and Shapely from pip gives you two GEOS builds in one process, and which one answers depends on load order. The symptom is usually not a crash but a disagreement — a fixture that is valid on one machine and invalid on another, or a repair whose output type differs. Install the whole compiled stack from one channel and assert the versions in CI.
# One channel for everything that links GEOS, GDAL or PROJ
conda install -c conda-forge "shapely>=2.0.3" "geopandas>=0.14" "pyogrio>=0.7" \
"pyproj>=3.6" "libgdal-core>=3.8"
import shapely
assert shapely.geos_version >= (3, 10, 0), (
f"GEOS {shapely.geos_version} cannot run make_valid(method='structure')"
)
print(shapely.geos_capi_version_string) # 3.12.1-CAPI-1.18.1
Pin the resolved versions in a lockfile rather than a floor in requirements.txt. A geometry repair is a data transformation, and an unpinned GEOS means the same input file can produce a different output file six months later with nothing in your own repository having changed.
Vectorized Operations & Core Workflow
is_valid: clean rows flow through untouched while defects are quarantined, repaired, and snapped — then both lanes must clear the 100% validity gate before export.The pipeline reads a layer, isolates invalid features, and quarantines them so the ETL never halts on a single bad record. Load with pyogrio for Arrow-backed speed, then split the frame on the is_valid mask — a vectorised call that runs entirely in GEOS.
import geopandas as gpd
# Ingest with the fast Arrow-backed engine
parcels = gpd.read_file("cadastre.gpkg", layer="parcels", engine="pyogrio")
# A CRS is mandatory before any metric topology test
if parcels.crs is None:
raise ValueError("Dataset lacks a CRS. Assign one before topology checks.")
# Vectorised validity check across the whole GeoSeries (runs in GEOS)
invalid_mask = ~parcels.geometry.is_valid
print(f"Invalid geometries: {int(invalid_mask.sum())} / {len(parcels)}")
# Quarantine the defects; keep the clean rows flowing
invalid_parcels = parcels[invalid_mask].copy()
valid_parcels = parcels[~invalid_mask].copy()
Flagging invalid geometries upfront lets you audit problematic features without blocking the clean majority. For files too large to hold in RAM, read in row-group chunks from GeoParquet or use the bounding-box filter on read_file, and apply the same mask per chunk — see Cloud-Native Geospatial Formats for windowed reads.
The mask alone tells you how many features are broken; the quarantine is only actionable if it also records how. shapely.is_valid_reason is the vectorized counterpart to explain_validity and returns the diagnostic string for a whole array in one compiled pass, which turns the quarantine frame into a report you can group, count, and route.
import shapely
import pandas as pd
reasons = shapely.is_valid_reason(invalid_parcels.geometry.values)
invalid_parcels["defect"] = pd.Series(reasons, index=invalid_parcels.index)
# The class, without the coordinate, is what you route on
invalid_parcels["defect_class"] = (
invalid_parcels["defect"].str.split("[", regex=False).str[0].str.strip()
)
print(invalid_parcels["defect_class"].value_counts())
# Self-intersection 184
# Ring Self-intersection 27
# Hole lies outside shell 6
# Too few points in geometry component 2
That distribution is the most useful thing this stage produces, and it should be a metric you watch across runs rather than a number you glance at once. A layer that is 99.4 % self-intersections is a digitising or coordinate-truncation problem and make_valid will handle all of it. A layer where Too few points in geometry component dominates has features with fewer than four coordinates in a ring — degenerate records that no repair can rescue, because there is no shape to recover. Those belong in a rejects file with their identifiers, not in the repair lane. Splitting the quarantine by defect class before repairing is what stops a pipeline from silently converting bad records into plausible-looking empty ones.
UNREPAIRABLE = {"Too few points in geometry component"}
rejects = invalid_parcels[invalid_parcels["defect_class"].isin(UNREPAIRABLE)]
repairable = invalid_parcels[~invalid_parcels["defect_class"].isin(UNREPAIRABLE)]
if not rejects.empty:
rejects.to_file("rejects.gpkg", layer="unrepairable", driver="GPKG")
print(f"{len(rejects)} features cannot be repaired; written to rejects.gpkg")
Throughput is dominated by the repair, not the check. is_valid costs roughly a microsecond per moderate polygon and scales linearly with vertex count, so validating ten million features is a matter of seconds; make_valid is one to two orders of magnitude dearer per geometry because it has to node every edge against every other edge in the ring. This is precisely why the mask comes first: on a typical municipal cadastre where a fraction of a percent of features are broken, repairing only the flagged subset turns an hour into a few seconds, and it has the second benefit that byte-identical clean geometry passes through untouched, so a diff of the output against the input shows only the rows you meant to change.
Geometry / Data Processing Details
Repair is deterministic once features are isolated. shapely.make_valid decomposes a self-intersecting polygon into a valid MultiPolygon or GeometryCollection without discarding coordinates. GEOS offers two algorithms: the default "linework" preserves all input vertices by rebuilding valid rings from the original edges, while "structure" (GEOS 3.10+) assumes rings describe area and returns a cleaner result that drops zero-area slivers — usually what you want for cadastral polygons.
import geopandas as gpd
import shapely
from shapely import make_valid
def repair_topology(gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
"""Repair invalid geometry, preserving attribute rows 1:1."""
repaired = gdf.copy()
# Vectorised repair across the array; "structure" drops slivers
geom = repaired.geometry.values
needs_fix = ~shapely.is_valid(geom)
geom[needs_fix] = make_valid(geom[needs_fix], method="structure")
repaired["geometry"] = geom
# A repair can turn one Polygon into a MultiPolygon or collection.
# Explode to restore a 1:1 row-to-geometry mapping, then drop any
# non-polygonal fragments (stray lines/points from a bowtie repair).
repaired = repaired.explode(index_parts=False).reset_index(drop=True)
repaired = repaired[repaired.geometry.geom_type == "Polygon"]
return repaired
repaired_parcels = repair_topology(invalid_parcels)
assert repaired_parcels.geometry.is_valid.all()
The explode step matters: because a repaired bowtie becomes two triangles, exploding keeps each polygon on its own row with the parent attributes copied, so your attribute table stays synchronized with the corrected geometry. Filtering on geom_type == "Polygon" discards the degenerate line and point fragments that a GeometryCollection repair can leave behind. For the edge cases specific to crossed rings, the full walkthrough is in Fixing Self-Intersecting Polygons Programmatically.
Duplicate vertices and near-coincident points are a second, subtler class of defect. Snap the whole layer to a precision grid with set_precision so that vertices that differ only in the tenth decimal collapse into one, which removes the micro self-intersections that break later overlays.
from shapely import set_precision
# Snap to a 1 mm grid (in the layer's metric CRS units)
repaired_parcels["geometry"] = set_precision(
repaired_parcels.geometry.values, grid_size=0.001
)
Duplicate vertices that are exactly coincident are a cheaper problem than near-coincident ones, and worth removing separately because they inflate every subsequent operation without changing the shape. shapely.remove_repeated_points deletes consecutive coordinates closer together than a tolerance, which on a layer that has been through a raster-to-vector conversion or a GPS trace can cut the vertex count by a third with no visible effect.
import shapely
before = int(shapely.get_num_coordinates(repaired_parcels.geometry.values).sum())
repaired_parcels["geometry"] = shapely.remove_repeated_points(
repaired_parcels.geometry.values, tolerance=0.0
)
after = int(shapely.get_num_coordinates(repaired_parcels.geometry.values).sum())
print(f"Vertices: {before:,} -> {after:,}")
# Vertices: 4,812,904 -> 3,190,551
Keep tolerance=0.0 for a pure duplicate removal. Anything larger starts deleting real vertices, which is a simplification decision with consequences for shared borders and belongs in a coverage-aware step rather than here.
Line layers need noding, not make_valid. A road, river, or utility network is a LineString layer, and its characteristic defects are not ring violations at all: they are undershoots (a segment that stops just short of the junction it should meet), overshoots (a dangle projecting past the junction), and crossings without a shared node (two lines that intersect geometrically but have no coordinate in common, so a routing graph sees no connection). None of these makes a geometry invalid, and make_valid will return every one of them unchanged. The repair is shapely.node, which inserts a vertex at every intersection, usually preceded by a snap to close the gaps that are smaller than the survey tolerance.
import shapely
import geopandas as gpd
streets = gpd.read_file("streets.gpkg", engine="pyogrio").to_crs("EPSG:32633")
# Close sub-metre gaps first, then node every crossing into a shared vertex
merged = shapely.union_all(streets.geometry.values) # also dissolves duplicates
noded = shapely.node(merged)
segments = gpd.GeoDataFrame(
geometry=list(shapely.get_parts(noded)), crs=streets.crs
)
print(len(streets), "->", len(segments), "segments after noding")
# 18422 -> 21107 segments after noding
The row count grows because every crossing splits the two lines that meet there, which is exactly what a graph builder needs — see Street Network Analysis with OSMNX for what consumes the result. Note that union_all discards attributes: if the segments carry names, speed limits, or ownership, join them back by spatial predicate afterwards rather than trying to preserve them through the union.
Multipart geometry is a separate decision from validity. A repair can turn one Polygon into a MultiPolygon, and explode restores one geometry per row — but exploding duplicates the attribute row, so a parcel with an assessed value of 400,000 becomes two rows each claiming 400,000. If the attributes are additive, split them proportionally by area before exploding; if they are identifiers, add a part index so the original feature can be reconstructed. Deciding this per layer, and recording the decision, is more important than the repair itself, because an unnoticed row duplication propagates into every downstream count.
CRS Alignment & Projection Pipeline
Geometric validity and precision snapping are only meaningful in a metric coordinate system. In geographic coordinates (WGS84, EPSG:4326) a grid_size is measured in degrees, and intersection tests accumulate floating-point error near the poles. Reproject to a local projected CRS before validation using the transformation patterns from Coordinate Reference System Transformations. Avoid Web Mercator (EPSG:3857) for this — its area distortion makes sliver detection unreliable; pick the UTM zone or national grid that covers your data.
import geopandas as gpd
# Choose a projected CRS whose units are metres (UTM 33N here).
# Do NOT use EPSG:3857 — its scale distortion corrupts area/sliver tests.
METRIC_CRS = "EPSG:32633"
parcels_metric = parcels.to_crs(METRIC_CRS)
assert parcels_metric.crs.is_projected, "Topology checks need a projected CRS"
# pyproj respects axis order; GeoPandas geometries are always (x, y) = (lon, lat)
# ordering internally, so no always_xy juggling is needed at this layer.
Run validation, repair, and precision snapping in this metric CRS, then reproject the cleaned layer back to your storage CRS only at export time. Reprojecting a geometry can itself reintroduce tiny invalidities, so if you must transform after repair, re-run the validity gate afterwards.
The reason reprojection can break what you just fixed is that it is a nonlinear warp applied to vertices only. PROJ transforms each coordinate independently and the straight segment between two transformed vertices is not the transform of the straight segment between the originals — the further apart the vertices, the larger the divergence. Two edges that merely touched can therefore end up crossing, and a shape that was convex can develop a fold. The mitigation is to add vertices before the transform so each segment is short enough that the error is below your precision grid. segmentize does this in one vectorized call and is cheap compared to the repair it saves you.
import shapely
# Densify long edges before a projection change: no segment longer than 100 m
dense = shapely.segmentize(parcels_metric.geometry.values, max_segment_length=100)
parcels_metric["geometry"] = dense
parcels_wgs = parcels_metric.to_crs("EPSG:4326")
assert parcels_wgs.geometry.is_valid.all(), "reprojection reintroduced invalidity"
The grid size has to be restated whenever the units change, and this catches people out on the way back to storage. A grid_size of 0.001 means a millimetre in UTM metres and roughly 111 metres in EPSG:4326 degrees — the same literal, four orders of magnitude apart in effect. If a layer must be snapped in geographic coordinates, 1e-9 degrees is about 0.1 mm at the equator and is the sensible analogue; anything coarser will weld together vertices you needed. Choosing the metric CRS from the data rather than pinning a constant avoids a second class of error entirely, since a hard-coded zone silently distorts any batch that drifts out of it — the estimator is covered in Choosing a UTM Zone Automatically in Python. Two extents defeat every projected CRS and need handling before this stage rather than inside it: data crossing the antimeridian, whose bounding box wraps and whose polygons acquire a 360-degree-wide sliver on export, and truly global layers, which should be validated in an equal-area projection chosen once for the whole dataset.
Production Export & Integration
Close the stage with a hard gate: no layer leaves topology validation unless every geometry is valid. Rebuild the spatial index once so downstream Spatial Joins & Merging and Geometric Intersections & Overlays start from a clean STRtree.
import geopandas as gpd
# Pipeline gate: refuse to export residual invalid geometry
assert repaired_parcels.geometry.is_valid.all(), (
"Pipeline halted: residual invalid geometries detected"
)
# Warm the spatial index for the join/overlay stage that follows
_ = repaired_parcels.sindex
# GeoPackage for interchange; GeoParquet for cloud-native analytics
repaired_parcels.to_file(
"validated_parcels.gpkg", layer="clean_parcels", driver="GPKG", engine="pyogrio"
)
repaired_parcels.to_parquet("validated_parcels.parquet", geometry_encoding="WKB")
For PostGIS-backed pipelines, enforce validity at the database boundary too — a CHECK (ST_IsValid(geom)) constraint, or an ST_MakeValid call in the load query, stops an unvalidated feature from ever landing in a shared table. See PostGIS Integration with Python for the round-trip pattern. On datasets of millions of rows, prefer the array-level shapely.make_valid over GeoSeries.apply; the vectorised call bypasses Python iteration and is typically 5–10× faster.
import shapely
# Array-level repair for very large layers — no per-row Python loop
geom_array = repaired_parcels.geometry.values
repaired_parcels["geometry"] = shapely.make_valid(geom_array, method="structure")
The output format enforces its own rules, and they are not GEOS's rules. A layer that passes is_valid can still be rejected, or silently altered, by the driver that writes it. GeoPackage is the most permissive: any geometry type, mixed types in one layer, and no opinion on ring orientation. Shapefile forces clockwise exterior rings, has no way to distinguish Polygon from MultiPolygon (everything becomes a multipart record), truncates field names to ten characters, and cannot store an empty geometry at all. GeoJSON written to RFC 7946 wants counter-clockwise exteriors and coordinates in EPSG:4326 only, and expects polygons crossing the antimeridian to be split. GeoParquet records the geometry types present in its metadata, so a GeometryCollection that survived your explode filter will show up there and can break a reader that trusted the declared schema. Assert the type set before writing, rather than discovering the mutation on read-back.
types = set(repaired_parcels.geom_type.unique())
assert types <= {"Polygon", "MultiPolygon"}, f"unexpected geometry types: {types}"
assert not repaired_parcels.geometry.is_empty.any(), "empty geometry would be dropped"
Make the gate a test, not a comment. The counts this stage produces — features in, features invalid, features repaired, features rejected, vertices before and after — are pipeline metrics with the same status as row counts, and a change in them between runs is the earliest warning that an upstream source has shifted. Emit them as a structured artefact and assert on thresholds rather than on zero, so that a single new bad record does not stop the nightly run while a source that suddenly delivers 40 % broken geometry does.
import json
metrics = {
"features_in": int(len(parcels)),
"invalid_in": int(invalid_mask.sum()),
"repaired": int(len(repaired_parcels)),
"rejected": int(len(rejects)),
"invalid_rate": round(float(invalid_mask.mean()), 5),
}
with open("topology_metrics.json", "w") as fh:
json.dump(metrics, fh, indent=2)
assert metrics["invalid_rate"] < 0.05, (
f"invalid rate {metrics['invalid_rate']:.1%} exceeds the 5% threshold — "
"investigate the source before repairing"
)
When the layer never needs to enter Python at all — a nightly refresh of a shapefile drop, say — GDAL can do the whole stage from the shell, and on a very large file it will do it with a fraction of the memory: ogr2ogr -makevalid -nlt PROMOTE_TO_MULTI -t_srs EPSG:32633 clean.gpkg raw.shp. The -makevalid flag calls the same GEOS entry point as shapely.make_valid, so the result is identical; what you give up is the per-feature reporting, which is usually the reason to stay in Python.
Windows / Platform Edge Cases & Debugging
- PROJ database not found on Windows. A
PROJ: proj_create_from_databaseerror atto_crstime means thePROJ_LIB/PROJ_DATApath is unset in a conda env — reinstallpyprojfrom conda-forge rather than mixing pip and conda GEOS/PROJ builds. - GEOS version mismatch between machines. If a fixture that is valid locally fails in CI, print
shapely.geos_versionon both — aMultiPolygonunder GEOS 3.12 can be aGeometryCollectionunder 3.8, changing your explode count. Pin GEOS in the lockfile. method="structure"raisesGEOSException. You are on GEOS < 3.10; drop the argument to use the default linework algorithm, or upgrade Shapely wheels.- Empty geometries survive the gate.
is_validreturnsTruefor an empty geometry, so also filter~repaired_parcels.geometry.is_emptybefore export or joins will silently drop rows. make_validreturns aGeometryCollectionyou did not expect. A repair that yields stray lines or points means the input ring had a spur; explode and filter ongeom_typeas shown above rather than casting blindly toMultiPolygon.ImportErrorabout NumPy ABI after an upgrade. Shapely 2.0.x wheels are built against NumPy 1.x; installing NumPy 2 alongside them raises at import with a message about a module compiled for a different version. Upgrade to Shapely 2.0.4 or newer, which ships NumPy 2-compatible wheels, rather than pinning NumPy back.- Reading a
.shpfails on Windows with a path error. Shapefile paths longer than 260 characters hit the legacyMAX_PATHlimit that GDAL cannot work around; enable long paths in the registry, or stage the file to a short directory before reading. This surfaces disproportionately in topology work because quarantine and reject files tend to be written into deeply nested run directories. - The same file gives different validity counts on two machines. Print
shapely.geos_versionfirst; if they match, check the reader.fionaandpyogriocan differ in how they surface a geometry with aNaNcoordinate — one raises, the other returns a geometry that fails validation — so name the engine explicitly onread_filerather than depending on which is installed. - A repaired layer looks correct but every area is zero. The geometry is 2.5D and something in the chain dropped to a degenerate plane, or the layer is in a geographic CRS where areas are square degrees rounded to zero at display precision. Check
parcels.crs.is_projectedandshapely.has_z(...).any()before suspecting the repair. - Repair succeeds locally and the CI container runs out of memory.
make_validallocates a noded edge graph proportional to the square of the vertex count of a single geometry, so one 400,000-vertex coastline polygon can dominate peak memory regardless of how many rows the file has. Process the worst offenders separately, or simplify them before repair if the detail is beyond what the analysis needs.
Frequently Asked Questions
Where does this stage belong in the pipeline? Immediately after parsing and before anything that consumes geometry — joins, overlays, dissolves, tiling. The reason is economic rather than stylistic: a defect that survives into an overlay does not stay contained, because GEOS propagates it into every output feature the bad input touched, so one broken parcel can corrupt a whole neighbourhood's worth of results and the failure surfaces far from its cause. Validating at ingest costs a single vectorized pass over the layer and pins the blame to a named feature in a named file. If the layer is later reprojected, run the gate again after the reprojection, since the transform can reintroduce defects.
How can a layer where every geometry is valid still be wrong?
Because is_valid is a per-geometry predicate and most real requirements are relationships between geometries. A set of administrative boundaries that should tile a region with no gaps and no overlaps can consist entirely of individually perfect polygons and still have hairline cracks along every shared border and strips claimed twice. Nothing in this stage will detect that; it needs a coverage check, which GEOS 3.12 exposes and which is the subject of Snapping and Simplifying Polygons Without Creating Gaps. Treat per-feature validity as the floor, not the specification.
Should I repair automatically or quarantine and review?
Both, split by defect class. Self-intersections, ring self-intersections and repeated points are mechanical defects with a deterministic fix and no judgement involved — repair them in the pipeline and log the counts. Holes outside their shell, rings with too few points, and anything whose repaired area differs materially from its input deserve a human, because they usually indicate that the source data means something other than what the geometry says. The routing shown in the workflow section, which groups the quarantine by is_valid_reason before deciding, is what makes this practical at scale rather than a per-feature judgement call.
Is buffer(0) acceptable if it is faster?
It is not consistently faster, and speed is the wrong axis. A zero-width buffer produces valid output as a side effect of the offsetting machinery, with no contract about what it discards — it can delete slivers, drop a lobe, and shift area without a warning. make_valid has a documented behaviour per algorithm and lets you diff input against output. The detailed comparison, including where buffer(0) still has a role as a logged fallback, is in Fixing Self-Intersecting Polygons Programmatically.
Should validation run in Python or in PostGIS?
Run it in Python at ingest and enforce it in PostGIS at the boundary. The Python pass is where you can report per-feature diagnostics, route by defect class, and write a rejects file — none of which a database constraint gives you. The database constraint is where you guarantee that no other process, script, or colleague can insert an unvalidated geometry into a shared table. A CHECK (ST_IsValid(geom)) on the column plus an ST_MakeValid in any bulk-load query is the belt-and-braces pattern; the connection mechanics are in Connecting GeoPandas to PostGIS with SQLAlchemy.
Does invalid geometry actually break a spatial join, or only an overlay?
Both, in different ways. An overlay computes new geometry and will raise a TopologyException outright, which at least fails loudly. A join evaluates a predicate, and predicates on invalid geometry return answers — they are simply not trustworthy, because GEOS is asking "is this point inside a ring" of a ring that crosses itself and has no consistent inside. The result is a join that completes, produces plausible counts, and is wrong for the features that were broken. That silence is the argument for gating before Spatial Joins & Merging rather than after the first crash.