Build and Use a Spatial Index in PyQGIS

Every spatial join written as a nested loop is quadratic. Two layers of ten thousand features means a hundred million exact predicate evaluations, each of which is real geometric work — the script does not crash, it simply appears to hang. A QgsSpatialIndex is an R-tree over the features' bounding boxes, and it turns that hundred million into ten thousand cheap box lookups followed by a handful of exact tests each.

This page is a focused recipe within Geometry Operations and Spatial Predicates in PyQGIS. It covers building an index, the two query methods, the memory trade-off of storing geometries inside it, and keeping the index correct when the underlying layer changes.

How an R-tree narrows the searchEight features are grouped into two bounding boxes. A query rectangle overlaps only the left group's box, so the search descends into that branch and tests its four leaves, skipping the right branch and its four features entirely without examining them.A branch whose box misses the query is skipped wholeroot nodebranch Abox overlaps the querybranch Bno overlap — skipped4 exact testsnever examinedquery rectangle

Prerequisites

  • QGIS 3.34 LTR (bundled Python 3.12) or newer.
  • Two layers in the same CRS — the index stores raw coordinates and does no transformation.
  • Enough memory for the index. A bare index is small; one storing geometries is roughly the size of the geometry data itself.

Build an index

The constructor takes a feature iterator and consumes it in one pass.

from qgis.core import QgsProject, QgsSpatialIndex

zones = QgsProject.instance().mapLayersByName("flood_zones")[0]

index = QgsSpatialIndex(zones.getFeatures())
print("indexed", zones.featureCount(), "feature(s)")

Breakdown: The iterator is consumed completely, so it cannot be reused afterwards — call getFeatures() again if you need a second pass. Building the index costs one full read of the layer, which is why it pays off only when you will query it more than a handful of times. For a single lookup, QgsFeatureRequest().setFilterRect() is simpler and cheaper.

By default the index stores only bounding boxes and feature IDs. That keeps it small but means every candidate requires a provider round-trip to fetch its geometry. For a tight loop, storing the geometries inside the index is usually the better trade:

index = QgsSpatialIndex(
    zones.getFeatures(), flags=QgsSpatialIndex.FlagStoreFeatureGeometries
)
geometry = index.geometry(some_feature_id)

Breakdown: FlagStoreFeatureGeometries keeps a copy of every geometry, so index.geometry(id) returns it with no provider access at all. The memory cost is proportional to the geometry data — trivial for points, significant for a national polygon coverage. The alternative, building a {fid: geometry} dictionary yourself, uses the same memory with more code, so prefer the flag.

Query by rectangle

intersects() on the index takes a QgsRectangle and returns candidate feature IDs whose boxes overlap it.

from qgis.core import QgsProject, QgsSpatialIndex

parcels = QgsProject.instance().mapLayersByName("parcels")[0]
zones = QgsProject.instance().mapLayersByName("flood_zones")[0]

index = QgsSpatialIndex(
    zones.getFeatures(), flags=QgsSpatialIndex.FlagStoreFeatureGeometries
)

matches = {}
for parcel in parcels.getFeatures():
    geometry = parcel.geometry()
    hits = [
        zone_id
        for zone_id in index.intersects(geometry.boundingBox())
        if geometry.intersects(index.geometry(zone_id))
    ]
    if hits:
        matches[parcel["parcel_id"]] = hits

print(len(matches), "parcel(s) matched at least one zone")

Breakdown: The two-stage structure is the whole point and must not be collapsed. index.intersects() is an approximation — it returns everything whose box overlaps, including features the exact geometry misses entirely — so the inner geometry.intersects() is what produces the correct answer. Skipping it gives wrong results that look plausible, which is the worst kind. Note the index is built on the smaller layer where possible: the layer you index is the one you look things up in, and the layer you iterate is the one you scan.

Query by proximity

nearestNeighbor() answers "what is closest to this" without any distance computation on your side.

schools = QgsProject.instance().mapLayersByName("schools")[0]

for school in schools.getFeatures():
    geometry = school.geometry()
    candidates = index.nearestNeighbor(geometry, 5)
    nearest = min(
        candidates, key=lambda fid: geometry.distance(index.geometry(fid)), default=None
    )
    if nearest is not None:
        print(school["name"], "->", nearest)

Breakdown: nearestNeighbor(geometry, 5) ranks by bounding-box proximity, which can mis-order long or irregular shapes — a long diagonal river has a huge box whose corner may be nearer than a compact feature that is genuinely closer. Requesting several candidates and picking the true minimum with distance() corrects that. min(..., default=None) handles the empty-index case without an extra branch. For a whole-layer result, native:joinbynearest does the same thing and writes an output layer — see Calculate the Distance Between Features.

Keep the index in step with edits

An index is a snapshot. Add, move or delete a feature in the layer and the index knows nothing about it, so queries silently return stale results.

