Spatial joins g3

PostgreSQL with PostGIS
2 hours

Prerequisites: Loading and validating shapes (g2)

Join tables by where things are, and check one source against another: ST_Contains, ST_Covers and ST_Intersects, which side a spatial join starts from, the records on an edge, and the records two overlapping outlines count twice. Everything runs in your own PostgreSQL with PostGIS, on the lab tables (the setup page).

Concept

How to join tables by where things are, and how to check one source against another:

A spatial predicate answers true or false about two shapes. Three cover most work:
- ST_Contains(a, b): b lies inside a, and some of b is in a's interior. A point on a's edge doesn't count.
- ST_Covers(a, b): no point of b lies outside a. A point on the edge counts.
- ST_Intersects(a, b): the two share at least one point, edges included. For a point and a polygon it answers the same as ST_Covers.
- ST_Within(b, a) is ST_Contains(a, b) the other way round, and ST_CoveredBy(b, a) is ST_Covers(a, b). Each predicate is a pattern on the DE-9IM matrix (Loading and validating shapes, g2): ST_Contains is T*****FF*, the interiors meet (T: a point, a line or an area) and nothing of b meets a's exterior.
A spatial join is a join whose ON is a predicate, ArcGIS Pro's Spatial Join tool written as a query: country_shapes AS c JOIN occurrences AS o ON ST_Covers(c.geom, o.geom) pairs each record with every outline that covers it. A record in no outline drops out of an inner join, and a record in two outlines comes back twice.
Which table the join starts from changes its speed. While consecutive tests use the same polygon, PostGIS keeps a prepared copy of it, with an index of its edges, so each further point costs little. Started from the outlines, each outline is prepared once and tested against the records its spatial index (Indexes, d3) finds inside its bounding box. Started from the records (occurrences LEFT JOIN country_shapes), each test meets another outline, and the join is about 25 times slower on this sample (7 s against 0.3 s on the test server; compare the two with EXPLAIN ANALYZE, as in d3). So join from the polygons, summarize per record, and LEFT JOIN the records to that summary to keep the ones that matched nothing.
A point on an edge belongs to both neighbours under ST_Covers and ST_Intersects, and to neither under ST_Contains. No predicate gives it to exactly one; when each record must count once, the query needs a rule of its own, such as the smallest outline (the Exercise).
Two sources, neither the truth. Every GBIF record carries a country code (country_code_gbif), given by its publisher or worked out by GBIF from the coordinates. The outlines are Natural Earth's, simplified to 1:50m. Joining one to the other, on the whole sample, shows where they agree, where they disagree and which records fall in no outline. A disagreement is a question about the record or the map, not an answer.
Overlapping outlines count twice. g2 found eight pairs of outlines where one lies inside the other (France holds French Guiana, Réunion, Guadeloupe, Martinique and Mayotte). A record inside both matches both, so a count per country adds up to more than the records.

ST_Contains leaves out a point on a polygon's edge, ST_Covers and ST_Intersects keep it; where outlines overlap, a spatial join returns one record twice

Worked example

The three predicates on three points: inside a square, on its edge and outside it.

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

g3-boundary · PostgreSQL with PostGIS

WITH square AS (SELECT ST_GeomFromText('POLYGON((0 0, 1 0, 1 1, 0 1, 0 0))', 4326) AS geom),
points (place, geom) AS (
  VALUES ('inside', ST_SetSRID(ST_Point(0.5, 0.5), 4326)),
         ('on the edge', ST_SetSRID(ST_Point(1, 0.5), 4326)),
         ('outside', ST_SetSRID(ST_Point(2, 0.5), 4326))
)
SELECT p.place, ST_Contains(s.geom, p.geom) AS contains, ST_Covers(s.geom, p.geom) AS covers,
       ST_Intersects(s.geom, p.geom) AS intersects
FROM square AS s, points AS p
ORDER BY p.place;

-- result: g3-boundary · 3 rows · sha256 d10ead98357a
place | contains | covers | intersects
inside | true | true | true
on the edge | false | true | true
outside | false | false | false

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

Only the point on the edge tells the predicates apart: ST_Contains leaves it out, ST_Covers and ST_Intersects keep it. Had a second square shared that edge, ST_Covers would put the point in both squares and ST_Contains in neither.

Lab

Does the choice matter for the sample? Every record–outline pair that shares a point, counted under each predicate:
g3-edges · PostgreSQL with PostGIS
```
SELECT count(*) FILTER (WHERE ST_Contains(c.geom, o.geom)) AS contains,
       count(*) FILTER (WHERE ST_Covers(c.geom, o.geom)) AS covers,
       count(*) AS intersects
FROM country_shapes AS c
JOIN occurrences AS o ON ST_Intersects(c.geom, o.geom);
```
```
-- result: g3-edges · 1 row · sha256 72a27c571dbc
contains | covers | intersects
40329 | 40329 | 40329
```
In g3.sql, with the blocks it runs after.
The three counts are equal: no record lies exactly on an outline's edge, so here the predicates give the same join. They part where points and lines share a grid, such as records rounded to 0.1° against outlines drawn along whole degrees. ST_Covers is the usual choice for points in polygons: it loses no point, and a point on a shared edge shows up twice, where a count can catch it. The rest of this module uses it. The 40,329 pairs are more than the 39,396 records that fall in some outline; step 4 finds the extra 933.

