Join Attributes by Nearest Neighbour in PyQGIS

A spatial join by location needs the features to touch. Most real questions do not come with touching geometry: which fire station is nearest each school, which weather station should each field use, which bus stop is closest to each new housing site, what is the nearest hydrant to every building. Those are nearest-neighbour joins, and they have their own traps — a nearest neighbour always exists unless you cap the distance, ties are resolved arbitrarily, and "nearest" in degrees is not nearest on the ground.

This recipe belongs to Geometry Operations and Spatial Predicates in PyQGIS. It runs the Processing algorithm, asks for several neighbours with a maximum distance, reproduces the join in Python with a spatial index for full control, and covers the checks that make the results defensible.

Nearest is not always nearThree schools as blue squares and fire stations as red triangles. The first school has one station 1.2 km away and is joined to it with distance recorded. The second school has two stations at 2.4 and 2.5 km; the join picks one, and asking for two neighbours reveals the near tie. The third school's nearest station is 31 km away, beyond a 10 km cap, so it is left unmatched rather than joined to something meaningless.Every input finds a neighbour — unless you set a limitclear matchstation 7 · 1.2 kmone obvious answernear tie2.4 km vs 2.5 kmask for k = 2 to see itbeyond the cap31 km > 10 km capleft unmatched

Prerequisites

  • QGIS 3.40 LTR or newer, or the QGIS 4 series.
  • Two layers in the same projected CRS. The algorithm measures in layer units, so degrees give degree distances and a nearest neighbour that is only nearest on a flat map.
  • A clear idea of the largest distance at which a match still means something for your question.

Run the nearest-neighbour join

native:joinbynearest copies attributes from the nearest feature in a second layer onto each input feature and adds the distance.

import processing
from qgis.core import QgsProject, QgsVectorLayer

schools = QgsProject.instance().mapLayersByName("schools")[0]
stations = QgsProject.instance().mapLayersByName("fire_stations")[0]
assert schools.crs() == stations.crs() and not schools.crs().isGeographic()

result = processing.run("native:joinbynearest", {
    "INPUT": schools,
    "INPUT_2": stations,
    "FIELDS_TO_COPY": ["station_id", "station_name", "crew_type"],
    "DISCARD_NONMATCHING": False,
    "PREFIX": "fs_",
    "NEIGHBORS": 1,
    "MAX_DISTANCE": 10000,
    "OUTPUT": "/data/work/schools_nearest_station.gpkg",
})
joined = QgsVectorLayer(result["OUTPUT"], "schools + nearest station", "ogr")
print(result["JOINED_COUNT"], "joined,", result["UNJOINABLE_COUNT"], "without a station in range")

Breakdown: FIELDS_TO_COPY keeps the output narrow; left empty, every field of the station layer is copied. PREFIX avoids collisions when both layers have a name field. MAX_DISTANCE is in layer units — 10,000 metres here — and is the most important parameter: without it every school is joined to some station, however far. With DISCARD_NONMATCHING false, schools with no station in range stay in the output with empty join fields, which is almost always what you want for a coverage analysis, because the unmatched ones are the finding. Besides the copied fields, the algorithm adds n (the neighbour rank), distance, and the coordinates of both the input and the matched feature.

The assertion at the top is cheap insurance. On a layer in EPSG:4326, a MAX_DISTANCE of 10,000 means ten thousand degrees — everything matches — and a degree of longitude shrinks towards the poles, so the nearest by degrees can be the wrong station.

Several neighbours, and ties

Asking for more than one neighbour answers different questions — the second-nearest station as a backup, the three nearest weather stations for an average — and it exposes ties the single-neighbour join hides.

k neighbours means k rows per inputAn output table for two schools with NEIGHBORS equal to 3 and a 10 km cap. School A appears three times with n equal to 1, 2 and 3 at distances 1.2, 4.8 and 7.9 km. School B appears twice, with n equal to 1 and 2 at 2.4 and 2.5 km, because its third nearest station is beyond the cap. The near-equal first and second distances for school B flag a tie worth reviewing.NEIGHBORS = 3 multiplies rows, not columnsschoolnfs_station_iddistanceA1FS071,210A2FS034,780A3FS117,905B1FS022,412B2FS092,498near tie

from qgis.core import QgsVectorLayer, QgsFeatureRequest

three = processing.run("native:joinbynearest", {
    "INPUT": schools, "INPUT_2": stations,
    "FIELDS_TO_COPY": ["station_id"], "PREFIX": "fs_",
    "NEIGHBORS": 3, "MAX_DISTANCE": 10000, "DISCARD_NONMATCHING": False,
    "OUTPUT": "TEMPORARY_OUTPUT",
})["OUTPUT"]

by_school = {}
for f in three.getFeatures(QgsFeatureRequest().addOrderBy("distance")):
    by_school.setdefault(f["school_id"], []).append(f["distance"])

ties = {sid: d for sid, d in by_school.items()
        if len(d) > 1 and d[1] - d[0] < 150}
print(len(ties), "schools whose two nearest stations are within 150 m of each other")

Breakdown: With NEIGHBORS above one, each input feature appears once per neighbour found, ranked by n, so the output has more rows than the input — plan any later aggregation accordingly. Ordering by distance and grouping by the school's own id gives a sorted list of distances per school. A gap under 150 m between first and second is well within the uncertainty of where a building's centroid is, so for those schools "the nearest station" is a coin toss; flagging them is more honest than reporting whichever one the algorithm picked. For response-time questions, straight-line distance is only a first cut anyway — network service areas answer the real one.

