critwire · Learn · Learn's data

Distance and overlay g4

Prerequisites: Spatial joins (g3)

Measure in metres and cut one layer by another: areas and lengths on the ellipsoid, ST_DWithin and buffers, ST_Intersection and ST_Difference, and how far two separately drawn maps are from fitting. Everything runs in your own PostgreSQL with PostGIS, on the lab tables (the setup page).

Concept

How to measure lengths, areas and distances in metres, and how to cut one layer by another:

Worked example

How much of Costa Rica does each ecoregion cover? Cut each ecoregion that touches the outline by it, measure the piece on the ellipsoid, and divide by the outline's area:

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

g4-ecoregions · PostgreSQL with PostGIS 

WITH cr AS (SELECT geom, ST_Area(geom::geography) AS m2 FROM country_shapes WHERE country_code = 'CR')
SELECT e.eco_id, e.name, b.name AS biome,
       round((ST_Area(ST_Intersection(e.geom, cr.geom)::geography) / 1e6)::numeric) AS km2,
       round((100 * ST_Area(ST_Intersection(e.geom, cr.geom)::geography) / cr.m2)::numeric, 1) AS pct
FROM ecoregions AS e
JOIN cr ON ST_Intersects(e.geom, cr.geom)
LEFT JOIN biomes AS b ON b.biome_id = e.biome_id
ORDER BY km2 DESC, e.eco_id;
-- result: g4-ecoregions · 7 rows · sha256 ed3884acd3d8
eco_id | name | biome | km2 | pct
470 | Isthmian-Atlantic moist forests | Tropical & Subtropical Moist Broadleaf Forests | 17288 | 33.7
471 | Isthmian-Pacific moist forests | Tropical & Subtropical Moist Broadleaf Forests | 9012 | 17.6
458 | Costa Rican seasonal moist forests | Tropical & Subtropical Moist Broadleaf Forests | 8400 | 16.4
506 | Talamancan montane forests | Tropical & Subtropical Moist Broadleaf Forests | 8311 | 16.2
527 | Central American dry forests | Tropical & Subtropical Dry Broadleaf Forests | 6062 | 11.8
… 2 more rows

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

Seven ecoregions reach into Costa Rica, led by the Isthmian-Atlantic moist forests (17,288 km², 33.7 %). Five are moist or dry broadleaf forest; the other two are mangroves along the coasts, Southern Mesoamerican Pacific (637 km²) and Mesoamerican Gulf-Caribbean (198 km²). ST_Intersects in the ON clause picks the ecoregions, with the index; ST_Intersection in the SELECT cuts only those seven. The shares add up to 97.3 %, not 100: the Lab asks why.

Lab

  1. Do the pieces add up to the country? The outline's area, the pieces' sum, and the part of the outline that no ecoregion covers (the outline less the union of the ecoregions):

    g4-gap · PostgreSQL with PostGIS 

    WITH cr AS (SELECT geom FROM country_shapes WHERE country_code = 'CR')
SELECT round((ST_Area(cr.geom::geography) / 1e6)::numeric) AS country_km2,
       round((sum(ST_Area(ST_Intersection(e.geom, cr.geom)::geography)) / 1e6)::numeric) AS ecoregion_pieces_km2,
       round((ST_Area(ST_Difference(cr.geom, ST_Union(e.geom))::geography) / 1e6)::numeric) AS in_no_ecoregion_km2
FROM cr
JOIN ecoregions AS e ON ST_Intersects(e.geom, cr.geom)
GROUP BY cr.geom;
    -- result: g4-gap · 1 row · sha256 5811a1ff6de4
