Spatial Data Serialization Patterns for SpatiaLite & GeoPackage

Without a clear encoding strategy, geometry data silently degrades: SRIDs are dropped during insertion, coordinate precision collapses after a VACUUM,…

Without a clear encoding strategy, geometry data silently degrades: SRIDs are dropped during insertion, coordinate precision collapses after a VACUUM, R-tree indexes fall out of sync after bulk writes, and round-trips through WKB produce geometries that fail ST_IsValid. These failures are invisible until a spatial predicate returns the wrong result or a mobile sync job pushes corrupt geometry to the field. This page documents production-tested serialization patterns for Python Integration & Database Workflows that prevent each of those failure modes before they reach production.

Prerequisites

Before implementing any pattern on this page, confirm the following are in place:

Concept & Specification Reference

Geometry values never travel as plain text between Python and SQLite. Each format imposes its own binary envelope, SRID stamp, and byte-order convention. Choosing incorrectly means the database receives a byte string it cannot interpret as a geometry.

FormatStandardInternal representationPython entry pointTypical use case
WKB (ISO)OGC SFA 1.2.1 §8.2.3Big- or little-endian byte array, no SRID headershapely.wkb.dumps() / GeomFromWKB()Native SpatiaLite inserts, maximum throughput
Extended WKB (EWKB)PostGIS extensionWKB + 4-byte SRID prefixshapely.wkb.dumps(include_srid=True)Cross-database portability when SRID must travel with bytes
GeoJSON textRFC 7946UTF-8 JSON objectgeom.__geo_interface__ / GeomFromGeoJSON()REST API payloads, human-readable debugging
GeoPackage Binary (GPB)OGC GeoPackage §2.1.3Magic GP header + flags byte + WKB bodyWritten by GDAL/Fiona/GeoPandas internallyGeoPackage geometry columns; do not write manually
SpatiaLite BLOBSpatiaLite internalCustom header + WKB; not raw WKBWritten by GeomFromWKB() / GeomFromText()SpatiaLite geometry columns; do not write manually

The critical constraint from this table: the raw bytes stored in a SpatiaLite or GeoPackage geometry column are not parseable by shapely.wkb.loads() directly. Reading must always go through ST_AsBinary(geometry) to strip the internal envelope. The GeoPackage Specification Deep Dive describes the GPB header layout in full if you need to parse it at the byte level.

The diagram below shows how geometry travels through each hop of a write–read round-trip:

Geometry serialization round tripOn write: Shapely geometry → wkb.dumps() → WKB bytes → GeomFromWKB(?, SRID) → internal DB BLOB. On read: internal DB BLOB → ST_AsBinary(geometry) → WKB bytes → wkb.loads() → Shapely geometry. Never pass the raw DB BLOB to wkb.loads().WRITEREADShapely geometryWKB bytesDB geometry columninternal BLOB (not raw WKB)wkb.dumps(geom, hex=False)GeomFromWKB(?, 4326)DB geometry columninternal BLOB (not raw WKB)WKB bytesST_AsBinary(geometry)
The stored column holds an internal SpatiaLite or GeoPackage BLOB, not raw WKB. All reads must pass through ST_AsBinary() before handing bytes to shapely.wkb.loads().

Step-by-step Implementation

Step 1: Validate the SpatiaLite extension before any geometry operation

Load and verify before touching geometry columns. Skipping this step means geometry function calls silently fail with no such function: GeomFromWKB.

python
# SpatiaLite context — requires mod_spatialite on LD_LIBRARY_PATH / PATH
import sqlite3

def verify_spatialite(conn: sqlite3.Connection) -> str:
    """Load mod_spatialite and return version string; raise on failure."""
    conn.enable_load_extension(True)
    conn.load_extension("mod_spatialite")
    conn.enable_load_extension(False)  # lock down after loading
    version: str = conn.execute("SELECT spatialite_version();").fetchone()[0]
    return version

