SpatiaLite & GeoPackage Automation with Python

How to Find Invalid Geometries in a GeoPackage

Run SELECT fid, STIsValidReason(geom) FROM layer WHERE geom IS NOT NULL AND STIsValid(geom) = 0 with modspatialite loaded — or, on a machine without the…

Run SELECT fid, ST_IsValidReason(geom) FROM layer WHERE geom IS NOT NULL AND ST_IsValid(geom) = 0 with mod_spatialite loaded — or, on a machine without the extension, decode each BLOB with Shapely and test is_valid in Python. Group the results by reason before doing anything else: the list of distinct reasons is short, and it tells you whether this is a data problem or a pipeline problem.

This page belongs to the Geometry Validity & Topology Repair guide, which covers the repair side. Here the scope is narrower: finding the rows, classifying them, and producing something you can act on.

Why This Matters

Invalid geometry does not announce itself. A container holding a thousand self-intersecting polygons opens without complaint, renders on a map, and returns rows from a bounding-box query — because the bounding box of an invalid polygon is perfectly well defined. What breaks is everything past the bounding box: the exact predicate, the area, the overlay. And because those all fail quietly, returning a plausible wrong answer rather than raising, the defect can travel a long way before anyone notices.

The scan is therefore worth running on any container you did not build yourself, and worth running as a gate on the ones you did. It is also the cheapest possible way to find out whether a producer upstream has a problem: a hundred distinct reasons across a hundred rows is data entry, and one reason across four thousand rows is a bug in whatever wrote them.

Prerequisites

  • Python 3.9+ with shapely 2.0+ for the extension-free path
  • mod_spatialite 5.0+ loadable from sqlite3 for the SQL path
  • A GeoPackage with at least one registered feature layer
  • Read access is sufficient — nothing here writes

Primary Method

python
# Scan every feature layer in a GeoPackage for invalid geometry
import sqlite3
from collections import Counter


def scan_invalid(gpkg_path: str) -> dict[str, Counter]:
    """
    Return {layer_name: Counter(reason -> count)} for every feature layer.
    Uses mod_spatialite, so the whole test runs inside SQLite.
    """
    conn = sqlite3.connect(f"file:{gpkg_path}?mode=ro", uri=True)
    conn.enable_load_extension(True)
    conn.load_extension("mod_spatialite")

    layers = conn.execute("""
        SELECT c.table_name, g.column_name
        FROM gpkg_contents c
        JOIN gpkg_geometry_columns g ON g.table_name = c.table_name
        WHERE c.data_type = 'features'
    """).fetchall()

    results = {}
    for table, geom_col in layers:
        rows = conn.execute(f"""
            SELECT ST_IsValidReason("{geom_col}") AS reason, count(*) AS n
            FROM "{table}"
            WHERE "{geom_col}" IS NOT NULL
              AND ST_IsValid("{geom_col}") = 0
            GROUP BY reason
        """).fetchall()
        if rows:
            results[table] = Counter({reason: n for reason, n in rows})

    conn.close()
    return results

Discovering the layers from gpkg_contents rather than hard-coding them matters: a container that gained a layer since the last scan is exactly the container most likely to have a new problem.

Scan path with and without mod_spatialiteTwo routes to the same answer. With the extension loaded, the validity test runs inside SQLite and only the summary crosses into Python, which is fast because no geometry is materialised as an object. Without it, every geometry blob is read into Python, stripped of its GeoPackage header, decoded by Shapely and tested there — the same answer at several times the cost, but available on any machine.the containergeometry BLOBswith mod_spatialitetest runs inside SQLiteonly the summary crosses outfastest by a wide marginShapely onlyevery blob decoded in Pythonheader stripped per rowworks on any machine
Both give the same verdict, because both call GEOS underneath. Only the amount of data crossing into Python differs.

Step-by-Step Walkthrough

1. Enumerate the layers to scan

python
layers = conn.execute("""
    SELECT c.table_name, g.column_name, g.geometry_type_name
    FROM gpkg_contents c
    JOIN gpkg_geometry_columns g ON g.table_name = c.table_name
    WHERE c.data_type = 'features'
    ORDER BY c.table_name
""").fetchall()
for t, col, gtype in layers:
    print(f"{t}.{col} ({gtype})")

