GeoPandas & GeoPackage Integration

Without explicit control over drivers, spatial indexing, and connection lifecycles, GeoPandas pipelines silently corrupt GeoPackage files, leak file descriptors, and produce data that fails OGC validation checks. These failures surface only at deployment time — on a field device with no internet connection or in a nightly batch job that overwrites its only backup copy. This guide details a production-tested approach to reading, transforming, and writing GeoPackage data with GeoPandas, covering every decision point where a naive implementation breaks.

GeoPandas sits in the Python Integration & Database Workflows stack as the highest-level abstraction: it wraps GDAL/OGR geometry I/O, coordinate reference system transformations via pyproj, and Shapely geometry operations into a DataFrame interface. Below it, the GeoPackage format depends on a correctly registered GDAL driver, an rtree_<layer>_rowid spatial index, and mandatory gpkg_contents and gpkg_geometry_columns metadata rows. When any of those lower layers misbehave, GeoPandas gives you no warning — you just get empty results or silently mismatched geometry counts.

Prerequisites #

Before implementing spatial I/O, verify that every layer of the stack is correctly installed and wired together. GeoPandas 1.0 switched from fiona to pyogrio as its default I/O backend; the old fiona path still works but is slower and lacks vectorized geometry transfer.

• Python 3.9 or higher
• geopandas>=1.0
• pyogrio>=0.7.2 (installs GDAL wheels automatically on most platforms)
• shapely>=2.0 (required for vectorized geometry predicates)
• GDAL 3.4 or higher with GeoPackage driver enabled
• sqlite3 compiled with SpatiaLite extensions (bundled with most GDAL distributions; needed only if you call ST_* functions directly — see Native sqlite3 Spatial Extensions)

Validate driver registration before writing any pipeline code:

import pyogrio
import geopandas as gpd

# Verify GPKG driver is registered in GDAL/OGR.
drivers = pyogrio.list_drivers()
assert "GPKG" in drivers, "GeoPackage driver not registered — reinstall pyogrio"
print(f"pyogrio {pyogrio.__version__} | geopandas {gpd.__version__}")

If the assertion fails, your GDAL build lacks SQLite/GeoPackage support. Reinstall pyogrio from the conda-forge channel (conda install -c conda-forge pyogrio) or build GDAL with --with-sqlite3 --with-spatialite. The GDAL GeoPackage Driver Reference documents all compile-time options.

Concept & Specification Reference #

Understanding the GeoPackage internal schema prevents the most common data-integrity mistakes.

When GeoPandas calls to_file() with SPATIAL_INDEX="YES", GDAL populates all five of these table groups automatically. Writing via raw sqlite3 bypasses GDAL and leaves gpkg_contents and the R-tree unpopulated — making the file unreadable by QGIS, ArcGIS, or any standards-compliant consumer. If you need raw SQLite access alongside GeoPandas, route writes through the Transaction Scoping & Rollback Strategies layer rather than bypassing GDAL.

The mandatory gpkg_geometry_columns.geometry_type_name must be one of the OGC geometry type names: POINT, LINESTRING, POLYGON, MULTIPOINT, MULTILINESTRING, MULTIPOLYGON, GEOMETRYCOLLECTION, or GEOMETRY (for mixed-type layers). Storing a MultiPolygon in a layer declared as POLYGON is a spec violation that causes silent read failures in some clients.

Step-by-Step Implementation #

1. Metadata Inspection and Layer Discovery #

Loading an entire dataset into memory solely to check its schema wastes RAM and increases cold-start latency on edge devices. Query layers and geometry types directly from the file header instead:

# -- GeoPackage context: field_survey.gpkg
import pyogrio

# list_layers returns an ndarray of [layer_name, geometry_type] rows.
layers = pyogrio.list_layers("field_survey.gpkg")
print("Available layers:", layers[:, 0].tolist())

# read_info returns schema without loading any geometry data.
for layer_name in layers[:, 0]:
    info = pyogrio.read_info("field_survey.gpkg", layer=layer_name)
    print(
        f"  {layer_name}: {info['geometry_type']} | "
        f"CRS: {info['crs']} | fields: {info['fields']}"
    )

This pattern enables dynamic pipeline routing: downstream processors can adapt to incoming geometry types or field schemas without hardcoding assumptions about layer names.

2. Filtered Reads and Memory Optimization #

Mobile and edge deployments frequently operate under strict memory constraints. pyogrio exposes native bounding-box and column filtering that pushes predicates down to the SQLite engine, bypassing the Python heap entirely for rows outside the filter:

import geopandas as gpd

# Define spatial extent as (minx, miny, maxx, maxy).
bbox = (-122.5, 37.7, -122.3, 37.9)

gdf = gpd.read_file(
    "field_survey.gpkg",
    layer="survey_points",
    bbox=bbox,
    columns=["asset_id", "condition_score", "geometry"],
    engine="pyogrio",
)
print(f"Loaded {len(gdf)} features within bbox")

When you need custom SQL predicates that the OGR filter layer cannot express — for example, filtering on a computed spatial relationship or a high-cardinality attribute with a covering index — drop to Native sqlite3 Spatial Extensions for direct SQL execution, then reconstruct a GeoDataFrame from the result rows.

3. Spatial Transformations and Geometry Validation #

Raw field data frequently contains topological errors, self-intersecting rings, or mismatched coordinate reference systems. Shapely 2.0 provides vectorized operations that execute at C-speed, but validation must happen before any write to prevent invalid geometry rows from being committed:

# Ensure consistent CRS — always reproject before any spatial join or write.
if gdf.crs.to_epsg() != 4326:
    gdf = gdf.to_crs("EPSG:4326")

# Identify and repair invalid geometries before writing.
invalid_mask = ~gdf.geometry.is_valid
if invalid_mask.any():
    print(f"Repairing {invalid_mask.sum()} invalid geometries via buffer(0)...")
    gdf.loc[invalid_mask, "geometry"] = (
        gdf.loc[invalid_mask, "geometry"].buffer(0)
    )
    # Verify repair succeeded — buffer(0) can fail on extreme degeneracy.
    still_invalid = ~gdf.geometry.is_valid
    if still_invalid.any():
        raise ValueError(
            f"{still_invalid.sum()} geometries could not be repaired. "
            "Inspect with shapely.validation.explain_validity()."
        )

The buffer(0) technique repairs the most common topological errors (self-intersecting rings, duplicate vertices). For production pipelines feeding regulatory or compliance systems, log invalid geometries to an audit table before repair rather than silently mutating them. Reserve silent repair only for pipelines where the source is trusted field collection equipment with known GPS jitter patterns.

4. Indexed Writes and Transaction Control #

Writing to GeoPackage requires explicit spatial index creation and correct layer creation options. GeoPandas’ to_file() abstracts most of this when you use the pyogrio engine and pass GDAL creation options as keyword arguments:

gdf.to_file(
    "processed_survey.gpkg",
    layer="validated_assets",
    driver="GPKG",
    engine="pyogrio",
    # GDAL layer creation options as keyword arguments (pyogrio >=0.7).
    SPATIAL_INDEX="YES",
    GEOMETRY_NAME="geom",
)

SPATIAL_INDEX="YES" triggers GDAL to build the rtree_validated_assets_geom virtual table during the write phase. Omitting this option produces a valid file that reads correctly but degrades to full-table scans on any spatial query, which is catastrophic for large layers on low-power field devices.

For multi-layer transactions — where you need to update several layers atomically or roll back the entire write on failure — combine to_file() with the explicit BEGIN/COMMIT/ROLLBACK patterns described in Transaction Scoping & Rollback Strategies. To reuse a pre-warmed connection across multiple writes without re-paying the extension-loading cost, apply the patterns from Connection Pooling & Lifecycle Management.

5. Post-Write Validation and Integrity Checks #

Never assume a write succeeded without re-reading and asserting correctness. Re-open the layer and compare row counts, geometry validity, and CRS identity:

verify_gdf = gpd.read_file(
    "processed_survey.gpkg",
    layer="validated_assets",
    engine="pyogrio",
)

assert len(verify_gdf) == len(gdf), (
    f"Row count mismatch: wrote {len(gdf)}, read back {len(verify_gdf)}"
)
assert verify_gdf.geometry.is_valid.all(), (
    "Invalid geometries persisted after write"
)
assert verify_gdf.crs.to_epsg() == 4326, (
    "CRS drift detected: expected EPSG:4326"
)
print(f"Write validated: {len(verify_gdf)} features, all valid, EPSG:4326")

This final checkpoint catches silent truncation (SQLite BLOB limit exceeded), CRS drift (missing gpkg_spatial_ref_sys row), and driver-level serialization errors before data propagates to downstream consumers or synchronizes to field devices.

Validation & Verification #

After writing a GeoPackage with GeoPandas, verify the OGC-mandated metadata tables are populated correctly using ogrinfo and direct SQLite inspection:

# Confirm layer appears in gpkg_contents and reports a valid bounding box.
ogrinfo -al -so processed_survey.gpkg validated_assets

