Spatial indexes g5

PostgreSQL with PostGIS
1.5 hours

Prerequisites: Spatial joins (g3), Indexes and query plans (d3)

See when a spatial index does the work and when it can't: nearest neighbours with <->, degrees that aren't metres, casts that hide a column from its index, and how to read a spatial plan. Everything runs in your own PostgreSQL with PostGIS, on the lab tables (the setup page).

Concept

When a GiST index does the work in a spatial query, and when it can't:

What the index holds (Indexes and query plans, d3). A GiST index on a geometry column keeps each shape's bounding box in a tree of boxes. It answers questions about boxes: && (two boxes overlap), @ and ~ (one box inside the other), and <-> (how far apart). ST_Intersects, ST_Covers and ST_DWithin add the box test for you, so the plan shows an Index Cond on boxes and a Filter (or a recheck) with the exact test on the rows the boxes let through.
Nearest neighbours. ORDER BY geom <-> point LIMIT k asks for the k shapes nearest a point. With the index, PostgreSQL walks the tree outward from the point, nearest boxes first, and stops after k rows: the plan is an Index Scan … Order By under a Limit, with no Sort and no search radius to guess (ArcGIS Pro's Near without one). The point must be a constant or come from an outer row (CROSS JOIN LATERAL), and <-> must compare the indexed column itself.
Degrees aren't metres. On SRID 4326, <-> measures in degrees, as if a degree of longitude were as long as a degree of latitude. It is only at the equator: at latitude φ a degree of longitude is cos φ as long on the ground, so <-> counts east–west distances 1 / cos φ too long (1.36 times at Hobart, 2 at 60°). The order it gives can differ from the order in metres. Either take more candidates by <-> and sort those by their distance in metres, or order by distance on geography, whose <-> measures on a sphere in metres.
A cast hides the column. An index serves the expression it was built on. ST_DWithin(geom::geography, …) compares geom::geography, not geom, so the index on geom can't help and every row is read (Distance and overlay, g4, met this at Costa Rica's borders). Build the index on the same expression, CREATE INDEX … USING gist ((geom::geography)), and the condition can use it. The other fix is a column of type geography with its own index.
When an index can't help: a test it can't turn into a box (ST_Distance(a, b) < d, d3); a cast or function around the column; a filter that keeps most of the table, which reads every page anyway; and statistics that mislead the planner (run ANALYZE after a load, d3).

Worked example

The five records nearest Hobart, Tasmania (147.33° E, 42.88° S), with their distance in kilometres on the ellipsoid:

Checked 2026-11-02 on PostgreSQL 16.15, PostGIS 3.4.2, GEOS 3.12.1, GDAL 3.8.4.

g5-nearest · PostgreSQL with PostGIS

WITH hobart AS (SELECT ST_SetSRID(ST_Point(147.33, -42.88), 4326) AS geom)
SELECT o.id, s.common_name,
       round((ST_Distance(o.geom::geography, h.geom::geography) / 1000)::numeric, 2) AS km
FROM occurrences AS o
JOIN species AS s ON s.id = o.species_id
CROSS JOIN hobart AS h
ORDER BY o.geom <-> h.geom, o.id
LIMIT 5;

-- result: g5-nearest · 5 rows · sha256 80d98a44e701
id | common_name | km
6004135246 | Tasmanian Devil | 2.44
6004145372 | Tasmanian Devil | 4.40
1647204917 | Tasmanian Devil | 5.57
1647151244 | Tasmanian Devil | 7.08
1647151259 | Tasmanian Devil | 7.08

In g5.sql, with the blocks it runs after.

The nearest records are Tasmanian devils, 2.44 to 7.08 km away. The last two lie at the same spot, with nine more records there, so o.id breaks the tie and the same five come back every time. <-> orders; ST_Distance on geography gives the distance in metres, worked out for the five rows only.

Lab

Each block runs in a transaction that is rolled back, as in d3: in psql, type BEGIN; before a block and ROLLBACK; after it. Each plan lab lists what its plan must show; the plan printed under it is an example, from the versions named above the first block.

Statistics first, as in d3: ANALYZE reading every row, so the planner's estimates, and the plans below, are the same on your server:

g5-analyze · PostgreSQL with PostGIS

