How to Create a Spatial Index in SQLite with Python

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.

The saving comes from pruning, not from faster comparisons. An R-tree groups nearby features under a shared bounding rectangle, and groups those rectangles under larger ones, all the way to a single root. A bounding-box query descends the tree and discards any branch whose rectangle does not intersect the search window — and discarding a branch discards everything beneath it without ever reading a geometry.

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 #

# 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:

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:

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 #

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:

-- 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)
);

Note the shape of that query: it is two filters, not one. The R-tree subquery is cheap and approximate — it returns everything whose bounding box overlaps the window, including features whose actual geometry does not. The exact predicate then runs only on those survivors. Skipping the first filter means running the expensive test on every row; skipping the second returns false positives.

5. The GeoPackage variant #

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

-- 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.

The two containers use the same underlying R-tree module, so the performance is identical — but almost every name around it differs, and code that hard-codes one set of names silently does nothing on the other format:

Detecting which wiring a file uses is a single query against sqlite_master: a container with a gpkg_contents table is a GeoPackage and wants the rtree_ path, while one with geometry_columns and spatial_ref_sys in the SpatiaLite layout wants CreateSpatialIndex. Files that carry both — a GeoPackage that has had InitSpatialMetaData run against it — are legal but ambiguous, and the GeoPackage path should win, because that is what other readers will look for.

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:

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:

-- 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.

Frequently Asked Questions #

Related #

• Native sqlite3 & Spatial Extensions — parent guide: loading and calling spatial SQL from the standard-library driver
• Using sqlite3 with SpatiaLite Functions — extension loading and platform library names
• How to Automate R-tree Index Rebuilds After Bulk Load — rebuild the index efficiently after large inserts
• GeoPackage Specification Deep Dive — how the OGC R-tree extension is defined and named