country_km2 | ecoregion_pieces_km2 | in_no_ecoregion_km2
51303 | 49908 | 2070

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

    The pieces add up to 49,908 km² of the outline's 51,303, and 2,070 km² lie in no ecoregion. Had the pieces been a partition, the two would add up to the outline exactly; they add up to 51,978, about 675 km² too many. So the pieces both leave gaps and overlap. About half of the gap lies within 5 km of the sea, where the ecoregions stop short of the outline's coast; the rest is slivers inland, where two neighbouring ecoregions both pulled back from their shared border.

  2. Where do the pieces overlap? Each pair of pieces whose interiors share an area (g2's '2********' test), with its DE-9IM matrix, and a total row first:

    g4-overlaps · PostgreSQL with PostGIS 

    WITH cr AS (SELECT geom FROM country_shapes WHERE country_code = 'CR'),
pieces AS (
  SELECT e.eco_id, ST_Intersection(e.geom, cr.geom) AS geom
  FROM ecoregions AS e
  JOIN cr ON ST_Intersects(e.geom, cr.geom)
),
shared AS (
  SELECT a.eco_id || '+' || b.eco_id AS pair, ST_Relate(a.geom, b.geom) AS de9im,
         ST_Area(ST_Intersection(a.geom, b.geom)::geography) AS m2
  FROM pieces AS a
  JOIN pieces AS b ON a.eco_id < b.eco_id AND ST_Relate(a.geom, b.geom, '2********')
)
SELECT CASE WHEN GROUPING(pair) = 1 THEN 'all' ELSE pair END AS pair, de9im,
       round((sum(m2) / 1e6)::numeric, 1) AS overlap_km2
FROM shared
GROUP BY ROLLUP ((pair, de9im))
ORDER BY GROUPING(pair) DESC, overlap_km2 DESC, pair;
    -- result: g4-overlaps · 11 rows · sha256 b0c2c01abe86
pair | de9im | overlap_km2
all | NULL | 678.2
470+506 | 212111212 | 181.7
458+527 | 212111212 | 143.8
458+506 | 212111212 | 76.3
458+471 | 212101212 | 62.0
… 6 more rows

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

    Ten pairs overlap, 678.2 km² in all: the 675 km² that step 1's pieces counted twice (a few square kilometres lie in three ecoregions at once, and come into three pairs). The largest are between neighbours along long borders: the Isthmian-Atlantic moist forests and the Talamancan montane forests (181.7 km²), the seasonal moist and the dry forests (143.8). Each is a thin strip along a border that each ecoregion, simplified on its own, drew in a slightly different place. Compare the matrices with g2's 2F2F11FF2, a territory wholly inside its country: here the third and seventh characters are 2 too, so each piece also reaches outside the other. These are partial overlaps between shapes that are each valid: the layer isn't a coverage, and nothing about one row says so. Pieces that overlap are counted twice by any sum over them, as step 1's was; a partition (each overlap given to one ecoregion, each gap to the nearest) would need rules of its own, and the honest report is the overlay's error: here about 1.3 % of the country in overlaps and 4.0 % in gaps.

  3. Costa Rica's land borders, their length, and the records within 10 km of each, by the code GBIF gave them:

    g4-border · PostgreSQL with PostGIS 

    WITH cr AS (SELECT geom FROM country_shapes WHERE country_code = 'CR'),
border AS (
  SELECT n.country_code AS neighbour,
         ST_LineMerge(ST_CollectionExtract(ST_Intersection(ST_Boundary(cr.geom), ST_Boundary(n.geom)), 2)) AS geom
  FROM cr
  JOIN country_shapes AS n ON n.country_code <> 'CR' AND ST_Touches(n.geom, cr.geom)
)
SELECT b.neighbour,
       round((ST_Length(b.geom::geography) / 1000)::numeric) AS border_km,
       count(o.id) AS records_within_10km,
       count(o.id) FILTER (WHERE o.country_code_gbif = 'CR') AS coded_cr,
       count(o.id) FILTER (WHERE o.country_code_gbif = b.neighbour) AS coded_neighbour
FROM border AS b
LEFT JOIN occurrences AS o
  ON o.country_code_gbif IN ('CR', b.neighbour) AND ST_DWithin(o.geom::geography, b.geom::geography, 10000)
GROUP BY b.neighbour, b.geom
ORDER BY b.neighbour;
    -- result: g4-border · 2 rows · sha256 539a1bf4a653
neighbour | border_km | records_within_10km | coded_cr | coded_neighbour
NI | 267 | 5 | 1 | 4
PA | 250 | 9 | 6 | 3

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

    ST_Touches finds the outlines that share a border with Costa Rica's but no interior: Nicaragua and Panama. The border is where their two boundaries meet: ST_Intersection of the boundaries returns the shared lines (and any lone points where they only touch), ST_CollectionExtract(…, 2) keeps the lines and ST_LineMerge joins them into one. The land borders are 267 km (Nicaragua) and 250 km (Panama). Within 10 km of them, on either side, lie 5 and 9 records coded for Costa Rica or the neighbour; 7 of the 14 are coded for the neighbour. So close to a line, the code and the outline can each be right or wrong: a 1:50m border is generalized and can be kilometres off, and the disagreements in Spatial joins, g3, are mostly records like these. The ON clause casts every record's geom to geography, so the index on occurrences.geom can't help, and the join reads all 40,395 records for each border (about 0.7 s); Spatial indexes, g5, shows the plan and an index that does help.

Result

Overlaying the ecoregions on Costa Rica gives seven pieces, led by the Isthmian-Atlantic moist forests (17,288 km², 33.7 %). The pieces add up to 49,908 km² of the outline's 51,303, with 2,070 km² in no ecoregion and about 675 km² counted twice: ten pairs of valid pieces overlap, 678.2 km² in all, along borders each ecoregion drew in its own place. A valid layer isn't a coverage, and an area from this overlay carries an error of a few per cent of the country. Costa Rica's land borders are 267 km (Nicaragua) and 250 km (Panama); 14 records lie within 10 km of them, half of them coded for the neighbour.

Exercise

How much of Costa Rica lies within 10 km of each land border, and which records are in that strip? Draw the strip with ST_Buffer on geography, cut it by Costa Rica's outline, and give its area, its records, how many of them GBIF codes for the neighbour, and, as a check on the buffer, how many of them ST_DWithin finds within 10 km of the border.

Show answer

g4-ex · PostgreSQL with PostGIS 

WITH cr AS (SELECT geom FROM country_shapes WHERE country_code = 'CR'),
border AS (
  SELECT n.country_code AS neighbour,
         ST_LineMerge(ST_CollectionExtract(ST_Intersection(ST_Boundary(cr.geom), ST_Boundary(n.geom)), 2)) AS geom
  FROM cr
  JOIN country_shapes AS n ON n.country_code <> 'CR' AND ST_Touches(n.geom, cr.geom)
),
strip AS (
  SELECT b.neighbour, b.geom AS line,
         ST_Intersection(ST_Buffer(b.geom::geography, 10000)::geometry, cr.geom) AS geom
  FROM border AS b, cr
)
SELECT s.neighbour,
       round((ST_Area(s.geom::geography) / 1e6)::numeric) AS strip_km2,
       count(o.id) AS records,
       count(o.id) FILTER (WHERE o.country_code_gbif = s.neighbour) AS coded_neighbour,
       count(o.id) FILTER (WHERE ST_DWithin(o.geom::geography, s.line::geography, 10000)) AS within_10km
FROM strip AS s
LEFT JOIN occurrences AS o ON ST_Covers(s.geom, o.geom)
GROUP BY s.neighbour, s.geom
ORDER BY s.neighbour;
-- result: g4-ex · 2 rows · sha256 e63983707c3d
neighbour | strip_km2 | records | coded_neighbour | within_10km
NI | 2673 | 2 | 1 | 2
PA | 2289 | 7 | 1 | 7

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

Check

The answer: g4-border's border lines, buffered by 10,000 m on geography, cast back to geometry and cut by the outline (ST_Intersection); then a spatial join of the records to each strip, as in g3.

The strips are 2,673 km² along Nicaragua and 2,289 km² along Panama, 5.2 % and 4.5 % of Costa Rica. A strip is close to its border's length times 10 km (2,670 and 2,500 km²): the buffer's round ends add a little, and the sea at the border's ends and the border's bends, where the buffer overlaps itself, take some away. The strips hold 2 and 7 records, against g4-border's 5 and 9 on both sides of the line: the other 3 and 2 lie outside Costa Rica's outline. One record in each strip is coded for the neighbour although it lies inside Costa Rica's outline, two of g3's 206 disagreements. within_10km equals records: the buffer's straight segments lost no record here. Had the cut been left out, the strip would reach 10 km into each neighbour as well, and a record there would count for Costa Rica.

Dialects

Lab schema: the PostGIS tables

critwire Learn lab database, PostGIS schema learn_curriculum

Every column, in words · The diagram alone