SET LOCAL default_statistics_target = 1000;
ANALYZE occurrences;
SELECT c.reltuples::bigint AS planner_rows, (SELECT count(*) FROM occurrences) AS table_rows
FROM pg_class AS c
WHERE c.oid = 'occurrences'::regclass;

-- result: g5-analyze · 1 row · sha256 2415e801456a
planner_rows | table_rows
40395 | 40395

In g5.sql, with the blocks it runs after.

The worked example's plan, without the names and distances. The plan must show the top node returning 5 rows, an index scan on occurrences_geom_idx, and its Order By:
g5-knn-plan · PostgreSQL with PostGIS

Runs after: g5-analyze
```
EXPLAIN (ANALYZE, COSTS OFF, TIMING OFF, SUMMARY OFF)
SELECT o.id
FROM occurrences AS o
ORDER BY o.geom <-> ST_SetSRID(ST_Point(147.33, -42.88), 4326)
LIMIT 5;
```
What the plan must show:
- the top node returns 5 rows (actual rows=5)
- a line with Index Scan using occurrences_geom_idx
- a line with Order By
An example plan, from the versions above; yours passes if it shows the lines listed, whatever else differs.
```
-- result: g5-knn-plan · 3 rows · sha256 c1ec5c12ff81
QUERY PLAN
Limit (actual rows=5 loops=1)
  ->  Index Scan using occurrences_geom_idx on occurrences o…
        Order By: (geom <-> '0101000020E6100000C3F5285C8F6A6…
```
In g5.sql, with the blocks it runs after.
The index scan reads the records in order of distance and stops when Limit has five: 5 rows read of 40,395, and no Sort.

Does the order in degrees match the order in metres? Rank every record both ways, and keep the rows where the ranks differ in either top 20:

g5-degrees · PostgreSQL with PostGIS

WITH hobart AS (SELECT ST_SetSRID(ST_Point(147.33, -42.88), 4326) AS geom),
ranked AS (
  SELECT o.id,
         ST_Distance(o.geom::geography, h.geom::geography) / 1000 AS km,
         o.geom <-> h.geom AS degrees,
         ST_X(o.geom) - ST_X(h.geom) AS east_deg,
         ST_Y(o.geom) - ST_Y(h.geom) AS north_deg,
         row_number() OVER (ORDER BY o.geom <-> h.geom, o.id) AS by_degrees,
         row_number() OVER (ORDER BY ST_Distance(o.geom::geography, h.geom::geography), o.id) AS by_metres
  FROM occurrences AS o
  CROSS JOIN hobart AS h
)
SELECT id, round(km::numeric, 3) AS km, round(degrees::numeric, 4) AS degrees,
       round(east_deg::numeric, 3) AS east_deg, round(north_deg::numeric, 3) AS north_deg,
       by_degrees, by_metres
FROM ranked
WHERE (by_degrees <= 20 OR by_metres <= 20) AND by_degrees <> by_metres
ORDER BY by_metres;

-- result: g5-degrees · 6 rows · sha256 f2b1681a3b45
id | km | degrees | east_deg | north_deg | by_degrees | by_metres
6004135264 | 7.997 | 0.0844 | -0.065 | -0.054 | 30 | 16
6004133344 | 8.003 | 0.0738 | 0.023 | -0.070 | 16 | 17
6004133858 | 8.003 | 0.0738 | 0.023 | -0.070 | 17 | 18
6004134874 | 8.003 | 0.0738 | 0.023 | -0.070 | 18 | 19
6004136872 | 8.003 | 0.0738 | 0.023 | -0.070 | 19 | 20
… 1 more row

In g5.sql, with the blocks it runs after.

Record 6004135264 is 7.997 km from the point, the 16th nearest in metres, but 30th by <->. It lies to the south-west, 0.065° west and 0.054° south: 0.0844° by <->. Thirteen records 8.003 km away, nearly due south (0.023° east, 0.070° south), are only 0.0738° away by <->, so they take ranks 16 to 28 ahead of it, and so does one 8.185 km away. At 42.88° south a degree of longitude is about 81.7 km on the ground and a degree of latitude about 111.1 km, so <-> counts the west part 1.36 times too long. The 20 nearest by degrees miss the 16th nearest in metres and take the 21st instead. The difference here is metres, but it grows towards the poles and with the distances involved; a list ranked by <-> on SRID 4326 is a list of candidates, not an answer in metres.

