Fixing Self-Intersecting Polygons Programmatically in Python

Self-intersecting "bowtie" polygons parse cleanly but crash overlays, corrupt area sums, and fail to render — this guide repairs them deterministically with make_valid while keeping the coordinate reference system honest end to end. It is for anyone loading digitized parcels, floodplain boundaries, or merged administrative layers into a Python pipeline and hitting TopologyException downstream. It sits under Topology Validation & Repair in Geospatial Data Ingestion & Processing Workflows, and pairs with Best Practices for GeoJSON Validation, which catches the same defect at the file boundary.

Why This Approach / What Goes Wrong

A self-intersection occurs when the edges of a single polygon ring cross one another — the classic figure-eight bowtie, but also unclosed rings, spikes, and zero-area slivers born of coordinate truncation. Such geometry violates the OGC Simple Features rule that a polygon boundary must be a simple, non-crossing ring, so GEOS-backed predicates like intersects, union, and area either raise or return silently wrong answers.

Two naive fixes are common and both are traps:

Three bowtie repair strategies scored against four guarantees A matrix comparing make_valid run in a UTM zone, a zero-width buffer, and make_valid run on raw degrees. Only make_valid in UTM satisfies all four guarantees: it splits at the exact crossing, its area change is zero by design, it returns a predictable Polygon or MultiPolygon, and it works on a true metric plane. buffer(0) fails the first three and depends on the CRS for the fourth. Running make_valid in degrees is unstable, shifts area silently, varies in output type, and treats degrees as a flat plane. Three ways to fix a bowtie — only one of them is auditable What you can guarantee make_valid in UTM recommended buffer(0) the folk remedy make_valid in degrees raw EPSG:4326 input Splits at the exact crossing Area change is observable Output type is predictable Stable near poles / antimeridian exact node zero by design Polygon or Multi true metric plane a side effect silent area shift may fracture depends on CRS unstable silent area shift varies not a flat plane buffer(0) can close a bowtie, but only make_valid in a metric CRS gives a repair you can defend later.
Both shortcuts do produce a valid shape sometimes — neither tells you what it changed, which is the property a pipeline audit actually needs.

The correct approach is to project to an appropriate metric CRS — a local UTM zone, never Web Mercator, whose area and distance distortion makes it unfit for any metric operation — run the deterministic make_valid in that planar space, verify, then reproject back to the caller's original CRS. make_valid splits each intersecting edge exactly at its crossing coordinate and returns a valid Polygon or MultiPolygon, which is auditable in a way that buffer(0) is not.

Repairing a self-intersecting bowtie polygon with make_valid in a metric CRS On the left, an invalid bowtie polygon whose two diagonal edges cross at a single self-intersection point. An arrow labelled make_valid, run in a metric UTM CRS, points to the right, where the shape has become a valid MultiPolygon of two separate triangles that meet at a gap. A band below shows the three-step pipeline: reproject to UTM, then make_valid plus set_precision, then reproject back to the source CRS. Deterministic bowtie repair: split the crossing, keep the CRS honest Invalid: self-intersecting edges cross at one point make_valid() run in UTM · metric CRS splits each edge at its crossing Valid: MultiPolygon two clean, disjoint triangles 1 · reproject to UTM 2 · make_valid + set_precision 3 · reproject to source CRS
make_valid splits the bowtie at its exact crossing coordinate, returning an auditable MultiPolygon — run in an estimated UTM zone, never Web Mercator, then reprojected back.

Prerequisites

python -m pip install "geopandas>=1.0.0" "shapely>=2.0.0" "pyproj>=3.4.0"

Shapely 2.0 is required: the top-level vectorized shapely.make_valid and the stable explain_validity signature used here do not exist in the 1.x API, a break detailed in Shapely 1.x vs Shapely 2.0 Vectorization.

Two later boundaries change what make_valid will accept and what it returns. Shapely 2.1 added the method and keep_collapsed arguments to the top-level function, so shapely.make_valid(geom, method="structure", keep_collapsed=False) is a syntax error on 2.0.x and raises TypeError rather than doing something subtly different — a clean failure, but one that only appears on the machine with the older wheel. Underneath, the "structure" algorithm itself requires GEOS 3.10 or newer; a Shapely 2.1 wheel linked against GEOS 3.8 will import fine and then raise GEOSException at call time. Print both numbers before you rely on either:

import shapely

print(shapely.__version__)      # 2.1.2
print(shapely.geos_version)     # (3, 13, 1)