The same join in Python with a spatial index

When the join is part of a larger script, or when you need behaviour the algorithm does not offer — skipping candidates with a certain attribute, distance to polygon edges rather than centroids, custom tie-breaking — use QgsSpatialIndex.nearestNeighbor directly.

from qgis.core import QgsSpatialIndex, QgsFeatureRequest

index = QgsSpatialIndex(stations.getFeatures(
    QgsFeatureRequest().setFilterExpression("\"crew_type\" = 'full-time'")),
    flags=QgsSpatialIndex.FlagStoreFeatureGeometries)

nearest = {}
for school in schools.getFeatures():
    geom = school.geometry()
    candidate_ids = index.nearestNeighbor(geom, 1, 10000)
    if not candidate_ids:
        nearest[school["school_id"]] = (None, None)
        continue
    best = min(candidate_ids, key=lambda fid: index.geometry(fid).distance(geom))
    nearest[school["school_id"]] = (best, round(index.geometry(best).distance(geom)))

unmatched = [k for k, v in nearest.items() if v[0] is None]
print(len(unmatched), "schools with no full-time station within 10 km")

Breakdown: Building the index from a filtered request restricts candidates to full-time stations without creating an intermediate layer. FlagStoreFeatureGeometries keeps geometries inside the index so exact distances can be computed without fetching features again. nearestNeighbor(geometry, k, maxDistance) may return more than k ids when several are equally near, which is why the code takes the minimum by exact distance rather than trusting the first id. Passing a geometry rather than a point means polygon inputs are measured edge to edge, not centroid to centroid — the right behaviour for "nearest hydrant to a building". The general pattern is covered in building and using a spatial index.

Write the result back and summarise coverage

The dictionary is only useful once it is on the map. Adding two fields and updating them in one edit session makes the result part of the school layer itself.

Distances become a coverage answerA stacked horizontal bar of 412 schools. 318 are within 3 km of a full-time station, 81 are between 3 and 10 km, and 13 have no full-time station within 10 km. The 13 are the policy finding, and they only appear because the distance cap was set and non-matching features were kept.412 schools by distance to a full-time station318 within 3 km81under 3 km77%3 – 10 km20%none within 10 km: 13the finding

from qgis.core import QgsField, edit
from qgis.PyQt.QtCore import QMetaType

with edit(schools):
    for name, kind in (("nearest_fs", QMetaType.Type.LongLong), ("fs_dist_m", QMetaType.Type.Int)):
        if schools.fields().indexOf(name) < 0:
            schools.addAttribute(QgsField(name, kind))
    schools.updateFields()
    i_fs = schools.fields().indexOf("nearest_fs")
    i_d = schools.fields().indexOf("fs_dist_m")
    for school in schools.getFeatures():
        fid, dist = nearest[school["school_id"]]
        schools.changeAttributeValues(school.id(), {i_fs: fid, i_d: dist})

bands = {"under 3 km": 0, "3–10 km": 0, "none within 10 km": 0}
for _, dist in nearest.values():
    key = ("none within 10 km" if dist is None
           else "under 3 km" if dist < 3000 else "3–10 km")
    bands[key] += 1
print(bands)

Breakdown: Adding fields and changing values inside one edit block commits them together or not at all, as described in editing features with transactions. changeAttributeValues takes a dictionary of field index to value, so both fields update in one call per feature. Writing None leaves the fields null for unmatched schools, which keeps them distinguishable from a real distance of zero. The banded counts are the one-line answer a report needs, and they only exist because the distance cap was set and unmatched features were kept.

QGIS version compatibility

native:joinbynearest has been available since QGIS 3.8, with NEIGHBORS and MAX_DISTANCE from the start. QgsSpatialIndex.nearestNeighbor accepts a geometry and a maximum distance from 3.8 as well; earlier releases take only a point and a count. QMetaType.Type field types are required on the QGIS 4 series; use QVariant.LongLong and QVariant.Int on releases before 3.38.

Troubleshooting

  • Every feature matches, including remote ones. No MAX_DISTANCE, or the layers are in degrees.
  • Distances are tiny decimals. Geographic CRS; reproject both layers to a projected CRS.
  • The output has more rows than the input. NEIGHBORS is above one; filter to n = 1 for one row per input.
  • Results differ between runs for a few features. Exact ties are resolved in index order; break ties explicitly by an attribute if it matters.
  • Polygons measured from centroids. Use the Python route with geometries in the index for edge-to-edge distance.

Conclusion

Put both layers in a projected CRS, always set a maximum distance, and keep unmatched features — they are usually the answer. Ask for two or three neighbours to reveal near ties, drop to QgsSpatialIndex.nearestNeighbor when you need filtering or custom tie-breaking, and write the nearest id and distance back to the layer so the result can be mapped and summarised.

Frequently Asked Questions

Is nearest-neighbour distance the same as travel distance? No. It is straight-line distance. Use network analysis when access follows roads.

Can I join lines or polygons, not just points? Yes. Distances are measured between geometries, so a polygon input measures to its nearest edge.

How do I find the nearest feature within the same layer? Use the same layer for both inputs with NEIGHBORS set to 2 and discard n = 1, which is each feature matching itself.

Why is the join slow on large layers? It is not usually — both routes use spatial indexes. Check that the second layer is not a remote service being fetched repeatedly; copy it locally first.

What distance cap should I choose? Choose it from the question, not the data. For emergency cover it is the distance that still meets a response standard; for assigning weather stations it is the distance over which conditions are plausibly similar. Record the cap in the output layer's metadata so readers know what "no match" means.