2. Count invalid rows per layer

Start with counts rather than rows. A layer with zero invalid geometries needs no further work, and finding that out costs one query.

sql
-- GeoPackage + mod_spatialite: how bad is it, per layer?
SELECT 'parcels' AS layer,
       count(*) AS total,
       sum(CASE WHEN ST_IsValid(geom) = 0 THEN 1 ELSE 0 END) AS invalid
FROM parcels WHERE geom IS NOT NULL;

3. Group by reason

The reason string is what turns a count into a diagnosis. GEOS reports a small vocabulary — self-intersection, hole outside shell, nested holes, ring not closed, too few points — and the distribution across those is the signal.

python
from collections import Counter

reasons = Counter()
for (reason,) in conn.execute("""
    SELECT ST_IsValidReason(geom) FROM parcels
    WHERE geom IS NOT NULL AND ST_IsValid(geom) = 0
"""):
    # strip the trailing coordinate so the reasons group
    reasons[reason.split("[")[0].strip()] += 1

for reason, n in reasons.most_common():
    print(f"{n:>6}  {reason}")

Splitting off the bracketed coordinate is what makes the grouping useful — GEOS appends the location of the fault, so without it every row is its own unique “reason”.

4. Produce a work list with identifiers

Once the shape of the problem is known, extract the rows. Include whatever business key the layer carries, not just the row identifier: fid is meaningless to whoever has to look at the feature on a map.

sql
-- The work list: business key, reason, and where GEOS found the fault
SELECT fid, parcel_ref, ST_IsValidReason(geom) AS reason
FROM parcels
WHERE geom IS NOT NULL AND ST_IsValid(geom) = 0
ORDER BY parcel_ref
LIMIT 50;

5. Scan without the extension

On a machine where mod_spatialite will not load — a restricted interpreter, a minimal container image, a locked-down device — the same scan runs in Python. The only extra step is stripping the GeoPackage Binary header before Shapely sees the bytes.

python
# Extension-free scan: strip the GPB header, then let Shapely judge
import sqlite3
from shapely import from_wkb

ENVELOPE_LEN = {0: 0, 1: 32, 2: 48, 3: 48, 4: 64}


def gpb_to_wkb(blob: bytes) -> bytes:
    if blob[:2] != b"GP":
        raise ValueError("not a GeoPackage geometry BLOB")
    envelope = ENVELOPE_LEN[(blob[3] >> 1) & 0x07]
    return blob[8 + envelope:]


conn = sqlite3.connect("file:field.gpkg?mode=ro", uri=True)
invalid = []
for fid, blob in conn.execute("SELECT fid, geom FROM parcels WHERE geom IS NOT NULL"):
    geom = from_wkb(gpb_to_wkb(blob))
    if not geom.is_valid:
        from shapely.validation import explain_validity
        invalid.append((fid, explain_validity(geom)))

print(f"{len(invalid)} invalid of the rows scanned")

explain_validity is Shapely’s equivalent of ST_IsValidReason and returns the same GEOS vocabulary, so a work list produced either way is directly comparable. The header-stripping logic is the same as in How to Read GeoPackage Geometry Blobs in Python.

Two distributions of failure reason, and what each one meansTwo scan results compared. In the first, one reason accounts for nearly all failures across thousands of rows, which points at a single upstream bug producing the same fault repeatedly. In the second, a handful of rows fail for several different reasons, which points at manual data entry. The two call for completely different responses: fix the producer, or fix the rows.one reason, thousands of rows"Self-intersection" × 4,182everything else: 3a producer is emitting bad geometryrepair the layer, then fix the producermany reasons, few rowsfive distinct reasons × 12 rowsno dominant patternmanual capture, individually wrongrepair the rows, review the workflowthe count alone does not distinguish these — the grouping does
Grouping costs one extra clause and changes what you do next, which is why it belongs in the scan rather than in the follow-up.

Verification