Each record against GBIF's code: an outline it falls in is GBIF's country (agree), none of them is (disagree), or it falls in no outline. The join starts from the outlines, as the Concept says, and the LEFT JOIN comes after, on the per-record summary:

g3-agree · PostgreSQL with PostGIS

WITH hits AS (
  SELECT o.id, count(*) AS outlines, bool_or(c.country_code = o.country_code_gbif) AS agrees
  FROM country_shapes AS c
  JOIN occurrences AS o ON ST_Covers(c.geom, o.geom)
  GROUP BY o.id
)
SELECT CASE WHEN h.id IS NULL THEN 'no outline' WHEN h.agrees THEN 'agree' ELSE 'disagree' END AS outcome,
       count(*) AS records,
       count(*) FILTER (WHERE h.outlines > 1) AS in_two_outlines
FROM occurrences AS o
LEFT JOIN hits AS h ON h.id = o.id
GROUP BY 1
ORDER BY 1;

-- result: g3-agree · 3 rows · sha256 a48ca6bc5bd6
outcome | records | in_two_outlines
agree | 39190 | 932
disagree | 206 | 1
no outline | 999 | 0

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

Of the sample's 40,395 records, 39,190 agree, 206 disagree and 999 lie in no outline; every record has a GBIF code. bool_or makes a record agree when any outline it falls in is its code, so the 932 records in two outlines that agree are counted once, as agreeing.

Where they disagree, GBIF's code against the outline's:

g3-disagree · PostgreSQL with PostGIS

WITH hits AS (
  SELECT o.id, o.country_code_gbif AS gbif_code,
         string_agg(c.country_code, '+' ORDER BY c.country_code) AS outline_code,
         bool_or(c.country_code = o.country_code_gbif) AS agrees
  FROM country_shapes AS c
  JOIN occurrences AS o ON ST_Covers(c.geom, o.geom)
  GROUP BY o.id, o.country_code_gbif
)
SELECT gbif_code, outline_code, count(*) AS records
FROM hits
WHERE NOT agrees
GROUP BY gbif_code, outline_code
ORDER BY records DESC, gbif_code, outline_code
LIMIT 10;

-- result: g3-disagree · 10 rows · sha256 ba081d2e5cbb
gbif_code | outline_code | records
US | MX | 78
BJ | BF | 28
GF | BR | 18
EC | PE | 15
GF | SR | 12
… 5 more rows

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

The disagreements are pairs of neighbours: 78 records coded US inside Mexico's outline, 28 coded Benin inside Burkina Faso's, 18 coded French Guiana inside Brazil's, 15 coded Ecuador inside Peru's, 12 coded French Guiana inside Suriname's. Each is a record a few kilometres over a border, a border the 1:50m outline draws a few kilometres off, or a publisher's country that isn't where the point is; the query can't say which. Distance and overlay, g4, measures how close to a border such records lie.