# Inspect mandatory metadata tables directly (no SpatiaLite needed).
sqlite3 processed_survey.gpkg \
  "SELECT table_name, data_type, srs_id, min_x, min_y, max_x, max_y
   FROM gpkg_contents;"

# Confirm spatial index virtual table exists.
sqlite3 processed_survey.gpkg \
  "SELECT name FROM sqlite_master WHERE type='table' AND name LIKE 'rtree_%';"

# Check geometry column registration.
sqlite3 processed_survey.gpkg \
  "SELECT table_name, column_name, geometry_type_name, srs_id
   FROM gpkg_geometry_columns;"

All four queries should return non-empty results with consistent srs_id values. If gpkg_contents is missing a row for validated_assets, GDAL silently dropped the write — re-check driver version compatibility.

Performance Notes #

WAL Mode for Offline-First Writes #

GeoPackage’s SQLite foundation makes it exceptionally portable but introduces a single-writer constraint. SQLite allows unlimited concurrent readers but serializes all writes to one connection at a time. For offline-first mobile platforms handling intermittent sync events, database is locked errors are the most frequent production failure mode.

Enable WAL (Write-Ahead Logging) mode before heavy write operations to allow concurrent reads during the write and to reduce lock-contention latency:

import sqlite3

# Enable WAL mode; survives file close/reopen.
conn = sqlite3.connect("processed_survey.gpkg")
result = conn.execute("PRAGMA journal_mode=WAL;").fetchone()
print(f"Journal mode: {result[0]}")  # Should print 'wal'
conn.close()

WAL mode enables readers to access the last committed snapshot while a writer is in progress. It does not eliminate the single-writer constraint, but it eliminates read-blocking by writers and substantially reduces SQLITE_BUSY frequency during mixed read/write field-sync workloads.

Page Cache Sizing and Batch Write Throughput #

SQLite’s default page cache is 2 MB, which is insufficient for batch writes of large geometry datasets. Increasing the cache size before a bulk import reduces the number of OS-level read-backs during index rebuild:

conn = sqlite3.connect("processed_survey.gpkg")
# Set 64 MB page cache (negative value = kibibytes).
conn.execute("PRAGMA cache_size = -65536;")
conn.execute("PRAGMA synchronous = NORMAL;")  # safe with WAL mode
conn.close()

R-tree Index Rebuild Cost #

Rebuilding an R-tree index after a large bulk insert is the most expensive single operation in a GeoPackage write pipeline. On a 500,000-feature polygon layer, index rebuild can take 30–90 seconds on an ARM SoC. To minimize the cost:

1. Write features in bounding-box order (sort by ST_MinX(geometry) before write).
2. Disable the spatial index during bulk insert, then rebuild once at the end — this is not directly exposed by to_file() but is achievable through GDAL layer creation with SPATIAL_INDEX="NO" followed by a manual CREATE VIRTUAL TABLE rebuild via sqlite3.
3. Run VACUUM after index rebuild if the file will be distributed to field devices — it compacts the R-tree pages and reduces read amplification.

For the patterns governing index rebuild after bulk inserts, consult the Native sqlite3 Spatial Extensions reference which covers DisableSpatialIndex / CheckSpatialIndex SpatiaLite functions and their GeoPackage equivalents.

Common Failure Modes & Fixes #

1. GPKG Driver Not Found #

Symptom: AssertionError: GeoPackage driver not registered in GDAL/OGR

Diagnosis:

import pyogrio
print([d for d in pyogrio.list_drivers() if "gpkg" in d.lower()])
# Returns empty list

Fix: Reinstall pyogrio with a prebuilt binary that includes GeoPackage support. Preferred path:

conda install -c conda-forge pyogrio

If conda is unavailable: pip install pyogrio pulls precompiled GDAL wheels on most platforms. For source builds, compile GDAL with --with-sqlite3 --enable-shared. Check gdalinfo --formats | grep GPKG to confirm the driver is available at the system level before reinstalling the Python package.

2. Silent CRS Drift After Write #

Symptom: gpkg_spatial_ref_sys contains an entry but downstream QGIS or ArcGIS opens the layer in an unknown CRS.

Diagnosis:

gdf = gpd.read_file("processed_survey.gpkg", layer="validated_assets", engine="pyogrio")
print(gdf.crs)  # None or unexpected EPSG

-- GeoPackage context
SELECT srs_id, organization, organization_coordsys_id, definition
FROM gpkg_spatial_ref_sys
WHERE srs_id = (
  SELECT srs_id FROM gpkg_geometry_columns WHERE table_name = 'validated_assets'
);