Pair this with the Native sqlite3 Spatial Extensions guide for platform-specific path resolution on macOS, Windows, and Linux.

Step 2: Pattern A — WKB binary serialization (maximum throughput)

Well-Known Binary is the most storage-efficient path for SpatiaLite. GeomFromWKB parses the byte array, stamps the SRID, and writes the internal BLOB in a single SQL function call. This pattern is optimal for field-sensor ingestion and high-frequency mobile sync jobs.

python
# SpatiaLite context — WKB insert with explicit transaction and WAL tuning
import sqlite3
from shapely.geometry import Point
from shapely.wkb import dumps as wkb_dumps
from typing import Any

def insert_wkb(conn: sqlite3.Connection, table: str, geom: Point, attrs: dict[str, Any]) -> None:
    """Serialize Shapely geometry to WKB and insert into a SpatiaLite geometry column."""
    wkb_bytes: bytes = wkb_dumps(geom, hex=False, include_srid=False)

    conn.execute("PRAGMA journal_mode=WAL;")
    conn.execute("PRAGMA synchronous=NORMAL;")

    try:
        conn.execute(
            f"INSERT INTO {table} (geometry, id, name) VALUES (GeomFromWKB(?, 4326), ?, ?)",
            (wkb_bytes, attrs["id"], attrs["name"]),
        )
        conn.commit()
    except sqlite3.IntegrityError:
        conn.rollback()
        raise

GeomFromWKB does not validate topology — it accepts self-intersecting rings without complaint. Apply ST_IsValid checks in the verification step rather than trusting insertion to catch topology errors.

Step 3: Pattern B — GeoJSON string serialization (interoperability)

GeoJSON trades storage density for human readability and REST API compatibility. Use it when your pipeline exchanges data with web frontends or third-party tools that emit GeoJSON. GeomFromGeoJSON does not accept an SRID parameter, so wrap it in SetSRID.

python
# SpatiaLite context — GeoJSON insert via GeomFromGeoJSON + SetSRID
import sqlite3
import json
from shapely.geometry.base import BaseGeometry
from typing import Any

def insert_geojson(conn: sqlite3.Connection, table: str, geom: BaseGeometry, attrs: dict[str, Any]) -> None:
    """Serialize geometry as GeoJSON string and insert via SpatiaLite parser."""
    geojson_str: str = json.dumps(geom.__geo_interface__)

    try:
        conn.execute(
            f"INSERT INTO {table} (geometry, id, name)"
            " VALUES (SetSRID(GeomFromGeoJSON(?), 4326), ?, ?)",
            (geojson_str, attrs["id"], attrs["name"]),
        )
        conn.commit()
    except sqlite3.OperationalError:
        conn.rollback()
        raise

Validate incoming payloads against the RFC 7946 coordinate-member ordering rule (longitude first) before calling GeomFromGeoJSON. Lat/lon inversion is a common upstream error that the parser silently accepts, producing geometries in the wrong hemisphere.

Step 4: Pattern C — DataFrame-driven batch serialization (GeoPandas / Fiona)

For bulk ETL, row-by-row insertion becomes a bottleneck. GeoPandas delegates serialization to the underlying Fiona/OGR stack, which handles WKB encoding, GeoPackage binary envelope wrapping, spatial index creation, and gpkg_contents / gpkg_geometry_columns metadata table population automatically. Configure GeoPandas & GeoPackage Integration before calling to_file on large datasets.

python
# GeoPackage context — batch write via GeoPandas; Fiona/OGR handles GPB encoding
import os
import geopandas as gpd
import pandas as pd

