SpatiaLite & GeoPackage Automation with Python

How to Create a Spatial Index in SQLite with Python

From a sqlite3 connection with modspatialite loaded, call SELECT CreateSpatialIndex('table', 'geomcolumn') to build an R-tree index, then verify the…

From a sqlite3 connection with mod_spatialite loaded, call SELECT CreateSpatialIndex('table', 'geom_column') to build an R-tree index, then verify the planner uses it by running EXPLAIN QUERY PLAN on a bounding-box query and looking for a scan of the idx_<table>_<geom> virtual table.

This page sits under the Native sqlite3 & Spatial Extensions guide, which covers loading mod_spatialite and calling spatial SQL directly from the standard-library driver. Here the focus is narrow: building an R-tree spatial index from Python and proving that queries actually hit it.

Why This Matters

Without a spatial index, every bounding-box or nearest-feature query is a full table scan: SQLite computes the geometry predicate for all rows. On a 500,000-row parcel layer that turns a sub-second map redraw into a multi-second stall — unacceptable on field hardware. An R-tree index stores each geometry’s bounding box in a compact tree so the planner can discard almost every non-matching row before touching the geometry BLOB.

SpatiaLite and GeoPackage both build on SQLite’s R-tree module but wire it up differently. SpatiaLite exposes CreateSpatialIndex(), which creates an idx_<table>_<geom> virtual table and keeps it in sync with triggers. GeoPackage defines its own gpkg_rtree_index extension with tables named rtree_<table>_<geom>. Knowing which one your file uses — and confirming the planner picks it up — is the difference between an index that speeds queries and dead weight the optimizer ignores.

Prerequisites

  • Python 3.9+ with a sqlite3 build where enable_load_extension is available
  • mod_spatialite installed (apt install libsqlite3-mod-spatialite, brew install libspatialite, or the conda-forge librspatialite package)
  • A table with a registered geometry column (SpatiaLite metadata initialised, or a GeoPackage feature table)
  • Familiarity with loading mod_spatialite from sqlite3
  • GDAL/OGR 3.4+ only if you plan to cross-check with ogrinfo

Primary Method

python
# sqlite3 + mod_spatialite — build and verify an R-tree spatial index
import sqlite3


def create_spatial_index(db_path: str, table: str, geom_col: str = "geom") -> None:
    conn = sqlite3.connect(db_path)
    try:
        conn.enable_load_extension(True)
        conn.load_extension("mod_spatialite")

        # Ensure the spatial metadata exists (safe no-op if already initialised).
        # WGS84_ONLY keeps the srs_ref table small; use 'NONE' then insert your own.
        already = conn.execute(
            "SELECT count(*) FROM sqlite_master WHERE name='spatial_ref_sys'"
        ).fetchone()[0]
        if not already:
            conn.execute("SELECT InitSpatialMetaData(1)")

        # Build the R-tree. Returns 1 on success, 0 if it already exists.
        ok = conn.execute(
            "SELECT CreateSpatialIndex(?, ?)", (table, geom_col)
        ).fetchone()[0]
        conn.commit()

        if ok != 1:
            # 0 usually means the index already exists — recover it to be sure.
            conn.execute("SELECT RecoverSpatialIndex(?, ?)", (table, geom_col))
            conn.commit()
    finally:
        conn.close()

CreateSpatialIndex builds the idx_<table>_<geom> R-tree virtual table and installs triggers that keep it current as rows are inserted, updated, or deleted. It returns 1 on success and 0 when the index already exists or the geometry column is not registered.

Step-by-Step Walkthrough

1. Load the extension

The standard-library driver ships extension loading disabled. Enable it, then load mod_spatialite by its bare name so the OS resolver finds the shared object:

python
import sqlite3

conn = sqlite3.connect("parcels.gpkg")
conn.enable_load_extension(True)
conn.load_extension("mod_spatialite")
print(conn.execute("SELECT spatialite_version()").fetchone()[0])

2. Confirm the geometry column is registered

CreateSpatialIndex only works on a geometry column SpatiaLite knows about. Check the registry — if the row is missing, register it with RecoverGeometryColumn:

python
row = conn.execute(
    "SELECT srid, geometry_type FROM geometry_columns "
    "WHERE f_table_name = ? AND f_geometry_column = ?",
    ("parcels", "geom"),
).fetchone()
print(row)   # e.g. (4326, 3) -> POLYGON in EPSG:4326

3. Create the R-tree index

python
ok = conn.execute(
    "SELECT CreateSpatialIndex('parcels', 'geom')"
).fetchone()[0]
conn.commit()
print("created" if ok == 1 else "already existed / not registered")

4. Query through the index with a bounding-box predicate

An R-tree only helps when the query expresses a bounding-box filter. In SpatiaLite the idiom is a subquery against the idx_ virtual table:

sql
-- SpatiaLite: parcels intersecting a search window, using the R-tree
SELECT p.id, p.name
FROM parcels AS p
WHERE p.ROWID IN (
    SELECT ROWID FROM SpatialIndex
    WHERE f_table_name = 'parcels'
      AND search_frame = BuildMbr(-1.6, 50.9, -1.4, 51.1, 4326)
);

5. The GeoPackage variant

A GeoPackage feature table uses the OGC R-tree extension. Create it from Python against the rtree_<table>_<geom> naming:

sql
-- GeoPackage: create the OGC-defined R-tree on a feature table
SELECT gpkgAddSpatialIndex('parcels', 'geom');

If the gpkgAddSpatialIndex helper is unavailable in your mod_spatialite build, GDAL creates the same index when you write the layer with the SPATIAL_INDEX=YES layer creation option.

6. Populate or rebuild after a bulk load

The index triggers fire per-row, which slows large inserts. For a bulk load, it is faster to load first and build the index afterward — the same pattern covered in automating R-tree index rebuilds after a bulk load.

Verification

The only proof that matters is whether the planner uses the index. Run EXPLAIN QUERY PLAN and look for the R-tree virtual table in the plan:

python
plan = conn.execute(
    "EXPLAIN QUERY PLAN "
    "SELECT ROWID FROM SpatialIndex "
    "WHERE f_table_name = 'parcels' "
    "  AND search_frame = BuildMbr(-1.6, 50.9, -1.4, 51.1, 4326)"
).fetchall()
for step in plan:
    print(step[-1])

A healthy plan mentions idx_parcels_geom (SpatiaLite) or rtree_parcels_geom (GeoPackage). If instead you see SCAN parcels, the index is not being used — the predicate is not expressed as a bounding-box search. You can also confirm the virtual table exists:

sql
-- Either format: confirm the spatial index virtual table was created
SELECT name FROM sqlite_master
WHERE type = 'table'
  AND (name LIKE 'idx_parcels_%' OR name LIKE 'rtree_parcels_%');

Alternative Approaches or Edge Cases

GDAL at write time. If you build the GeoPackage with Fiona or ogr2ogr, passing the SPATIAL_INDEX=YES layer creation option creates the R-tree during the write, so no separate step is needed. This is the least error-prone path when you control the export.

Recovering a stale index. If rows were bulk-loaded with triggers disabled, or the index was created before the data, it can be out of sync. SELECT RecoverSpatialIndex('parcels', 'geom') rebuilds the R-tree from the current geometries without dropping and recreating the virtual table.

Troubleshooting

sqlite3.OperationalError: no such function: CreateSpatialIndex

Cause: mod_spatialite was not loaded on this connection, or extension loading is disabled. Fix: Call conn.enable_load_extension(True) then conn.load_extension("mod_spatialite") before the index SQL — see using sqlite3 with SpatiaLite functions for platform-specific library names.

CreateSpatialIndex returns 0 and no index appears

Cause: The geometry column is not registered in geometry_columns, so SpatiaLite refuses to index it. Fix: Run SELECT RecoverGeometryColumn('parcels', 'geom', 4326, 'POLYGON') to register it, then re-run CreateSpatialIndex.

EXPLAIN QUERY PLAN shows SCAN parcels despite the index existing

Cause: The query filters with a function like ST_Intersects directly instead of a bounding-box search_frame, so the planner cannot use the R-tree. Fix: Add a ROWID IN (SELECT ... FROM SpatialIndex WHERE search_frame = BuildMbr(...)) subquery to pre-filter by bounding box, then apply the exact predicate to the survivors.