Spatial Indexing and Query Optimization in PyQGIS

When processing municipal parcel datasets, environmental monitoring grids, or large infrastructure networks, naive feature iteration becomes a hard computational bottleneck — $O(n)$ per query means milliseconds become minutes at scale. This page, part of the PyQGIS Core Architecture & Data Handling guide, covers the QgsSpatialIndex R-tree, two-phase spatial query execution, QgsFeatureRequest filter composition, and advanced patterns for nearest-neighbour searches and hybrid attribute-spatial filtering. By the end, you will have production-ready code that reduces query latency from minutes to milliseconds regardless of dataset size.

Prerequisites Checklist #

Before implementing spatial indexing patterns, confirm your environment meets these baseline requirements:

• QGIS 3.28+ (LTR recommended): API stability and optimised C++ backend bindings. QgsSpatialIndex.FlagStoreFeatureGeometries is available from 3.20 onward.
• Python 3.8+: Required for the type hints used in every code sample below.
• Core API familiarity: Working knowledge of QgsVectorLayer, QgsFeature, and QgsGeometry is assumed; see vector and raster data access patterns if you need a refresher.
• Validation datasets: Medium-to-large vector layers (100 000+ features) are needed to observe meaningful performance differences; use OS OpenData boundaries or Natural Earth admin layers if you lack production data.
• CRS alignment: All geometries passed to the index must share the same coordinate reference system. If your pipeline mixes projections, resolve that first — the coordinate transformations and CRS handling guide details how to normalise CRS across layers before any spatial operation.

Spatial indexing is a runtime acceleration layer, not a substitute for clean data modelling. It performs best when paired with valid geometries, consistent attribute schemas, and a well-defined CRS.

How the QGIS R-tree Index Works #

QgsSpatialIndex wraps libspatialindex’s R-tree implementation behind the PyQGIS SIP bindings. The R-tree partitions geometric space into a hierarchy of minimum bounding rectangles (MBRs). Each leaf node stores the MBR and feature ID of one geometry; internal nodes store the union of their children’s MBRs. A query descends the tree, pruning subtrees whose MBR cannot overlap the query envelope, and returns the feature IDs of leaf nodes whose MBR intersects.

This structure delivers $O(\log n)$ candidate retrieval instead of $O(n)$ full-table iteration. The key trade-off: the index answers bounding-box overlap, not exact geometry predicates. A narrow sliver polygon and a large rotated rectangle may share overlapping MBRs even when their actual geometries are disjoint. That is why every correct spatial index workflow has two phases: index-based candidate retrieval followed by exact geometry validation.

The diagram below shows the full data flow from a query geometry through the index to the validated result set.

The index itself is entirely in-memory and does not automatically track edits, feature additions, or deletions on the underlying QgsVectorLayer. For static or snapshot-based workflows, build a QgsSpatialIndex once per session. For interactive editing sessions, either rely on the layer’s built-in spatial index (activated by layer.dataProvider().createSpatialIndex()) or rebuild your custom index after each commit cycle. The official QgsSpatialIndex API documentation details every constructor flag and memory allocation behaviour.

Step-by-Step Implementation #

Step 1 — Verify the Layer and Align CRS #

Always validate the layer and confirm CRS alignment before touching the index. A silent CRS mismatch produces bounding boxes in different coordinate spaces — the index will return wrong candidates without raising any exception.

from qgis.core import (
    QgsVectorLayer,
    QgsSpatialIndex,
    QgsFeatureRequest,
    QgsGeometry,
    QgsCoordinateTransform,
    QgsProject,
)

def _assert_same_crs(layer: QgsVectorLayer, query_geom: QgsGeometry, query_crs_authid: str) -> None:
    """
    Raises ValueError if the layer CRS does not match the query CRS.
    Call this before index construction to catch projection mismatches early.
    """
    layer_crs = layer.crs().authid()
    if layer_crs != query_crs_authid:
        raise ValueError(
            f"CRS mismatch: layer is {layer_crs}, query geometry is {query_crs_authid}. "
            "Reproject one of them before building the spatial index."
        )

Step 2 — Build and Populate the Index #

Construct the index in a single feature pass. Using setNoAttributes() strips attribute data from the iteration, cutting RAM consumption by 40–60 % on wide tables.

def build_spatial_index(layer: QgsVectorLayer) -> QgsSpatialIndex:
    """
    Build a QgsSpatialIndex from a QgsVectorLayer in one memory-efficient pass.

    Args:
        layer: A valid QgsVectorLayer whose features have geometries.

    Returns:
        A populated QgsSpatialIndex ready for bounding-box queries.

    Raises:
        ValueError: If the layer is invalid or contains no features.
    """
    if not layer.isValid():
        raise ValueError(f"Layer '{layer.name()}' is invalid and cannot be indexed.")

    index = QgsSpatialIndex()
    # setNoAttributes() avoids fetching attribute data — geometry only is needed
    request = QgsFeatureRequest().setNoAttributes()

    count = 0
    for feature in layer.getFeatures(request):
        # Skip null or empty geometries — they corrupt the R-tree leaf list
        if feature.hasGeometry() and not feature.geometry().isEmpty():
            index.addFeature(feature)
            count += 1

    if count == 0:
        raise ValueError(f"Layer '{layer.name()}' produced zero indexable features.")

    return index

