Implementing Shapely Geometry Checks in Python

You have a GeoDataFrame loaded from a GeoPackage, Shapefile, or GeoJSON source, and you need to know — with certainty — which features are geometrically invalid before passing data downstream. This page covers the exact Shapely 2.0+ API calls, their execution order, how to interpret diagnostic output, and how to apply automated repair within the Building Rule Engines with GeoPandas workflow. All patterns here are vectorized: no row-by-row Python loops.

Prerequisites

• Python 3.10+ — needed for X | Y union types in function signatures and improved dataclass semantics.
• Shapely 2.0+ (pip install shapely>=2.0) — Shapely 2.x re-built the API around NumPy-backed vectorized functions. Confirm with import shapely; print(shapely.__version__).
• GeoPandas 0.14+ — required for stable Shapely 2.x interop and chunksize support in read_file.
• PyPROJ 3.4+ — needed if you apply coordinate reference system (CRS) reprojection before running area or distance rules.
• Before running geometry checks, confirm your CRS is correctly defined. Geometry validity is evaluated in the XY plane, so mixed or undefined projections produce meaningless topology results. See coordinate reference system precision standards for projection best practices.
• Data volume: this page targets single-node, in-memory workloads up to roughly 500,000 features. Above that threshold, see the escalation guidance at the bottom of this page.

Step-by-Step Procedure

Step 1 — Load and confirm geometry column integrity

import geopandas as gpd
import shapely
import pandas as pd

gdf = gpd.read_file("parcels.gpkg", layer="parcels")

# Confirm CRS is defined before any spatial operation
assert gdf.crs is not None, "Input GeoDataFrame has no CRS — set it before validation"

# Detect missing (None/NaN) geometries — shapely functions will raise on these
null_geom_mask = gdf.geometry.isna()
print(f"Null geometries: {null_geom_mask.sum()}")

Verification: null_geom_mask.sum() should print 0 if the source data is complete. Records where this is True must be handled before passing to Shapely — they are separate from empty geometries (GEOMETRYCOLLECTION EMPTY) which Shapely can process.

Step 2 — Run vectorized structural validity and diagnostic checks

# Vectorized — operates on all geometries in one C-level GEOS call
is_valid = shapely.is_valid(gdf.geometry)
validity_reasons = shapely.is_valid_reason(gdf.geometry)

print(f"Valid: {is_valid.sum()} / {len(gdf)}")
print("Sample failure reasons:")
print(validity_reasons[~is_valid].value_counts().head(10))

shapely.is_valid_reason() returns strings from the GEOS library directly — for example "Self-intersection[404732.5 6789234.1]" or "Ring Self-intersection". These strings are the primary diagnostic signal for routing failures to the correct repair strategy.

Verification: Any row where is_valid is False and validity_reason is "Valid Geometry" indicates a Shapely/GEOS version mismatch — upgrade both libraries.

Step 3 — Check topological simplicity and empty geometries

# is_simple catches self-intersecting lines and degenerate rings
is_simple = shapely.is_simple(gdf.geometry)

# is_empty catches GEOMETRYCOLLECTION EMPTY and point/line collections with no coordinates
is_empty = shapely.is_empty(gdf.geometry)

print(f"Non-simple: {(~is_simple).sum()}")
print(f"Empty:      {is_empty.sum()}")

Note that is_simple() is distinct from is_valid(). A LineString that crosses itself fails is_simple() but may pass is_valid() under the Open Geospatial Consortium (OGC) Simple Features Access standard, because that standard only mandates self-intersection rules for ring-bearing geometries (polygons). Run both checks in all pipelines that handle mixed geometry types.

Step 4 — Enforce coordinate bounds

from shapely.geometry import box

# Define allowed extent — e.g. continental US in EPSG:4326
ALLOWED_BOUNDS = (-125.0, 24.0, -66.0, 50.0)
bbox = box(*ALLOWED_BOUNDS)

within_bounds = gdf.geometry.intersects(bbox)
print(f"Outside expected bounds: {(~within_bounds).sum()}")

Verification: Features from a different coordinate system will often appear at coordinates like (0, 0) or extreme latitudes. Bounds enforcement catches CRS errors that is_valid() will miss entirely.

Step 5 — Compile the validation report

def build_validation_report(gdf: gpd.GeoDataFrame, bounds: tuple | None = None) -> pd.DataFrame:
    """
    Return a per-feature validation report as a DataFrame.
    Requires Shapely 2.0+ and GeoPandas 0.14+.
    """
    null_mask = gdf.geometry.isna()

    # Operate only on non-null geometries; fill results for nulls afterward
    geom_array = gdf.geometry.values

    is_valid = shapely.is_valid(geom_array)
    validity_reasons = shapely.is_valid_reason(geom_array)
    is_simple = shapely.is_simple(geom_array)
    is_empty = shapely.is_empty(geom_array)

    if bounds:
        bbox = box(*bounds)
        within_bounds = gdf.geometry.intersects(bbox)
    else:
        within_bounds = pd.Series(True, index=gdf.index)

    report = pd.DataFrame({
        "null_geometry":    null_mask.values,
        "geometry_valid":   is_valid,
        "validity_reason":  validity_reasons,
        "topology_simple":  is_simple,
        "is_empty":         is_empty,
        "within_bounds":    within_bounds.values,
        "geometry_type":    gdf.geometry.geom_type.values,
    }, index=gdf.index)

    report["validation_failed"] = (
        report["null_geometry"] |
        ~report["geometry_valid"] |
        ~report["topology_simple"] |
        report["is_empty"] |
        ~report["within_bounds"]
    )

    return report