Three ways to keep an index correctA layer feeds an index. When the layer is edited the index becomes stale. Three remedies are shown: calling addFeature and deleteFeature to update it incrementally, rebuilding the whole index from scratch, or deferring index construction until after all edits are committed, which is labelled as the simplest option.An index does not watch the layer — you keep them in steplayer editedindex now staleindex.addFeature / deleteFeaturecheap, but every edit must be mirroredrebuild from getFeatures()simple, costs one full readbuild after commitChanges()no synchronisation problem at allprefer thisread-only index, no drift

from qgis.core import QgsSpatialIndex

index = QgsSpatialIndex(zones.getFeatures())

new_feature = make_zone()                     # your own construction
zones.dataProvider().addFeatures([new_feature])
index.addFeature(new_feature)                 # mirror the change

zones.dataProvider().deleteFeatures([old_id])
index.deleteFeature(old_zone_feature)         # takes the feature, not the id

Breakdown: addFeature() and deleteFeature() update the tree incrementally, which is far cheaper than rebuilding. The asymmetry catches people out: deleteFeature() needs the whole feature — including its geometry, so the tree can find the right node — not just the ID, which means you have to hold on to it before deleting from the provider. In practice the simplest correct approach is to finish all edits, commit, and only then build the index; a read-only index cannot drift.

When an index is not the answer

Building an index costs a full read of the layer. That is only worth it when you will query it many times, and there are two situations where something else is simply better.

Choosing the right tool for a spatial lookupStarting from a spatial lookup, three branches. If only a handful of queries are needed, a single feature request with a filter rectangle wins. If many queries are needed against a file-based layer, a spatial index wins. If both layers live in the same PostGIS database, a SQL spatial join beats anything done in Python.How many lookups, and where does the data live?a spatial lookupfewer than ~100 queriessetFilterRect() on a requestmany queries, file layerQgsSpatialIndexboth layers in PostGISa SQL spatial joinno build cost to amortiseone read, then thousands of lookupsthe database already has an index

For a handful of lookups, skip the index entirely:

from qgis.core import QgsFeatureRequest

request = QgsFeatureRequest().setFilterRect(query_geometry.boundingBox())
candidates = [
    f for f in zones.getFeatures(request)
    if query_geometry.intersects(f.geometry())
]

Breakdown: setFilterRect() pushes the box test to the provider, which for a GeoPackage or PostGIS layer uses the file or database's own spatial index — one that already exists and cost you nothing to build. For a small number of queries this beats building a QgsSpatialIndex outright, and the two-stage exact test is identical.

When both layers live in the same PostGIS database, neither approach is right: the join belongs in SQL, where the planner can use both tables' indexes and never transfer a row into Python. Load the result as a layer with a query-backed URI rather than iterating in PyQGIS at all — the same principle as pushing filters down described in Working with QGIS Expressions.

QGIS version compatibility

The examples target QGIS 3.34 LTR (Python 3.12).

QGIS versionPythonNotes
3.28 LTR3.9Identical API. FlagStoreFeatureGeometries available since 3.4.
3.34 LTR3.12Baseline for this page.
3.40 / 3.443.12QgsSpatialIndexKDBush added for point-only layers — faster and smaller where it applies.

QgsSpatialIndex, intersects() and nearestNeighbor() are unchanged across 3.x. On a points-only layer in 3.40 or newer, QgsSpatialIndexKDBush is worth benchmarking.

Troubleshooting

  • Results include features that clearly do not match. The exact predicate step was skipped. index.intersects() returns box candidates, not answers.
  • Results are missing features that should match. The two layers are in different CRSs, or the index was built before the features were added.
  • The index appears empty. The feature iterator was already consumed. getFeatures() returns a single-use iterator — call it fresh for the constructor.
  • index.geometry(id) raises or returns nothing. The index was built without FlagStoreFeatureGeometries. Either add the flag or keep your own {fid: geometry} dictionary.
  • Building the index takes longer than the query saved. You are querying too few times. Below roughly a hundred lookups, a plain setFilterRect() request is faster overall.
  • Nearest neighbour picks something visibly further away. Only one candidate was requested. Ask for five and confirm with an exact distance().

Conclusion

A spatial index converts a quadratic scan into a tree descent plus a few exact tests. Build it on the layer you look things up in, add FlagStoreFeatureGeometries when you will fetch geometries in a loop, always follow a box query with the exact predicate, and build the index after edits are committed rather than trying to keep a live one in step.

Frequently Asked Questions

Does the index give me the exact answer? No. It returns candidates whose bounding boxes match, which is an approximation containing false positives. Always follow it with the exact predicate — intersects(), contains() or distance() — on each candidate.

Should I store geometries in the index? Yes, if you will fetch candidate geometries inside a loop. FlagStoreFeatureGeometries avoids a provider round-trip per candidate at the cost of memory proportional to the geometry data.

Which layer should I index? The one you look things up in, not the one you iterate. If you are asking "which zone does each parcel fall in", index the zones and iterate the parcels.

What happens to the index when I edit the layer? Nothing — it silently goes stale. Mirror each change with addFeature() and deleteFeature(), or simply build the index after committing your edits so it never needs synchronising.

Does the index work across coordinate reference systems? No. It stores raw coordinates, so a query geometry in a different CRS matches nothing. Reproject one side before building or querying the index.