Keeping these identical between your workstation and CI matters more here than in most stages, because the shape of the repair output is version-dependent: a bowtie whose two lobes touch only at the crossing point comes back as a MultiPolygon on GEOS 3.12 and can come back as a GeometryCollection on 3.8, which changes how many rows a subsequent explode emits and therefore the row count of the file you ship.

Step-by-Step Implementation

1. Pick a metric CRS from the data, not a constant. Hard-coding one projected CRS distorts geometry that lands outside its zone. Let GeoPandas estimate the UTM zone from the layer's own bounds so the repair runs in a true metric space — the same auto-selection covered in Choosing a UTM Zone Automatically in Python.

import geopandas as gpd

parcels = gpd.read_file("parcels.gpkg")
if parcels.crs is None:
    raise ValueError("Layer has no CRS; set one before repair (e.g. parcels.set_crs('EPSG:4326')).")

# estimate_utm_crs returns the metric UTM zone covering the layer's centroid
metric_crs = parcels.estimate_utm_crs()   # e.g. EPSG:32633 (UTM 33N)

2. Reproject into that metric CRS. GeoPandas transforms with pyproj under the hood and honours axis order internally, so to_crs needs no always_xy flag — but keep the original CRS to return to.

source_crs = parcels.crs
parcels_utm = parcels.to_crs(metric_crs)

3. Flag the invalid geometries with GEOS. is_valid is vectorized; repairing only the invalid rows leaves clean geometry byte-for-byte untouched and makes the before/after count auditable.

invalid_mask = ~parcels_utm.is_valid
print(f"{invalid_mask.sum()} of {len(parcels_utm)} geometries are invalid.")

4. Repair with vectorized make_valid. Call the top-level Shapely function on the geometry array of the invalid rows. It runs in compiled C, splitting each self-intersection at its exact crossing coordinate.

import shapely

parcels_utm.loc[invalid_mask, "geometry"] = shapely.make_valid(
    parcels_utm.loc[invalid_mask, "geometry"].values
)

The default algorithm is "linework", and for a bowtie it does exactly what the figure above shows: it noded the crossing, rebuilds rings from the original edges, and keeps every input coordinate. The alternative, "structure", treats each ring as a description of area rather than a list of edges, and resolves the overlap by taking the symmetric result — which for a true bowtie is the same two lobes, but for the other defects hiding in the same column is not. A polygon with a hole that pokes outside its shell repairs to shell-minus-hole under "structure" and to a two-part MultiPolygon under "linework". A ring with a zero-width spike keeps the spike as a degenerate line under "linework" (inside a GeometryCollection) and drops it entirely under "structure". Neither is correct in the abstract; the question is whether your downstream consumer is a geometry(Polygon) column that will reject a collection, or an audit that must account for every coordinate the surveyor recorded.

from shapely import make_valid
from shapely.geometry import Polygon

# A shell with a spike: the last two vertices retrace the same segment
spiked = Polygon([(0, 0), (10, 0), (10, 10), (5, 10), (5, 16), (5, 10), (0, 10)])

print(make_valid(spiked).geom_type)                        # GeometryCollection
print(make_valid(spiked, method="structure").geom_type)    # Polygon
print(round(make_valid(spiked, method="structure").area, 1))   # 100.0

For cadastral and administrative polygons — anything where the ring is meant to bound a region and slivers are digitising noise — method="structure" with the default keep_collapsed=True is the pragmatic choice, because it returns a single typed polygon that a database column, a tile builder, or a spatial join will accept without a filtering step. Set keep_collapsed=False when a feature that degenerates to a line or a point should disappear rather than survive as a zero-area fragment.

4b. Audit what the repair changed before you accept it. make_valid is deterministic, but "deterministic" is not "harmless". The area of an invalid ring is not a meaningful number to begin with: GEOS evaluates the shoelace formula over whatever vertex sequence it was handed, and in a bowtie the two lobes are traversed in opposite directions, so their contributions partly cancel. A perfectly symmetric bowtie reports an area of zero before repair and its true total afterwards. That is a repair moving the number toward the truth, but it is still a number that changed on a feature someone may have signed off on, so capture the delta per row while both versions are in memory and still in metres.

import pandas as pd
import shapely

before = parcels_utm.loc[invalid_mask, "geometry"].values   # captured pre-repair
after = shapely.make_valid(before, method="structure")

