Using KNN operators and lateral joins in PostGIS
This page builds a nearest-neighbour model the way PostGIS intends: a LATERAL subquery per row, ordered by the index-backed <-> operator, limited to k, with exact distances computed afterwards on the handful of survivors.
When to use this approach
- Each row needs its own k nearest matches. “The three closest depots to every ping” is not a filter, it is a per-row ranking, and a
WHEREclause cannot express it. - A proximity join is timing out. The usual cause is
ORDER BY ST_Distanceover a cross join, which computes every pair before sorting. - A radius is the wrong shape for the question.
ST_DWithinanswers “within 500 m”; if the answer must exist even when nothing is within 500 m, you want k-nearest instead. The broader topic is optimizing proximity joins.
Prerequisites
- PostGIS 2.2+ for the index-backed
<->operator on geometry (recentre versions also support it on geography). - A GiST index on the candidate side’s geometry column — the operator is only fast because of it.
- Both sides in the same SRID, projected if distances must be metres.
ANALYZErun after loading, so the planner costs the lateral correctly.
Step-by-step instructions
1. Write the lateral, ordered by the operator
-- models/intermediate/int_ping_nearest_depots.sql
{{ config(
materialized = 'table',
post_hook = "ANALYZE {{ this }}"
) }}
select
p.ping_id,
n.depot_id,
n.rank,
round(st_distance(p.geom::geography, n.geom::geography)::numeric, 1) as distance_m
from {{ ref('stg_trip_pings') }} as p
cross join lateral (
select
d.depot_id,
d.geom,
row_number() over () as rank
from {{ ref('stg_depots') }} as d
order by d.geom <-> p.geom -- index-backed nearest-neighbour ordering
limit 3
) as n
Three details make this work. CROSS JOIN LATERAL lets the inner query reference p, which is what allows a per-row LIMIT. The <-> operator in ORDER BY is what PostGIS recognises as a KNN search and answers from the GiST index. And the exact distance is computed outside the ordering, on three rows per ping rather than on every candidate.
Verify the plan shows an index scan driven by the operator:
explain (analyze, buffers)
select ... ; -- the compiled model SQL
Nested Loop (actual time=0.041..812.3 rows=3000000 loops=1)
-> Seq Scan on stg_trip_pings p (rows=1000000)
-> Limit (actual rows=3 loops=1000000)
-> Index Scan using stg_depots_geom_idx on stg_depots d
Order By: (geom <-> p.geom)
Order By: on the index scan line is the confirmation. If it says Sort followed by a sequential scan, the operator is not being used.
2. Understand what the operator actually orders by
<-> on geometry returns the distance between bounding-box centroids at index level and the true distance between geometries at the recheck. For points these are the same thing, and the ordering is exact. For polygons and lines they are not, so the index returns candidates in approximate order and the exact ordering must be re-established.
-- Correct for polygon candidates: over-fetch, then rank exactly
select
p.ping_id, n.zone_id, n.distance_m
from {{ ref('stg_trip_pings') }} as p
cross join lateral (
select
z.zone_id,
st_distance(p.geom, z.geom) as distance_m
from {{ ref('stg_zones') }} as z
order by z.geom <-> p.geom
limit 10 -- over-fetch candidates
) as c
cross join lateral (
select c.zone_id, c.distance_m
order by c.distance_m
limit 3 -- exact top 3
) as n
Over-fetching by a factor of three and re-ranking exactly is the standard remedy, and it is still far cheaper than a cross join.
Verify the approximation matters or does not, for your data, by comparing the two orderings on a sample:
-- Rows where the index order and the exact order disagree in the top 3
select count(*) from ... where index_rank <> exact_rank;
3. Keep distances in a unit someone can act on
round(st_distance(p.geom::geography, n.geom::geography)::numeric, 1) as distance_m
Casting to geography for the final distance gives geodesic metres regardless of the storage SRID, which is almost always the number a consumer wants. Do not put the cast in the ORDER BY — it would defeat the geometry index that makes the search fast. Order with the operator on geometry, measure with geography afterwards.
Verify the distances are plausible and in metres:
select min(distance_m), percentile_cont(0.5) within group (order by distance_m), max(distance_m)
from {{ ref('int_ping_nearest_depots') }};
4. Guard the model’s grain and cost
models:
- name: int_ping_nearest_depots
description: Exactly k rows per ping, ranked nearest first, distance in geodesic metres.
tests:
- dbt_utils.expression_is_true:
expression: "count(*) = (select count(*) * 3 from {{ ref('stg_trip_pings') }})"
columns:
- name: rank
tests:
- accepted_values:
values: [1, 2, 3]
- name: distance_m
tests:
- not_null
- dbt_utils.accepted_range:
min_value: 0
max_value: 200000
The row-count assertion is the useful one: a lateral with LIMIT 3 produces exactly three rows per input unless the candidate table is smaller than k, which is a real and easily-missed edge case.
5. Fall back to a radius when k is the wrong question
-- "the nearest depot, but only if it is within 5 km"
select p.ping_id, n.depot_id, n.distance_m
from {{ ref('stg_trip_pings') }} as p
left join lateral (
select d.depot_id, st_distance(p.geom::geography, d.geom::geography) as distance_m
from {{ ref('stg_depots') }} as d
where st_dwithin(p.geom::geography, d.geom::geography, 5000)
order by d.geom <-> p.geom
limit 1
) as n on true
LEFT JOIN LATERAL … ON true keeps the ping when nothing qualifies, with nulls in the depot columns — the difference between “no depot nearby” and “this ping does not exist”, which a CROSS JOIN LATERAL erases.
Verify how many inputs fall through:
select count(*) filter (where depot_id is null) as pings_with_no_depot_within_5km,
count(*) as total
from {{ ref('int_ping_nearest_depot_bounded') }};
Configuration reference
| Element | Where | Note |
|---|---|---|
<-> |
ORDER BY inside the lateral |
The only form PostGIS answers from the index |
CROSS JOIN LATERAL |
join clause | Drops inputs with no candidates |
LEFT JOIN LATERAL … ON true |
join clause | Keeps them, with nulls |
LIMIT k |
inside the lateral | Per-row limit; the whole reason for the lateral |
| over-fetch factor | inside the lateral | 3× k for polygon candidates, then re-rank exactly |
::geography |
final distance only | Geodesic metres; never in the ORDER BY |
Gotchas & edge cases
ORDER BY ST_Distance(...)is not KNN. It is a sort over the full candidate set. Only the<->operator triggers the index-backed search.- A
WHEREclause inside the lateral can defeat the index scan if it is selective enough that the planner prefers a filter-then-sort. Check the plan after adding one. <->on geography behaves differently by version. Older PostGIS returns sphere-based distances for the index ordering; verify with a known pair before relying on exact ordering.- The lateral runs once per input row. With a million inputs the constant factors matter, so keep the inner query narrow — select the id and geometry, nothing else.
- Ties are arbitrary without a tie-break. Add
, d.depot_idto the innerORDER BYso repeated runs agree, exactly as with the grain policies in spatial joins and predicate tuning.
FAQ
Why is <-> faster than sorting by distance?
Because it is not a sort. The GiST index is a hierarchy of bounding boxes, and the KNN search walks it outward from the query point, expanding the nearest box first and stopping once k results are certain. A sort must produce every distance before it can order them, so it cannot stop early.
Can this run on DuckDB or a cloud warehouse?
The lateral can; the index-backed operator cannot, because those engines have no persistent spatial index. There, nearest-neighbour work is usually restructured as a grid-key candidate search followed by an exact ranking — the approach in discrete global grid macros.
Should k-nearest results be materialized or computed on demand?
Materialized, if more than one consumer needs them or if the input table is large — the lateral is fast per row but still runs once per row. Materialize with a rank column so consumers can take the top one or the top three without re-running the search.
How do I handle candidates that should be excluded per input row?
Put the condition inside the lateral, referencing the outer row — that is precisely what a lateral allows. Watch the plan afterwards, because a correlated filter can change the planner’s mind about using the KNN index scan, and switch to over-fetch-then-filter if it does.
Related
- Optimizing Proximity Joins — the topic this implements.
- Speeding Up Nearest-Neighbor Joins in PostGIS — complementary tuning for the same workload.
- Spatial Joins & Predicate Tuning — when a radius predicate is the better tool.
Up: Part of Optimizing Proximity Joins.