The records within 20 km of Hobart, measured on the ellipsoid. The cast hides geom from its index, so every row is read. The plan must show 55 rows, a Seq Scan, and no index:
g5-cast-plan · PostgreSQL with PostGIS

Runs after: g5-analyze
```
EXPLAIN (ANALYZE, COSTS OFF, TIMING OFF, SUMMARY OFF)
SELECT o.id
FROM occurrences AS o
WHERE ST_DWithin(o.geom::geography, ST_SetSRID(ST_Point(147.33, -42.88), 4326)::geography, 20000);
```
What the plan must show:
- the top node returns 55 rows (actual rows=55)
- a line with Seq Scan on occurrences
- no line with Index
An example plan, from the versions above; yours passes if it shows the lines listed, whatever else differs.
```
-- result: g5-cast-plan · 3 rows · sha256 05bf25e03416
QUERY PLAN
Seq Scan on occurrences o (actual rows=55 loops=1)
  Filter: st_dwithin((geom)::geography, '0101000020E6100000C…
  Rows Removed by Filter: 40340
```
In g5.sql, with the blocks it runs after.
PostgreSQL measured the distance from all 40,395 records to keep 55.
An index on the expression itself fixes it. The plan must show the same 55 rows and a Bitmap Index Scan on the new index:
g5-geography-index · PostgreSQL with PostGIS

Runs after: g5-analyze
```
CREATE INDEX occurrences_geography_idx ON occurrences USING gist ((geom::geography));
EXPLAIN (ANALYZE, COSTS OFF, TIMING OFF, SUMMARY OFF)
SELECT o.id
FROM occurrences AS o
WHERE ST_DWithin(o.geom::geography, ST_SetSRID(ST_Point(147.33, -42.88), 4326)::geography, 20000);
```
What the plan must show:
- the top node returns 55 rows (actual rows=55)
- a line with Bitmap Index Scan on occurrences_geography_idx
An example plan, from the versions above; yours passes if it shows the lines listed, whatever else differs.
```
-- result: g5-geography-index · 6 rows · sha256 00b2f96fdafb
QUERY PLAN
Bitmap Heap Scan on occurrences o (actual rows=55 loops=1)
  Filter: st_dwithin((geom)::geography, '0101000020E6100000C…
  Rows Removed by Filter: 11
  Heap Blocks: exact=19
  ->  Bitmap Index Scan on occurrences_geography_idx (actual…
… 1 more row
```
In g5.sql, with the blocks it runs after.
The index finds 66 records whose boxes come within 20 km of the point, the exact test keeps the same 55, and only the pages holding the 66 are read (Heap Blocks), not the whole table. The query didn't change: only the index did.

The 999 records in no outline in Spatial joins, g3: how far is each from the nearest outline? Nearest-neighbour search works on polygons too, once per record through CROSS JOIN LATERAL:

g5-coast · PostgreSQL with PostGIS

WITH inside AS (
  SELECT DISTINCT o.id FROM country_shapes AS c JOIN occurrences AS o ON ST_Covers(c.geom, o.geom)
), nearest AS (
  SELECT o.id, o.country_code_gbif, n.country_code,
         ST_Distance(o.geom::geography, n.geom::geography) / 1000 AS km
  FROM occurrences AS o
  CROSS JOIN LATERAL (SELECT c.country_code, c.geom FROM country_shapes AS c ORDER BY c.geom <-> o.geom LIMIT 1) AS n
  WHERE o.id NOT IN (SELECT id FROM inside)
)
SELECT CASE WHEN km < 1 THEN 'under 1 km' WHEN km < 5 THEN '1 to 5 km' WHEN km < 25 THEN '5 to 25 km' ELSE '25 km or more' END AS distance,
       count(*) AS records,
       count(*) FILTER (WHERE country_code = country_code_gbif) AS nearest_is_gbif_country,
       round(max(km)::numeric, 1) AS max_km
FROM nearest
GROUP BY 1
ORDER BY min(km);

-- result: g5-coast · 4 rows · sha256 ef446e9ae320
distance | records | nearest_is_gbif_country | max_km
under 1 km | 236 | 226 | 1.0
1 to 5 km | 200 | 193 | 4.5
5 to 25 km | 529 | 529 | 24.8
25 km or more | 34 | 31 | 343.9

In g5.sql, with the blocks it runs after.