report = pd.DataFrame({
    "parcel_id": parcels_utm.loc[invalid_mask, "parcel_id"].values,
    "type_after": [g.geom_type for g in after],
    "area_before_m2": shapely.area(before).round(2),
    "area_after_m2": shapely.area(after).round(2),
    "area_delta_m2": (shapely.area(after) - shapely.area(before)).round(2),
    "vertices_before": shapely.get_num_coordinates(before),
    "vertices_after": shapely.get_num_coordinates(after),
})

print(report.sort_values("area_delta_m2", key=abs, ascending=False).head(3))
#   parcel_id    type_after  area_before_m2  area_after_m2  area_delta_m2  ...
# 3    A-4417  MultiPolygon            0.00        8421.55        8421.55
# 7    A-5102       Polygon         3310.08        3310.11           0.03
# 1    A-0088       Polygon         1902.44        1902.44           0.00

Read the tail of that table, not the head. A large positive delta on a MultiPolygon row is a genuine bowtie whose lobes were cancelling — expected, explainable, and worth a line in the changelog. A delta of a few square centimetres on a row that stayed a Polygon is a sliver being cleaned up. What should alarm you is a negative delta of any size on a Polygon row: that means method="structure" decided part of the ring described a hole rather than land, and you have lost real area. Write this frame out beside the cleaned layer; on a cadastral or flood-designation dataset it is the difference between a repair and an unexplained change to a legal boundary.

5. Reproject back and pin precision. Return to the caller's CRS, then snap vertices to a fixed grid so re-truncation cannot reintroduce micro-intersections downstream. set_precision with a metre-scale grid is applied before the final reprojection while units are still metres.

# 1 mm grid in the metric CRS removes floating-point slivers deterministically
parcels_utm["geometry"] = shapely.set_precision(
    parcels_utm.geometry.values, grid_size=0.001
)
parcels_fixed = parcels_utm.to_crs(source_crs)
How a one millimetre precision grid quantises vertex coordinates A magnified axis shows grid nodes one millimetre apart along the easting 691240 point something metres. Four raw vertices A, B, C and D sit at arbitrary float positions above the axis, and arrows drop each one onto its nearest node. B and C, which differ by only a quarter of a millimetre, land on the same node and become a single vertex. A table below lists each raw coordinate, its snapped value, and how far it moved. What set_precision(grid_size=0.001) does to a vertex A B C D B and C collapse onto one node .996 .997 .998 .999 1.000 1.001 1.002 1.003 grid nodes 1 mm apart · easting 691240.xxx m in the metric CRS vertex raw x, metres snapped x result A 691240.99722 691240.997 moved 0.22 mm B 691240.99904 691240.999 moved 0.04 mm C 691240.99928 691240.999 same node as B, duplicate gone D 691241.00131 691241.001 moved 0.31 mm Coordinates that differ below the grid size stop being separate vertices, so they can no longer form a micro-crossing.
Snapping runs while the units are still metres: vertices closer together than the grid size merge into one node, which is exactly how the sub-millimetre self-intersections disappear.

6. Wrap it as a reusable, logged function. A pipeline needs the counts and the failure reasons, not just the output layer.

import logging
import shapely
import geopandas as gpd

logging.basicConfig(level=logging.INFO, format="%(levelname)s: %(message)s")


def repair_bowties(gdf: gpd.GeoDataFrame, grid_size: float = 0.001) -> gpd.GeoDataFrame:
    """Repair self-intersecting polygons in a metric UTM CRS, then reproject back.

    Runs make_valid in an estimated UTM zone (not Web Mercator), snaps to a
    fixed grid, and returns geometry in the caller's original CRS.
    """
    if gdf.crs is None:
        raise ValueError("GeoDataFrame must have a defined CRS before repair.")

    source_crs = gdf.crs
    work = gdf.to_crs(gdf.estimate_utm_crs())

    invalid_mask = ~work.is_valid
    if not invalid_mask.any():
        logging.info("All %d geometries already valid; nothing to repair.", len(work))
        return gdf

    logging.info("Repairing %d invalid geometries.", int(invalid_mask.sum()))
    work.loc[invalid_mask, "geometry"] = shapely.make_valid(
        work.loc[invalid_mask, "geometry"].values
    )
    work["geometry"] = shapely.set_precision(work.geometry.values, grid_size=grid_size)

    still_invalid = ~work.is_valid
    for idx in work.index[still_invalid]:
        logging.warning("Row %s still invalid: %s",
                        idx, shapely.validation.explain_validity(work.geometry.loc[idx]))

    return work.to_crs(source_crs)

Verification

Build a deliberate bowtie, repair it, and assert the result is valid and area-consistent. The bowtie below is the canonical self-intersecting square; make_valid returns a two-triangle MultiPolygon.