Verification: Call report["validation_failed"].value_counts() — if every row is False, your dataset is clean. Any True rows are candidates for repair or quarantine.

Step 6 — Apply targeted repair with make_valid()

def repair_geometries(
    gdf: gpd.GeoDataFrame,
    report: pd.DataFrame,
    max_area_loss_fraction: float = 0.05,
) -> tuple[gpd.GeoDataFrame, pd.DataFrame]:
    """
    Apply make_valid() to invalid rows and re-validate.
    Returns (repaired_gdf, updated_report).
    """
    failed_mask = report["validation_failed"] & ~report["null_geometry"]
    if not failed_mask.any():
        return gdf.copy(), report.copy()

    repaired_gdf = gdf.copy()
    original_areas = gdf.loc[failed_mask].geometry.area

    repaired_geoms = shapely.make_valid(gdf.geometry[failed_mask].values)
    repaired_gdf.loc[failed_mask, "geometry"] = repaired_geoms

    # Area-loss guard: reject repairs that shrink features more than the threshold
    new_areas = repaired_gdf.loc[failed_mask].geometry.area
    area_loss = (original_areas - new_areas) / original_areas.clip(lower=1e-9)
    excessive_loss = area_loss > max_area_loss_fraction
    if excessive_loss.any():
        n = excessive_loss.sum()
        print(f"WARNING: {n} repaired feature(s) lost >{max_area_loss_fraction*100:.0f}% area — quarantine recommended")

    # Re-validate repaired rows
    post_valid = shapely.is_valid(repaired_gdf.loc[failed_mask].geometry.values)
    updated_report = report.copy()
    updated_report.loc[failed_mask, "geometry_valid"] = post_valid
    updated_report.loc[failed_mask, "validation_failed"] = ~post_valid

    fixed = post_valid.sum()
    still_invalid = (~post_valid).sum()
    print(f"Repair complete: {fixed} fixed, {still_invalid} still invalid (quarantine these)")

    return repaired_gdf, updated_report

The area-loss guard is critical in cadastral and infrastructure datasets where geometry is legally or financially significant. make_valid() can quietly drop sliver polygons or collapse near-degenerate rings — both of which may represent real features.

Interpreting Results

The validity_reason column maps directly to GEOS error codes. Use this table to select the appropriate repair strategy:

When make_valid() changes the geometry type (e.g. Polygon → MultiPolygon), downstream spatial joins and area calculations may produce unexpected results. For categorizing and prioritizing spatial errors, flag type-changed features as warnings even if the repair technically succeeded.

Gotchas & Edge Cases

Floating-point coordinate drift causes false-positive self-intersections. Coordinates exported from CAD systems or legacy shapefiles sometimes suffer from sub-millimetre floating-point noise that creates apparent self-intersections. Apply shapely.set_precision(geom, grid_size=0.001) before validation to snap coordinates to a stable grid. Confirm the grid size matches the tolerance documented in your data’s spatial data quality policies.

is_empty() and isna() catch different categories of missing geometry. shapely.is_empty() returns True for GEOMETRYCOLLECTION EMPTY and similar WKT forms — the geometry object exists but contains no coordinates. GeoDataFrame.geometry.isna() catches Python None and NumPy NaN. Shapely functions will raise a TypeError or produce NaN results if passed None values, so always run the isna() check first and exclude those rows.

3D coordinates are invisible to validity checks. is_valid() operates purely on XY projection. If your pipeline ingests LiDAR-derived features or survey data with Z coordinates, add a separate bounds rule using shapely.get_coordinates(geom, include_z=True)[:, 2] to verify elevation ranges.

Mixed geometry types require type-filtered rules. A GeoDataFrame containing both Polygon and LineString features will return mixed geometry_type values. Rules like minimum-area or ring-orientation checks are only meaningful for polygon types. Partition the DataFrame by geom_type and apply type-specific rules to each partition, then merge reports before routing to remediation.

make_valid() on large datasets can split features across the anti-meridian. For global datasets that cross the 180° meridian (e.g. Pacific Ocean boundaries), make_valid() may produce geometries that wrap unexpectedly. Clip to a projected CRS that avoids the anti-meridian before running validation, then reproject to the storage CRS afterward.

When to Escalate

This Shapely-based approach is appropriate for in-memory, single-node workloads. Escalate to a more powerful method in these situations:

• Feature count exceeds 500,000: Shapely 2.x vectorized functions remain single-threaded at the Python layer. Move to scaling GeoPandas validation with Dask or Apache Sedona for distributed execution.
• Cross-feature topology rules are required: Checks like “no two parcels overlap” or “road network must be fully connected” require spatial joins across features — not just per-feature validity. PostGIS ST_IsValid combined with ST_Intersects on a spatially indexed table outperforms GeoPandas for these predicates at scale.
• Continuous ingestion with sub-second latency requirements: Synchronous in-process validation cannot meet real-time SLAs. Move to an asynchronous validation workflow using Celery or a message queue to decouple ingestion throughput from validation latency.
• Audit trail and lineage are compliance requirements: Shapely checks produce in-memory results. For regulated datasets (INSPIRE, ISO 19157), integrate with a lineage-aware pipeline that records which rule version ran against which dataset version. See the compliance framework alignment guidance for audit trail requirements.

Related:

• Building Rule Engines with GeoPandas — the parent rule engine pattern this page extends
• Geometry Validity Checks for Vector Data — OGC validity concepts and GEOS constraint definitions
• Categorizing and Prioritizing Spatial Errors — severity tiers and routing logic for the failures this page detects
• Scaling GeoPandas Validation with Dask — what to do when datasets outgrow single-node Shapely checks

Back to Building Rule Engines with GeoPandas