Fix: Always set gdf.crs explicitly before calling to_file(). Never rely on CRS being inferred from geometries. If the source layer had None CRS, assign before any transformation:

if gdf.crs is None:
    gdf = gdf.set_crs("EPSG:4326")

3. Database Is Locked #

Symptom: sqlite3.OperationalError: database is locked during concurrent field-sync writes.

Diagnosis: Multiple Python processes (or a Python process plus QGIS) have the file open for writing simultaneously.

Fix: Enable WAL mode (see Performance Notes above), set a busy timeout, and serialize writes through a queue:

conn = sqlite3.connect("processed_survey.gpkg", timeout=30.0)
conn.execute("PRAGMA journal_mode=WAL;")

For persistent concurrent write workloads, implement an exclusive write lock at the application layer. The pooling patterns in Connection Pooling & Lifecycle Management show how to serialize writers while allowing concurrent readers through the WAL snapshot mechanism.

4. R-tree Out of Sync After Bulk Insert #

Symptom: Spatial queries return no results or wrong results after a large insert. ogrinfo reports correct feature count but QGIS renders the layer empty.

Diagnosis:

-- GeoPackage context: check for R-tree virtual table
SELECT name FROM sqlite_master
WHERE type = 'table' AND name LIKE 'rtree_validated_assets%';

Fix: Drop and recreate the spatial index:

-- GeoPackage context
DROP TABLE IF EXISTS rtree_validated_assets_geom;
-- Re-trigger GDAL index creation by opening and re-saving the layer,
-- or use the GPKG extension trigger mechanism directly.
SELECT gpkgAddSpatialIndex('validated_assets', 'geom');

The gpkgAddSpatialIndex() function is available in recent GDAL SQLite installations. Alternatively, re-write the layer through GeoPandas with SPATIAL_INDEX="YES".

5. Geometry Type Mismatch in gpkg_geometry_columns #

Symptom: WARNING 1: Geometry to be inserted is of type Multi Polygon, whereas the column type is Polygon. Written features silently downcast or are rejected.

Diagnosis:

-- GeoPackage context
SELECT geometry_type_name FROM gpkg_geometry_columns
WHERE table_name = 'validated_assets';

Fix: Pass geometry_type="MULTIPOLYGON" (or "GEOMETRY" for mixed-type layers) in the to_file() call, or normalize all geometries to a single type before writing:

# Explode any multi-part geometries to single-part before writing a POLYGON layer.
gdf = gdf.explode(index_parts=False).reset_index(drop=True)

Legacy Migration and Interoperability #

Many organizations still maintain shapefile archives that require modernization. Shapefiles lack transaction support, enforce 10-character field name limits, and scatter geometry and attributes across at least three separate files. Migrating to GeoPackage resolves these constraints while preserving attribute fidelity and enabling multi-layer consolidation in a single portable file.

The child page Converting Shapefiles to GeoPackage with GeoPandas covers schema normalization, field name mapping, CRS harmonization, and batch conversion patterns. During migration, always audit field name truncation and encoding shifts — particularly for UTF-8 metadata or extended ASCII characters that shapefiles store as Latin-1.

For teams maintaining hybrid environments where pyogrio encounters unsupported layer types (WebP raster tiles in GPKG, non-standard extensions), the Fiona & OGR Driver Configuration guide provides fallback driver registration, environment variable tuning, and CPL_DEBUG=ON tracing patterns.

Child Pages #

• Converting Shapefiles to GeoPackage with GeoPandas — automated schema normalization, CRS harmonization, and field name mapping for batch shapefile-to-GPKG migration pipelines.

Related #

• Python Integration & Database Workflows — parent overview: connection lifecycles, serialization patterns, and concurrency control for SpatiaLite and GeoPackage.
• Native sqlite3 Spatial Extensions — drop below the GDAL layer to run custom SQL predicates, rebuild spatial indexes, and call SpatiaLite geometry functions directly.
• Connection Pooling & Lifecycle Management — serialize multi-writer workloads, pre-warm connections with mod_spatialite already loaded, and avoid SQLITE_BUSY in field-sync scenarios.
• Transaction Scoping & Rollback Strategies — explicit BEGIN/COMMIT/ROLLBACK control for atomic multi-layer GeoPackage updates.
• Fiona & OGR Driver Configuration — driver registration, CPL_DEBUG tracing, and fallback I/O paths when pyogrio encounters non-standard GeoPackage extensions.