The records inside two outlines at once:
g3-two-outlines · PostgreSQL with PostGIS
```
WITH hits AS (
  SELECT o.id, o.country_code_gbif AS gbif_code, string_agg(c.country_code, '+' ORDER BY c.country_code) AS outlines
  FROM country_shapes AS c
  JOIN occurrences AS o ON ST_Covers(c.geom, o.geom)
  GROUP BY o.id, o.country_code_gbif
  HAVING count(*) > 1
)
SELECT outlines, gbif_code, count(*) AS records
FROM hits
GROUP BY outlines, gbif_code
ORDER BY records DESC, outlines, gbif_code;
```
```
-- result: g3-two-outlines · 4 rows · sha256 f16f1fc4d7c9
outlines | gbif_code | records
FR+GF | GF | 923
FR+RE | RE | 6
FR+MQ | MQ | 3
FR+GF | SR | 1
```
In g3.sql, with the blocks it runs after.
933 records sit in two outlines (932 that agree, 1 that doesn't), every one inside France's and an overseas region's: French Guiana 924, Réunion 6, Martinique 3. Natural Earth's France includes its overseas regions, which also have outlines of their own (g2's g2-coverage), so a spatial join counts these records twice unless the outlines are a partition. The one coded Suriname sits inside French Guiana's outline, and is also one of the 206 disagreements.

Optional appendix: records against a range you downloaded yourself

IUCN's range maps may be downloaded with a free account for your own non-commercial use, and nothing made from them may be shared; critwire never holds them, so this block is never run here. Load your shapefile into a schema of your own, which the lab loader never touches (CREATE SCHEMA my_iucn, then QGIS, ogr2ogr or shp2pgsql into my_iucn.ranges, one multipolygon per taxon, the SIS taxon id in iucn_id). IUCN's files call that id id_no and hold several rows per taxon, by presence, origin and season: dissolve them with ST_Union, or a record counts once per row. The key is that id, which critwire publishes as species.iucn_id, never a name: IUCN renamed its name column in 2023. Records outside a range aren't errors by themselves: captive animals, misidentifications and incomplete ranges all put them there.

g3-iucn · PostgreSQL with PostGIS

SELECT s.common_name, count(*) AS records,
       count(*) FILTER (WHERE ST_Covers(r.geom, o.geom)) AS inside_range
FROM my_iucn.ranges AS r
JOIN species AS s ON s.iucn_id = r.iucn_id
JOIN occurrences AS o ON o.species_id = s.id
GROUP BY s.common_name
ORDER BY s.common_name;

Not run here: reads your own IUCN download, which critwire never holds.

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

Result

On the square, only ST_Covers and ST_Intersects count the point on the edge. On the sample the three predicates give the same 40,329 pairs, as no record lies on an edge. Of the 40,395 records, 39,190 agree with an outline, 206 disagree and 999 lie in no outline. The disagreements are border pairs (78 records coded US inside Mexico's outline, 28 coded Benin inside Burkina Faso's): a record a few kilometres over a border, a generalized border, or a publisher's country that isn't where the point is. 933 records sit in two outlines at once, every one inside France's and an overseas region's (French Guiana 924, Réunion 6, Martinique 3), so a per-country count adds up to 933 more than the records in outlines. Spatial indexes, g5, finds where the 999 in no outline are.

Exercise

Count the records in each country's outline, the usual spatial-join report. Which countries' counts are too high because outlines overlap, and by how much does the total exceed the records the outlines hold? Then give each record one country, the smallest outline that covers it, and count again.

Show answer

g3-ex · PostgreSQL with PostGIS

WITH outlines AS MATERIALIZED (
  SELECT country_code, geom, ST_Area(geom::geography) AS m2 FROM country_shapes
), hits AS (
  SELECT c.country_code, row_number() OVER (PARTITION BY o.id ORDER BY c.m2) AS rank
  FROM outlines AS c
  JOIN occurrences AS o ON ST_Covers(c.geom, o.geom)
)
SELECT CASE WHEN GROUPING(country_code) = 1 THEN 'all' ELSE country_code END AS country,
       count(*) AS every_outline, count(*) FILTER (WHERE rank = 1) AS smallest_outline
FROM hits
GROUP BY ROLLUP (country_code)
HAVING count(*) <> count(*) FILTER (WHERE rank = 1)
ORDER BY GROUPING(country_code), country_code;

-- result: g3-ex · 2 rows · sha256 3d08155603d3
country | every_outline | smallest_outline
FR | 944 | 11
all | 40329 | 39396

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

Check

The answer: one join from the outlines; row_number() ranks each record's outlines by area, smallest first, and ROLLUP adds a total row. Only the rows whose two counts differ are kept. The outline areas are worked out once, in a MATERIALIZED CTE.

Only France's count changes. Its outline holds 944 records, but 933 of them are also inside an overseas region's outline (step 4), so given one country each, France keeps 11: the 11 records GBIF codes FR. French Guiana, Réunion and Martinique keep their counts, since their outlines are the smaller ones. The total falls from 40,329 pairs to 39,396 records, 39,190 that agree and 206 that don't (step 2). The smallest outline is the right rule here because each overlap is a territory drawn inside its country (g2); for outlines that overlap along a border, no rule of area says which country a record is in. Without MATERIALIZED, PostgreSQL may inline the CTE and work out ST_Area(geom::geography) for every joined row, 40,329 times instead of 247: the same answer, two orders of magnitude slower (about 82 s against 0.45 s on the test server).

Dialects

PostgreSQL: ST_Contains, ST_Covers, ST_Intersects, ST_Within, ST_CoveredBy and ST_Relate with a pattern; each uses a GiST index on either argument by itself, through its bounding box. ST_ContainsProperly(a, b) also refuses a b that touches a's boundary at all, and can be answered faster.
DuckDB: the spatial extension (not in this course's browser) has the same predicates by the same names, and plans a join on one of them as a spatial join of its own, with no index to create. Without it, the browser's lab can join records to a box of degrees, as in Joins, f3, but not to an outline.
SQL Server: methods on the shape: a.STContains(b), a.STIntersects(b), b.STWithin(a), and a.STRelate(b, 'T*****FF*'); there is no STCovers, so STIntersects serves for points in polygons. A spatial index is used only when the method's result is compared with = 1 in WHERE or ON, as in ON c.geom.STContains(o.geom) = 1.