Nearest neighbours in PostGIS

“Find the five places closest to this point” sounds simple, and the obvious way to write it is slow. Sorting every row by ST_Distance makes PostgreSQL compute the distance to every row in the table, then sort them all, to return five. PostGIS has an operator that lets the index do the ordering instead.

The <-> operator and a GiST index

<-> returns the distance between two geometries. Its point is that a GiST index can return rows already in order of that distance, nearest first, so a query with LIMIT 5 reads only a handful of rows from the index and stops:

CREATE INDEX places_geom_idx ON places USING gist (geom);

SELECT id, name
FROM places
ORDER BY geom <->
  ST_SetSRID(ST_MakePoint(-0.1276, 51.5072), 4326)
LIMIT 5;

Two conditions decide whether the index is used:

• The operator has to be in ORDER BY. In a WHERE clause or a SELECT list it’s an ordinary distance function.
• One side has to be fixed for the whole scan: a literal, a parameter, or a scalar subquery that returns one value. A column from another table in a plain join doesn’t count.

EXPLAIN shows which you got. The fast plan is an Index Scan using places_geom_idx with an Order By: (geom <-> ...) line under it. A Sort node over a sequential scan means the index isn’t doing the ordering.

That rule seems to rule out “the nearest hospital to each school”, but a LATERAL join gets around it. For each outer row, the inner query runs with that row’s geometry as a fixed value, so every inner scan is an index-ordered LIMIT 1:

SELECT s.name, h.name AS nearest_hospital
FROM schools s
CROSS JOIN LATERAL (
  SELECT name
  FROM hospitals
  ORDER BY hospitals.geom <-> s.geom
  LIMIT 1
) h;

One more trap sits in that first query: ST_MakePoint takes longitude first, then latitude. Coordinates pasted from a map are usually the other way round, and the query still runs, just somewhere else on the planet.

Geometry or geography

What <-> measures depends on the column’s type.

On geometry, it’s the planar distance in the units of the coordinate system. For longitude and latitude in SRID 4326, that unit is degrees, and a degree of longitude shrinks from about 111 km at the equator to nothing at the poles. Ordering by degrees is close to ordering by metres near the equator and increasingly wrong away from it, so the “nearest” row can be the second or third nearest on the ground. That includes the London query above.

There are two ways out:

• A geography column. Distances are in metres. <-> on geography uses a sphere, which is fast and good for ordering. ST_Distance on geography uses the spheroid by default, which is more exact. Order with <->, and show ST_Distance when the number itself matters; the two can disagree slightly on near ties.
• A projected geometry column in a metric coordinate system that fits the area, such as the local UTM zone. Planar maths is then both fast and accurate, as long as the data stays inside that area.

For data spread across a country or the world, geography is the simpler choice. For a city or a region, a projected geometry is the usual one.

ST_DWithin as a prefilter

“Nearest” with no limit on distance will happily return a place 800 km away when nothing is closer. Usually the real question is “the nearest few, within 5 km”. The tempting way to write the radius uses ST_Distance in the WHERE clause, which computes every distance again. ST_DWithin asks the same question in a way the index can answer, because it starts with a bounding-box test:

-- $1 is the search point, as geography; places has a GiST index on geog
SELECT id, name, ST_Distance(geog, $1) AS metres
FROM places
WHERE ST_DWithin(geog, $1, 5000)
ORDER BY geog <-> $1
LIMIT 10;

On geography the distance is in metres, and the spheroid is used by default. On a geometry column in SRID 4326 the same 5000 means 5,000 degrees, which covers the whole planet and filters nothing. It returns results, just not the ones asked for.

Sources: the PostGIS reference for <-> and ST_DWithin, and the PostGIS workshop’s nearest-neighbour chapter.