import geopandas as gpd
from shapely.geometry import Polygon

# A bowtie: the ring crosses itself between the 2nd and 4th vertices
bowtie = Polygon([(0, 0), (2, 2), (2, 0), (0, 2), (0, 0)])
assert not bowtie.is_valid   # GEOS flags the self-intersection

parcels = gpd.GeoDataFrame(
    {"parcel_id": ["A-1"]}, geometry=[bowtie], crs="EPSG:32633"
)

fixed = repair_bowties(parcels)

print(fixed.geometry.iloc[0].geom_type)   # MultiPolygon
print(fixed.is_valid.all())               # True
print(round(fixed.geometry.iloc[0].area, 3))  # 2.0  (two 1x1 triangles)

assert fixed.is_valid.all()
assert fixed.geometry.iloc[0].geom_type == "MultiPolygon"

Expected console output:

INFO: Repairing 1 invalid geometries.
MultiPolygon
True
2.0

Edge Cases & Debugging

Frequently Asked Questions

Should I use make_valid or buffer(0)? Use make_valid. buffer(0) works by constructing a zero-width offset curve and returning whatever the offsetting machinery produces, which happens to be valid — the validity is a by-product, not the contract, and nothing in the API tells you what it discarded. It routinely deletes zero-area slivers and can drop an entire lobe of a bowtie without a warning. make_valid has a documented behaviour per algorithm, splits at the exact crossing coordinate, and lets you diff the input against the output as shown in the audit step. The only reason to keep buffer(0) in a modern codebase is as an explicitly logged last resort for a geometry that make_valid returns empty for.

Which algorithm should I pass — linework or structure? Pick by what consumes the output. If the next stop is a typed database column, a vector-tile build, or anything that expects one polygon per row, use method="structure": it returns polygonal output and drops the degenerate lines and points that would otherwise arrive inside a GeometryCollection. If the next stop is an audit, a conflation against the original survey, or any process that must account for every recorded coordinate, use the default "linework", which preserves all input vertices and merely re-organises them into valid rings. When in doubt, run both on the invalid subset and compare the area report — a large disagreement between the two is itself a signal that the feature needs a human.

Why repair in UTM rather than just leaving the data in EPSG:4326? Because every decision make_valid makes is a planar computation. It solves line-segment intersections, compares areas to decide which side of a ring is inside, and snaps coordinates that fall within its precision model — all in a flat Cartesian space that assumes the two axes carry the same unit. In degrees they do not: near 60° latitude a degree of longitude is half a degree of latitude on the ground, so a defect that is symmetric in reality is skewed in the coordinate space where the repair is computed, and the crossing point GEOS finds is not the one on the ground. It also makes every tolerance meaningless, which is the specific reason set_precision in step 5 runs before the reprojection back. The zone selection itself is automatic, so there is no cost to doing it correctly.

Can I skip set_precision if make_valid already returned valid geometry? You can, and many pipelines do, but you are leaving the door open for the defect to come back. A geometry that is valid at full float64 precision can become invalid the moment anything rounds it — a Shapefile write that truncates, a GeoJSON export at six decimal places, a tile builder quantising to integer tile coordinates, or a reprojection. Snapping to an explicit grid while you are still in metres means you have already chosen where the rounding happens and verified that the result survives it. On a layer that never leaves GeoPackage or GeoParquet at full precision, the step is optional; on anything destined for a web map, it is not.

Some geometries are still invalid after the repair. What now? That is rare, and it almost always means the input was not what you think. Run shapely.validation.explain_validity on the survivors and read the coordinate it reports, remembering that it is in the working CRS — a coordinate near 691240, 5334100 is UTM metres, one near 13.4, 52.5 means the reprojection never happened. The two common causes are a geometry that was already a GeometryCollection before repair, which make_valid will not flatten for you, and rings whose coordinates contain NaN, which no repair can resolve. Quarantine those rows to a separate file with their identifiers rather than dropping them, and gate the export on the count of quarantined rows as described on the parent topic.

Does repairing polygons fix gaps and overlaps between neighbouring features? No, and this is the most common misunderstanding of the validity gate. is_valid is evaluated on one geometry at a time and knows nothing about the other rows in the table, so a layer in which every polygon is individually perfect can still have hairline gaps along shared borders and strips claimed by two features at once. Those are coverage defects, they will pass this gate untouched, and they need a different tool — see Snapping and Simplifying Polygons Without Creating Gaps.

Related