PyQGIS Engineering & Automation

Nearest-Neighbour Search with QgsSpatialIndex

Find the genuinely nearest feature in PyQGIS: why nearestNeighbor ranks bounding box centroids rather than shapes, how many candidates to request, verifying…

TL;DR: QgsSpatialIndex.nearestNeighbor() ranks candidates by the distance between bounding box centroids, not between geometries — so ask for several candidates and measure the true distance yourself before deciding which is nearest. This page is part of the spatial indexing and query optimization guide.

Complete Runnable Code

python
"""Find the genuinely nearest feature, not merely the nearest centroid."""
from qgis.core import (QgsFeature, QgsFeatureRequest, QgsGeometry,
                       QgsSpatialIndex, QgsVectorLayer)


def build_index(layer: QgsVectorLayer) -> QgsSpatialIndex:
    """Populate an index in one pass, without fetching attributes."""
    if not layer.isValid():
        raise ValueError("layer %r is not valid" % layer.name())
    request = QgsFeatureRequest().setNoAttributes()
    return QgsSpatialIndex(layer.getFeatures(request))


def nearest_feature(layer: QgsVectorLayer, index: QgsSpatialIndex,
                    target: QgsGeometry, candidates: int = 10) -> tuple[QgsFeature, float] | None:
    """Return the nearest feature to `target` and its true distance.

    The index ranks by envelope centroid, which is only an approximation of
    shape distance, so `candidates` ids are fetched and measured properly.
    Returns None when the layer has no features.
    """
    ids = index.nearestNeighbor(target, candidates)
    if not ids:
        return None

    request = QgsFeatureRequest().setFilterFids(ids)
    best: QgsFeature | None = None
    best_distance = float("inf")
    for feature in layer.getFeatures(request):
        distance = feature.geometry().distance(target)      # true geometric distance
        if distance < best_distance:
            best, best_distance = QgsFeature(feature), distance

    return (best, best_distance) if best is not None else None
The two phases of a correct nearest-neighbour query A four-step pipeline: ask the index for k candidates ranked by centroid distance, fetch those features by id, measure true geometry distance, and take the closest — because the index ranks envelopes rather than shapes. nearestNeighbor(k) centroid ranking setFilterFids fetch candidates geometry.distance() true distance min() the real nearest ids features measured

Architecture Breakdown

What the index actually ranks

nearestNeighbor() walks the R-tree and returns feature ids ordered by the distance between the query point and the centroid of each candidate’s bounding box. For compact, similarly sized features that ordering matches shape distance closely. For anything long, thin or irregular it does not.

A river polygon spanning a county has a centroid tens of kilometres from most of its bank. A point standing on that bank is much closer to the river than to a small pond ten metres away, and the index will rank the pond first every time.

Why you ask for more than one

Because the correction is cheap and the error is not detectable otherwise. Fetching k candidates by id is a single indexed request, and measuring true distance on a handful of geometries costs almost nothing next to the query that produced them.

How many candidates you need to ask for A bar chart of how often the true nearest feature is missed when asking the index for 1, 3, 10 and 25 candidates, over a dataset of elongated polygons. k = 1 18% wrong k = 3 6% wrong k = 10 1.2% wrong k = 25 0.2% wrong Indicative miss rates over a layer of elongated polygons, where a centroid is a poor proxy for the shape. On compact features k=3 is usually ample.

How large k needs to be depends entirely on how well a centroid represents your features. Compact parcels need very little; elongated or multipart features need considerably more. Measuring the miss rate on a sample of your own data is a twenty-line script and settles the question for that dataset.

CRS and units

distance() returns a value in the layer’s coordinate reference system units. On a geographic CRS that is degrees, which is not a distance anyone can use and varies with latitude. Every proximity workflow should either run in a projected CRS or convert deliberately — the coordinate transformations guide covers doing that once rather than per feature.

Choosing the Right Tool for the Question

Proximity covers several genuinely different questions, and only one of them needs an index at all.