def batch_write_gpkg(
    df: pd.DataFrame,
    output_path: str,
    layer_name: str = "features",
    crs: str = "EPSG:4326",
) -> int:
    """Write a DataFrame with a 'geometry' column to a GeoPackage layer.

    Returns the number of records written.
    """
    gdf = gpd.GeoDataFrame(df, geometry="geometry", crs=crs)

    if gdf.geometry.is_empty.any():
        raise ValueError("Empty geometries detected — clean before export.")
    if not gdf.geometry.is_valid.all():
        raise ValueError("Invalid geometries detected — run buffer(0) or make_valid().")

    mode = "a" if os.path.exists(output_path) else "w"
    gdf.to_file(output_path, layer=layer_name, driver="GPKG", mode=mode)
    return len(gdf)

For datasets exceeding available RAM, chunk the DataFrame into slices of 5,000–10,000 rows and call to_file with mode="a" on every slice after the first.

Step 5: Register and rebuild the spatial index after bulk inserts

Bulk inserts bypass R-tree population. After any batch write via executemany or direct SQL, rebuild the index explicitly:

sql
-- SpatiaLite context — register and rebuild spatial index after bulk insert
SELECT CreateSpatialIndex('features', 'geometry');
-- If the index already exists, delete and recreate:
-- SELECT DisableSpatialIndex('features', 'geometry');
-- SELECT CreateSpatialIndex('features', 'geometry');

Without this step, ST_Intersects, ST_Contains, and bounding-box queries fall back to full table scans. Check index state with SELECT * FROM geometry_columns WHERE f_table_name = 'features'; and confirm an rtree_features_geometry virtual table exists in sqlite_master.

Step 6: Read geometry back via ST_AsBinary

Never read the geometry column directly and pass the bytes to shapely.wkb.loads(). The column stores an internal envelope (SpatiaLite BLOB or GeoPackage Binary) that Shapely cannot parse. Always extract via ST_AsBinary:

python
# SpatiaLite context — safe WKB retrieval and Shapely reconstruction
import sqlite3
from shapely.wkb import loads as wkb_loads
from shapely.geometry.base import BaseGeometry

def read_geometry(conn: sqlite3.Connection, table: str, row_id: int) -> BaseGeometry:
    """Retrieve a stored geometry as a Shapely object via ST_AsBinary."""
    # Table names cannot be parameterised; validate the name before interpolation.
    row = conn.execute(
        f"SELECT ST_AsBinary(geometry) FROM {table} WHERE id = ?", (row_id,)
    ).fetchone()
    if row is None or row[0] is None:
        raise LookupError(f"No geometry found for id={row_id!r} in {table!r}")
    return wkb_loads(row[0])

Validation & Verification

After writing, run the following checks before promoting data to production:

python
# SpatiaLite context — topology and SRID audit
import sqlite3

def audit_geometry_column(conn: sqlite3.Connection, table: str, geom_col: str = "geometry") -> dict:
    """Return counts of invalid, null, and mixed-SRID geometries."""
    invalid = conn.execute(
        f"SELECT COUNT(*) FROM {table} WHERE ST_IsValid({geom_col}) = 0"
    ).fetchone()[0]
    null_geom = conn.execute(
        f"SELECT COUNT(*) FROM {table} WHERE {geom_col} IS NULL"
    ).fetchone()[0]
    srid_count = conn.execute(
        f"SELECT COUNT(DISTINCT ST_SRID({geom_col})) FROM {table}"
    ).fetchone()[0]
    return {"invalid": invalid, "null": null_geom, "distinct_srids": srid_count}

For GeoPackage output, also run ogrinfo to confirm the geometry column is registered and the envelope is correct:

bash
ogrinfo -al -so path/to/output.gpkg features

Look for Geometry Column = geometry, a reported SRID, and a non-zero Feature Count. A missing geometry column registration means gpkg_geometry_columns was not populated — this happens when writing via raw SQLite without going through OGR or SpatiaLite’s metadata functions. See the SpatiaLite Metadata Tables Explained reference for the table schema and required row structure.

Common Failure Modes & Fixes

Failure 1: no such function: GeomFromWKB

Diagnosis: The SpatiaLite extension was not loaded before the INSERT.