Prove the scan can actually detect something before trusting a clean result. Insert a known-bad geometry, confirm the scan finds it, and remove it.

python
# The scan must reject a geometry you know is invalid
from shapely.geometry import Polygon
from shapely import to_wkb

bowtie = Polygon([(0, 0), (2, 2), (2, 0), (0, 2), (0, 0)])
assert not bowtie.is_valid, "fixture is not actually invalid"

conn.execute("BEGIN IMMEDIATE")
conn.execute(
    "INSERT INTO parcels (parcel_ref, geom) VALUES (?, GeomFromWKB(?, 4326))",
    ("TEST-INVALID", to_wkb(bowtie)),
)
found = conn.execute("""
    SELECT count(*) FROM parcels
    WHERE parcel_ref = 'TEST-INVALID' AND ST_IsValid(geom) = 0
""").fetchone()[0]
conn.execute("ROLLBACK")

assert found == 1, "the scan did not detect a known-invalid geometry"
print("scan verified")

Rolling back rather than deleting means the container is untouched even if the assertion fails partway.

Verifying the scan against a known-bad fixtureA bow-tie polygon is constructed in Python and asserted invalid before it goes anywhere near the database. It is inserted inside a transaction, the scan is run, and the scan must report exactly one invalid row. The transaction is then rolled back, so the container is unchanged whether the assertion passed or failed.build a bow-tieassert it is invalidinsert in a txnnever committedrun the scanmust find exactly 1ROLLBACKcontainer intacta scan that has never rejected anything has not been shown to detect anything
Rolling back rather than deleting means the container is untouched even when the assertion fails partway.

Alternative Approaches or Edge Cases

Scanning incrementally. On a large layer the full scan is a topology test per row, which is expensive enough that it should not run every cycle. Where the layer carries a change-log version, scan only rows above the last checked watermark — the pattern from Tracking Row Changes for Incremental GeoPackage Sync applies unchanged.

Non-geometry values in a geometry column. SQLite’s typing is per value, so a text value can sit in a BLOB column and ST_IsValid returns NULL rather than 0 for it. Add typeof(geom) = 'blob' to the scan and report anything else separately; it is a different defect with a different fix.

Layers with mixed geometry types. A layer declared GEOMETRY can hold points and lines alongside polygons, and validity is nearly meaningless for the first two. Filter on ST_GeometryType where the layer is mixed, so the work list contains only the rows a repair could act on.

Troubleshooting

sqlite3.OperationalError: no such function: ST_IsValid

Cause: mod_spatialite is not loaded on this connection, or extension loading is disabled in this interpreter build. Fix: Call enable_load_extension(True) then load_extension("mod_spatialite") before the query — see Using sqlite3 with SpatiaLite Functions for platform names, or use the extension-free path above.

ST_IsValidReason returns a different string for every row

Cause: GEOS appends the coordinate of the fault, so the strings are unique by construction. Fix: Split the reason at the opening bracket before grouping, as in step 3.

The scan reports zero invalid rows on a layer that renders wrongly

Cause: The problem is not validity. Ring orientation, duplicate vertices and coordinate precision all produce visibly odd results while passing ST_IsValid. Fix: Check orientation and vertex duplication separately — validity is one of several properties, and the parent guide covers where the boundary lies.

Frequently Asked Questions

How long should a full scan take?

Roughly the cost of reading and decoding every geometry once, plus a topology test per row — so it scales with total vertex count rather than with row count. A layer of simple parcels runs at tens of thousands of rows per second; one of detailed coastline polygons with tens of thousands of vertices each is far slower. If the scan is unexpectedly slow, check vertex counts before suspecting the query.

Should the scan run on a copy?

It does not need to — the scan only reads, and opening with mode=ro guarantees that. Working on a copy is worthwhile for a different reason: it lets you experiment with repairs on the same file you scanned without holding a connection to the live container or risking a stray write.

Can I scan a SpatiaLite container the same way?

Yes, with one change: enumerate the layers from geometry_columns rather than gpkg_contents, since the registries differ. The validity functions themselves are identical, because both containers call the same GEOS. Detecting which registry to read is covered in Reading Spatial Metadata with Python.