Which proximity tool the question needs A decision tree over proximity questions: a radius query is a rectangle filter plus a distance test, a single nearest is a k-nearest query with verification, and an all-pairs join belongs in a Processing algorithm. What is the proximity question? within X metres rect + distance no index needed the nearest one nearestNeighbor then verify nearest for every row native:joinbynearest it is optimised

Not every proximity question needs an index. “Everything within 500 metres” is a rectangle filter followed by an exact distance test, which the provider can answer directly and which stays correct as the data changes. Building an index for that is work you do not need.

The index earns its place for genuine nearest-neighbour ranking, which no request filter can express, and for repeated queries against unchanging data. For the all-pairs case — the nearest feature for every row of another layer — reach for native:joinbynearest, which does the whole join in optimised C++ and handles the verification internally.

Using It Inside a Plugin

python
class NearestFinder:
    """Caches one index per layer for the life of a plugin session."""

    def __init__(self):
        self._indexes: dict[str, QgsSpatialIndex] = {}

    def index_for(self, layer: QgsVectorLayer) -> QgsSpatialIndex:
        key = layer.id()
        if key not in self._indexes:
            self._indexes[key] = build_index(layer)
            layer.editingStopped.connect(lambda lid=key: self._indexes.pop(lid, None))
        return self._indexes[key]

The connection is the important half. An index is a snapshot, so caching one without invalidating it on edit gives you fast answers that are quietly wrong — arguably worse than no cache at all.

Keeping the Index Honest Over Time

The failure that matters in production is not a wrong ranking but a stale index. Because the structure is a snapshot of envelopes taken at build time, every edit to the underlying layer makes it a little less true, and nothing about the API will tell you so.

There are two workable disciplines. In a batch job, build the index immediately before the query phase and discard it afterwards, so its lifetime is shorter than the interval between any two writes. In an interactive plugin, cache it but invalidate on the layer’s editing signals, as the class above does.

What does not work is rebuilding on a timer or on a count of queries. Both are approximations of “has the data changed”, and the layer already answers that question exactly.

Production Best Practices

  • Always verify with geometry().distance(). The index ranks envelopes.
  • Ask for enough candidates. Measure the miss rate on your data rather than guessing.
  • Work in a projected CRS so distances are in metres.
  • Rebuild the index after edits, or connect to editingStopped and drop it.
  • Use setNoAttributes() when building, which cuts index construction memory substantially.
  • Prefer native:joinbynearest for all-pairs work; it is faster and already correct.

Frequently Asked Questions

Can I get the distance from the index directly?

No — nearestNeighbor() returns ids only. That is deliberate: the ordering is based on an approximation, so returning a number would invite treating it as the answer. Fetch the candidates and measure.

What does FlagStoreFeatureGeometries change?

It makes the index keep a copy of each geometry, so geometry() can be answered from the index without going back to the layer. That saves a fetch at the cost of holding every geometry in memory twice.

It is worth it when you query the same index constantly and the layer is expensive to read; it is not worth it in the standard two-phase pattern, where fetching by id is already cheap.

How many candidates should I request?

Enough that the true nearest is essentially always in the set. Three is plenty for compact, evenly sized features; ten is a reasonable default for mixed data; anything long and thin may need several dozen. The cost of a larger k is one slightly bigger fetch, so erring high is inexpensive.

Does it work with points, lines and polygons together?

Yes — the index stores envelopes and does not care about geometry type, and distance() is defined between any two geometries. The centroid approximation simply gets worse as the features get less compact, which is the argument for a larger k rather than for a different tool.

Is the index thread-safe?

It is safe to read from several threads once population is complete, and not safe to write to while anything reads. Build it fully before sharing it, and never insert into a shared index from a worker thread.

What if two features are equidistant?

The loop as written keeps the first one it measured, and the order of getFeatures() on a filtered id set is not guaranteed. If ties matter — and for anything auditable they do — sort the candidates by distance and then by feature id, so the same input always produces the same answer.