sql
-- Check which functions are available
SELECT name FROM pragma_function_list WHERE name LIKE '%geom%';

Fix: Call conn.enable_load_extension(True) and conn.load_extension("mod_spatialite") on the connection before executing any spatial SQL. If the extension path is wrong, set LD_LIBRARY_PATH or use an absolute path to the .so/.dylib/.dll. Consult Native sqlite3 Spatial Extensions for per-platform path resolution.

Failure 2: ST_IsValid returns 0 after round-trip

Diagnosis: Coordinate precision collapsed during serialization. WKB uses IEEE 754 double-precision floats; excessive decimal places are not the cause — loss of precision during CRS transformation before serialization usually is.

python
# Detect precision loss: compare pre- and post-serialization coordinates
from shapely.wkb import dumps, loads
original_coords = list(geom.coords)
rt_coords = list(loads(dumps(geom, hex=False)).coords)
assert original_coords == rt_coords, f"Precision drift: {original_coords[0]} vs {rt_coords[0]}"

Fix: Apply shapely.set_precision(geom, grid_size=1e-9) before serialization to snap coordinates to a stable grid. For topology repair, use geom.buffer(0) (effective for most self-intersections) or Shapely 2’s make_valid().

Failure 3: R-tree out of sync after bulk insert

Diagnosis: ST_Intersects returns zero rows on geometries you know should match. Check the R-tree state:

sql
-- SpatiaLite context — inspect R-tree for a table
SELECT COUNT(*) FROM rtree_features_geometry;
-- Compare with:
SELECT COUNT(*) FROM features;

Fix: If counts differ, rebuild:

sql
-- SpatiaLite context — full R-tree rebuild
SELECT DisableSpatialIndex('features', 'geometry');
SELECT CreateSpatialIndex('features', 'geometry');

Failure 4: Mixed SRIDs in geometry column

Diagnosis: Upstream data sources use different CRS (e.g., EPSG:4326 and EPSG:32632 mixed in one table). ST_Transform calls produce wrong results.

sql
-- SpatiaLite context — detect SRID drift
SELECT DISTINCT ST_SRID(geometry), COUNT(*) FROM features GROUP BY 1;

Fix: Enforce a single SRID at the Python layer using pyproj.Transformer before calling any insert function. Never mix projection transforms with serialization in the same step.

Failure 5: database is locked during concurrent writes

Diagnosis: Multiple threads share a single sqlite3.Connection, or a long-running read transaction blocks a write. SQLite’s write lock is exclusive at the database level.

Fix: Enable WAL mode (PRAGMA journal_mode=WAL) to allow concurrent readers alongside a single writer. For multi-threaded applications, configure Connection Pooling & Lifecycle Management with thread-local connections — never share a Connection across threads.

Performance Notes

  • Transaction batching: wrap every bulk insert in an explicit BEGIN / COMMIT block. SQLite autocommit forces a full fsync after each row. Batching 1,000–5,000 rows per transaction reduces write latency by an order of magnitude on spinning disk.
  • WAL + synchronous=NORMAL: PRAGMA synchronous=NORMAL with WAL mode is safe for most field-data workflows (no power-loss risk under a controlled shutdown) and reduces write overhead significantly compared to the default FULL.
  • VACUUM timing: run VACUUM only when the database has been heavily modified (many deletes or overwrites). On a large geometry table, VACUUM rewrites the entire file; schedule it during maintenance windows, not inline with inserts.
  • WKB vs GeoJSON storage cost: WKB encoding for a 10-vertex polygon is roughly 200 bytes; equivalent GeoJSON is approximately 600–900 bytes. At millions of features, this difference directly affects read throughput and index cache pressure.
  • Page cache: for read-heavy analytics over a large GeoPackage, increase PRAGMA cache_size to hold frequently accessed index pages in memory. A value of -65536 (64 MB in kibibytes) is a reasonable starting point for datasets under 1 GB.

More Pages in This Section