Key points:

1. setNoAttributes() is not optional for large datasets — skip it and you load every attribute value into Python memory during iteration.
2. hasGeometry() guards against null-geometry features that some PostGIS tables legitimately contain.
3. For layers exceeding 500 000 features, use QgsSpatialIndex(layer.getFeatures(request)) — the constructor overload accepts an iterator and inserts in bulk, bypassing per-feature Python call overhead.

Step 3 — Execute the Two-Phase Query #

Query the index for candidates, then validate with exact geometry predicates.

def query_exact_intersections(
    index: QgsSpatialIndex,
    layer: QgsVectorLayer,
    query_geom: QgsGeometry,
) -> list:
    """
    Return features from `layer` that exactly intersect `query_geom`.

    Phase 1: index.intersects() — O(log n) bounding-box candidate retrieval.
    Phase 2: geometry.intersects() — exact DE-9IM predicate to discard false positives.

    Args:
        index:      A QgsSpatialIndex built from `layer`.
        layer:      The source QgsVectorLayer (same CRS as query_geom).
        query_geom: The geometry to test against.

    Returns:
        List of QgsFeature objects that pass the exact intersection test.
    """
    # Phase 1: bounding-box overlap — fast, returns a superset
    candidate_ids = index.intersects(query_geom.boundingBox())
    if not candidate_ids:
        return []

    # Phase 2: exact geometric validation — eliminates false positives
    request = QgsFeatureRequest().setFilterFids(candidate_ids)
    return [
        feature
        for feature in layer.getFeatures(request)
        if feature.geometry().intersects(query_geom)
    ]

Step 4 — Wire Into a Plugin or Script #

The following shows how to use the index as a session-level cache within a QGIS plugin or headless script, avoiding redundant index rebuilds on repeated calls.

class SpatialQueryEngine:
    """
    Session-scoped engine that builds a QgsSpatialIndex once and reuses it.
    Invalidate and rebuild by calling reset() after any write transaction.

    See also:
        - Layer lifecycle: /pyqgis-core-architecture-data-handling/working-with-qgsproject-and-layer-registry/
        - Memory ownership: /pyqgis-core-architecture-data-handling/memory-management-and-garbage-collection-for-gis-objects/
    """

    def __init__(self, layer: QgsVectorLayer) -> None:
        self._layer = layer
        self._index: QgsSpatialIndex | None = None

    def _ensure_index(self) -> QgsSpatialIndex:
        if self._index is None:
            self._index = build_spatial_index(self._layer)
        return self._index

    def reset(self) -> None:
        """Drop the cached index. Call after any edit session that alters geometry."""
        self._index = None

    def intersects(self, query_geom: QgsGeometry) -> list:
        """Return features exactly intersecting query_geom."""
        return query_exact_intersections(self._ensure_index(), self._layer, query_geom)

    def nearest(self, point: "QgsPointXY", n: int = 1) -> list:
        """Return the n nearest features to point by centroid distance."""
        from qgis.core import QgsPointXY
        idx = self._ensure_index()
        candidate_ids = idx.nearestNeighbor(point, n)
        request = QgsFeatureRequest().setFilterFids(candidate_ids)
        return list(self._layer.getFeatures(request))

Advanced Patterns #

Attribute + Spatial Hybrid Filtering #

Combining spatial and attribute filters at the QgsFeatureRequest level lets the QGIS query engine evaluate both conditions inside the C++ data provider before loading features into Python. This is significantly faster than filtering in Python after the fact.

from qgis.core import QgsFeatureRequest, QgsGeometry, QgsVectorLayer

def active_parcels_within(layer: QgsVectorLayer, query_geom: QgsGeometry) -> list:
    """
    Return 'ACTIVE' parcels larger than 1 000 m² that intersect query_geom.
    The spatial filter executes at the provider level; the attribute expression
    runs on the already-reduced candidate set.
    """
    request = (
        QgsFeatureRequest()
        .setFilterRect(query_geom.boundingBox())
        .setFilterExpression("status = 'ACTIVE' AND area > 1000")
        .setNoAttributes()  # drop attributes we don't need in Python
    )
    return [
        f for f in layer.getFeatures(request)
        if f.geometry().intersects(query_geom)
    ]

setFilterRect combined with setFilterExpression is the most efficient path for single-pass filtering when you do not need nearest-neighbour semantics. Note that the QGIS engine evaluates the rect filter first, then the expression — so always structure expressions to reduce the candidate set as aggressively as possible.

Memory-Efficient Streaming with Generators #

