Calculate the Distance Between Features in PyQGIS

"How far apart are these?" has three different correct answers depending on what you mean and where the data is. QgsGeometry.distance() measures in map units on a flat plane. QgsDistanceArea measures on the ellipsoid in real metres. And the distance between two features is the shortest distance between any pair of their points, which for polygons is zero as soon as they touch — not the distance between their centres.

This page is a focused recipe within Geometry Operations and Spatial Predicates in PyQGIS. It covers the planar and ellipsoidal measurements, what "distance between features" actually returns for lines and polygons, and how to find nearest neighbours without a quadratic scan.

Which distance does the question mean?The same two polygons are shown three times. In the first, a line joins their centroids, giving the longest of the three values. In the second, a shorter line joins their nearest edges, which is what QgsGeometry.distance returns. In the third the polygons share a boundary and the distance is zero even though their centroids are far apart.distance() is edge to edge — not centre to centrecentroid to centroidyou must compute this yourselfnearest edgeswhat distance() returnstouchingdistance() returns 0.0

Prerequisites

  • QGIS 3.34 LTR (bundled Python 3.12) or newer.
  • Both layers in the same CRS. distance() compares raw coordinates and does not transform.
  • A projected CRS if you want the planar answer to mean metres — see Coordinate Reference Systems in PyQGIS.

The planar measurement

distance() is the shortest distance between any point of one geometry and any point of the other, expressed in map units.

from qgis.core import QgsProject

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

school = next(schools.getFeatures()).geometry()
park = next(parks.getFeatures()).geometry()

print(round(school.distance(park), 1), "map units")
print(round(school.centroid().distance(park.centroid()), 1), "centroid to centroid")

Breakdown: For two points this is ordinary Euclidean distance. For a point and a polygon it is the distance to the nearest boundary point, and zero if the point is inside. For two polygons that touch or overlap it is exactly zero. If what you meant was centre to centre, take the centroids first — the second line shows how, and the two values are usually very different.

To know where the closest approach is, shortestLine() returns the connecting segment, which is useful both for measurement and for drawing:

connector = school.shortestLine(park)
print(connector.asWkt(precision=1))
print(round(connector.length(), 1), "same value as distance()")

Breakdown: shortestLine() returns a two-point QgsGeometry whose length equals distance(). Adding it to a memory layer makes the measurement visible on the map, which is how the built-in measure tool presents its result.

The ellipsoidal measurement

Planar distance is only meaningful in a projected CRS, and even there it ignores the curvature of the Earth over long spans. QgsDistanceArea does the real calculation.

from qgis.core import QgsDistanceArea, QgsProject, QgsUnitTypes

calculator = QgsDistanceArea()
calculator.setSourceCrs(schools.crs(), QgsProject.instance().transformContext())
calculator.setEllipsoid(QgsProject.instance().ellipsoid())

metres = calculator.measureLine(school.asPoint(), park.centroid().asPoint())
km = calculator.convertLengthMeasurement(metres, QgsUnitTypes.DistanceUnit.DistanceKilometers)
print(f"{metres:,.1f} m  ({km:.3f} km)")

Breakdown: setSourceCrs() tells the calculator what the incoming coordinates mean, and setEllipsoid() is what switches it from planar to ellipsoidal arithmetic — without that second call it silently returns the same planar number as distance(). measureLine() takes QgsPointXY objects and always returns metres regardless of the source CRS, which is exactly what makes it safe to use on lat/lon data. convertLengthMeasurement() handles unit conversion so no magic constants appear in your code.

The difference matters at scale. Over a few hundred metres in a well-chosen UTM zone, planar and ellipsoidal agree to within centimetres. Over hundreds of kilometres, or near the edge of a projection's valid band, the divergence becomes metres and then tens of metres.

Find the nearest feature

The naive nearest-neighbour search compares every pair. A QgsSpatialIndex answers it directly, because nearest-neighbour is one of the queries an R-tree is built for.

Brute-force scan versus nearestNeighbor()On the left, a query point is joined by lines to all fourteen candidate features, each requiring a distance calculation. On the right, the spatial index returns three nearest candidate identifiers immediately, and only those three are measured exactly to pick the winner.The index narrows the field before any distance is computedevery pair measured14 distance calculationsindex.nearestNeighbor(pt, 3)3 distance calculations

from qgis.core import QgsSpatialIndex, QgsProject

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

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

for school in schools.getFeatures():
    geometry = school.geometry()
    candidates = index.nearestNeighbor(geometry, 3)
    best_id, best_distance = None, float("inf")
    for park_id in candidates:
        d = geometry.distance(index.geometry(park_id))
        if d < best_distance:
            best_id, best_distance = park_id, d
    print(school["name"], "->", best_id, round(best_distance, 1))

Breakdown: FlagStoreFeatureGeometries makes the index keep the geometries so index.geometry(id) can return them without a second provider read — a large saving inside a loop, at the cost of memory. nearestNeighbor(geometry, 3) returns candidate IDs ranked by bounding-box proximity, which is an approximation, so the exact distance() still decides the winner. Asking for three rather than one is the guard against that approximation: the box-nearest feature is not always the truly nearest, particularly for long thin shapes.

For a whole-layer answer, the Processing algorithm does all of this and writes a joined output:

import processing