For a polygon, <-> is the distance to its nearest edge, not to its box, so LIMIT 1 gets the nearest outline, by degrees. Step 3's caution applies: here the nearest outline by degrees is also the nearest in metres for all 999 records (checked against each record's next two). 436 records lie within 5 km of an outline: beaches, mangroves and coasts the 1:50m outlines simplify. In nearly all of them the nearest outline is GBIF's country. 529 lie 5 to 25 km out, mostly on islands too small for 1:50m: 346 bilby records on one island off South Australia, 133 pig-nosed turtle records on the Gulf of Papua's delta. 34 are 25 km or more from any outline, up to 343.9 km, and are worth reading one by one.

Result

The nearest records to Hobart are Tasmanian devils, 2.44 to 7.08 km away, found by an Index Scan … Order By that reads 5 rows. <-> on SRID 4326 ranks in degrees: at 43° south it put the 16th nearest record in metres 30th. With the cast to geography, a 20 km search reads all 40,395 rows to keep 55; an index on geom::geography finds 66 candidates and keeps the same 55. Of g3's 999 records in no outline, 436 lie within 5 km of one, 529 more within 25 km, mostly on small islands, and 34 farther.

Exercise

Find the 20 records nearest Hobart in metres, with an index doing the work: no Sort and no Seq Scan. Show the plan; then run the query and find g5-degrees' record 6004135264 in it.

Show answer

g5-ex · PostgreSQL with PostGIS

Runs after: g5-analyze

CREATE INDEX occurrences_geography_idx ON occurrences USING gist ((geom::geography));
EXPLAIN (ANALYZE, COSTS OFF, TIMING OFF, SUMMARY OFF)
SELECT o.id
FROM occurrences AS o
ORDER BY o.geom::geography <-> ST_SetSRID(ST_Point(147.33, -42.88), 4326)::geography
LIMIT 20;

What the plan must show:

the top node returns 20 rows (actual rows=20)
a line with Index Scan using occurrences_geography_idx
a line with Order By: ((geom)::geography <->
no line with Sort
no line with Seq Scan

An example plan, from the versions above; yours passes if it shows the lines listed, whatever else differs.

-- result: g5-ex · 3 rows · sha256 ce122751eb76
QUERY PLAN
Limit (actual rows=20 loops=1)
  ->  Index Scan using occurrences_geography_idx on occurren…
        Order By: ((geom)::geography <-> '0101000020E6100000…

In g5.sql, with the blocks it runs after.

Check

The answer: an index on geom::geography, and ORDER BY the same expression <-> the point as geography. <-> on geography measures on a sphere, in metres, and the index walks outward from the point as in step 2.

Compare your plan with the answer's: a Limit over an Index Scan on your index with actual rows=20, its Order By on (geom)::geography <->, and no Sort or Seq Scan anywhere. Run the query itself: record 6004135264 is 16th, as in g5-degrees' by_metres. On a sphere the distances differ from ST_Distance's on the ellipsoid by a few metres in 7 km (7,088 m against 7,081 m for the eleven records at one spot), which can swap two records whose distances are that close; sort the 20 by ST_Distance if that matters. Ordering by ST_Distance(geom::geography, …) instead gives the right rows too, but measures all 40,395 and sorts them. Taking 40 candidates by <-> on geom and sorting them by metres also finds the 16th here, through the existing index, but nothing in the query says 40 is enough.

Dialects

PostgreSQL: GiST on geometry and geography (and SP-GiST and BRIN on geometry); <-> orders nearest first through GiST, in the column's units on geometry and in metres on a sphere on geography; <#> orders by the distance between boxes. An expression index serves only the same expression; ST_DWithin, ST_Intersects and the other predicates use either argument's index by themselves.
DuckDB: the spatial extension (not in this course's browser) builds an R-tree with CREATE INDEX … USING RTREE (geom), used to filter a table by a constant shape (ST_Intersects, ST_Within and the like). For the nearest rows, ORDER BY ST_Distance(geom, point) LIMIT k measures every row, in the coordinates' units; ST_Distance_Spheroid measures in metres (latitude first, g1).
SQL Server: spatial indexes are grids, not trees; one on geometry needs a BOUNDING_BOX. A nearest-neighbour query uses one when it is SELECT TOP (k) … WHERE geom.STDistance(@p) IS NOT NULL ORDER BY geom.STDistance(@p); in metres on geography.