When processing millions of features, avoid materialising a full result list. Use a generator to stream results through a processing pipeline one feature at a time. This keeps peak memory consumption flat regardless of result-set size, which is critical when paired with memory management for GIS objects.

from collections.abc import Generator
from qgis.core import QgsSpatialIndex, QgsVectorLayer, QgsGeometry, QgsFeature, QgsFeatureRequest

def stream_intersections(
    index: QgsSpatialIndex,
    layer: QgsVectorLayer,
    query_geom: QgsGeometry,
) -> Generator[QgsFeature, None, None]:
    """
    Yield features that exactly intersect query_geom without loading the full result set.
    Suitable for processing pipelines that write results incrementally (e.g. to GeoPackage).
    """
    for fid in index.intersects(query_geom.boundingBox()):
        feature = layer.getFeature(fid)
        if feature.hasGeometry() and feature.geometry().intersects(query_geom):
            yield feature

Nearest-Neighbour Queries #

QgsSpatialIndex.nearestNeighbor() returns a ranked list of feature IDs by bounding-box centroid distance. It is the only scalable approach for proximity analysis — iterating and computing QgsGeometry.distance() for every candidate is $O(n)$ and unacceptable beyond a few thousand features.

from qgis.core import QgsSpatialIndex, QgsVectorLayer, QgsPointXY, QgsFeatureRequest

def find_nearest_features(
    index: QgsSpatialIndex,
    layer: QgsVectorLayer,
    point: QgsPointXY,
    n: int = 5,
) -> list:
    """
    Return the n features geometrically closest to point.
    Distances are planar (units match the layer CRS); convert to ellipsoidal
    using QgsDistanceArea if absolute metric distance is required.
    """
    candidate_ids = index.nearestNeighbor(point, n)
    request = QgsFeatureRequest().setFilterFids(candidate_ids)
    return list(layer.getFeatures(request))

Important caveat: nearestNeighbor ranks by MBR-centroid distance, not exact geometry distance. For large or concave features, the ranking can differ from the true geometric nearest-neighbour order. Post-sort results by feature.geometry().distance(QgsGeometry.fromPointXY(point)) when exact ranking matters.

Pitfalls and Debugging #

• Stale index after edits: QgsSpatialIndex does not observe QgsVectorLayer edit signals. If an external ETL process or direct database write alters geometry after the index is built, all subsequent queries return stale candidates. Implement a rebuild trigger tied to layer.editingStopped or rebuild unconditionally before each processing batch.

• Invalid geometries entering the index: Self-intersecting polygons or zero-length lines produce malformed MBRs that corrupt the R-tree’s node-splitting logic. Run processing.run("native:fixgeometries", {"INPUT": layer, "OUTPUT": "memory:"}) or QgsGeometryValidator upstream before index construction. Check QgsGeometry.isGeosValid() per feature if you cannot control the input data.

• Projection drift (silent mismatch): Mixing a WGS 84 query geometry with a layer in EPSG:27700 produces bounding boxes in different numeric ranges — the index happily returns no candidates or all candidates, neither of which raises an exception. Always call layer.crs().authid() and compare it to your query geometry’s CRS before any index operation. Use QgsCoordinateTransform to align them, as described in coordinate transformations and CRS handling.

• High false-positive ratio: A candidate-to-result ratio above 10:1 signals poor index selectivity. Common causes: large bounding boxes on elongated or highly concave features, or a query geometry whose MBR is much larger than its actual area (e.g. a thin diagonal line). Mitigate by breaking the query geometry into smaller sub-geometries, or switch to setFilterExpression with a spatial function for pre-filtering.

• Layer invalid during index population: Always call layer.isValid() before iterating. Invalid layers return empty feature iterators without raising Python exceptions — your index silently builds empty, and all subsequent queries return no results.

• Thread safety: QgsSpatialIndex is not thread-safe for concurrent writes. If you populate the index on a background thread (e.g. inside QgsTask), complete population before any query thread accesses it. Use a threading lock or rebuild the index on the main thread after the background population task completes. The signal and slot event handling guide covers safe inter-thread communication patterns in QGIS.

Conclusion #

QgsSpatialIndex delivers $O(\log n)$ candidate retrieval and reduces spatial query latency from minutes to milliseconds at the cost of a one-time build pass and a validation step to discard false positives. Combining it with QgsFeatureRequest filter composition — spatial rect, attribute expression, and field subset — gives you the most efficient query path available in PyQGIS. Always align CRS before index construction, validate geometries upstream, and invalidate the index after any write cycle. With these patterns in place, your spatial queries remain responsive and memory-efficient as dataset complexity scales.

Related

• PyQGIS Core Architecture & Data Handling — parent guide covering the full data-handling stack
• Coordinate Transformations and CRS Handling — aligning CRS before index construction
• Vector and Raster Data Access Patterns — efficient feature iteration foundations
• Memory Management and Garbage Collection for GIS Objects — controlling peak memory in large spatial pipelines
• Working with QgsProject and Layer Registry — safe layer loading and lifecycle management before index population