joined = processing.run("native:joinbynearest", {
    "INPUT": "/data/schools.gpkg|layername=schools",
    "INPUT_2": "/data/parks.gpkg|layername=parks",
    "FIELDS_TO_COPY": ["park_name"],
    "NEIGHBORS": 1,
    "MAX_DISTANCE": 2000,
    "OUTPUT": "/data/output/schools_nearest_park.gpkg",
})["OUTPUT"]

Breakdown: native:joinbynearest adds a distance field to the output automatically, along with the copied attributes. MAX_DISTANCE caps the search, which both speeds it up and leaves genuinely isolated features unjoined rather than matched to something absurdly far away. NEIGHBORS above 1 produces one output feature per neighbour, which is how you build a "three nearest" table.

Measure the length of a line

Line length is the same question in one dimension, and it shows the planar-versus-ellipsoidal divergence more clearly than point distance because the error accumulates along every segment.

Where planar length stops being good enoughTwo rows. A two kilometre local line shows planar and ellipsoidal lengths agreeing to within a few centimetres. A twelve hundred kilometre route shows the planar figure reading several kilometres short of the ellipsoidal one, with the gap marked as accumulated error.The error is proportional to the span, not the vertex counta 2 km survey lineplanar (distance / length)2 041.36 mellipsoidal (QgsDistanceArea)2 041.41 m5 cm apart — ignore ita 1 200 km routeplanar (distance / length)1 194.2 kmellipsoidal (QgsDistanceArea)1 201.7 km7.5 km short

from qgis.core import QgsDistanceArea, QgsProject, QgsUnitTypes

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

calculator = QgsDistanceArea()
calculator.setSourceCrs(routes.crs(), QgsProject.instance().transformContext())
calculator.setEllipsoid(QgsProject.instance().ellipsoid())

for feature in routes.getFeatures():
    geometry = feature.geometry()
    planar = geometry.length()
    ellipsoidal = calculator.measureLength(geometry)
    km = calculator.convertLengthMeasurement(
        ellipsoidal, QgsUnitTypes.DistanceUnit.DistanceKilometers
    )
    drift = abs(ellipsoidal - planar) / ellipsoidal if ellipsoidal else 0
    print(f"{feature['route_id']}: {km:.2f} km  (planar differs by {drift:.2%})")

Breakdown: geometry.length() sums the straight-line distances between consecutive vertices in map units — fast, and correct only in a projected CRS over a modest span. measureLength() does the same walk on the ellipsoid and always returns metres. Reporting the relative drift rather than the absolute difference makes the number comparable across routes of different lengths, and gives you a threshold to act on: under about 0.1 % the planar figure is fine for most purposes, above 1 % it is misleading. Note this measures the length of the geometry as digitised — it is not a travel distance along a network, which needs a routing algorithm rather than a measurement.

QGIS version compatibility

The examples target QGIS 3.34 LTR (Python 3.12).

QGIS versionPythonNotes
3.28 LTR3.9Identical API. QgsSpatialIndex.FlagStoreFeatureGeometries available since 3.4.
3.34 LTR3.12Baseline for this page.
3.40 / 3.443.12QgsUnitTypes.DistanceUnit moved under Qgis; the old spelling still resolves.

QgsDistanceArea and the distance() / shortestLine() methods are unchanged across 3.x.

Troubleshooting

  • Distances look like small decimals. The layer is in a geographic CRS, so the answer is in degrees. Reproject, or switch to QgsDistanceArea with an ellipsoid set.
  • QgsDistanceArea returns the same number as distance(). setEllipsoid() was never called, so it stayed in planar mode. Pass QgsProject.instance().ellipsoid().
  • Two adjacent polygons report a distance of zero. That is correct — they touch. Use centroids if you wanted centre-to-centre separation.
  • The nearest neighbour is obviously wrong. Only one candidate was requested from the index, and the bounding-box ranking picked a long thin feature. Ask for three to five and measure them exactly.
  • Everything is far apart. The two layers are in different CRSs. Compare crs().authid() on both before measuring.
  • A nearest join left rows unmatched. MAX_DISTANCE excluded them, which is usually intentional. Raise it or remove it if every feature must match.

Conclusion

Pick the measurement that matches the question: distance() for planar edge-to-edge in map units, QgsDistanceArea with an ellipsoid for real-world metres, and centroids only when you explicitly want centre-to-centre. For nearest-neighbour work, let a spatial index shortlist the candidates and confirm the winner with an exact measurement.

Frequently Asked Questions

Why is the distance between two adjacent parcels zero? Because distance() measures the shortest gap between the geometries, and two polygons sharing a boundary have no gap. If you want the separation of their centres, measure a.centroid().distance(b.centroid()) instead.

When do I need QgsDistanceArea rather than distance()? Whenever the data is in a geographic CRS, or when the span is large enough that the curvature of the Earth matters. distance() is planar and returns map units; QgsDistanceArea with an ellipsoid returns real metres.

Does nearestNeighbor give the truly nearest feature? Not guaranteed. It ranks candidates by bounding-box proximity, which can mis-order long or irregular shapes. Request several candidates and compare them with an exact distance() call.

Does distance() work between different geometry types? Yes. Point to polygon, line to point and polygon to polygon all work, and each returns the shortest gap between any point of one and any point of the other — zero if they touch or overlap.

How do I see where the closest approach is?a.shortestLine(b) returns the connecting two-point geometry. Its length equals distance(), and adding it to a memory layer draws the measurement on the map.

Can I measure distance along a line rather than straight through? Yes — lineLocatePoint() gives the distance along a line to the nearest point on it, which is what linear referencing needs. Straight-line distance ignores the route entirely.