critwire · Learn · Learn's data

Zonal statistics s6

Prerequisites: Hot and cold spots (s5), Rasters (g7)

Do some zones hold more records than their size predicts, and what is left once recording effort is accounted for? Records and area by ecoregion in SQL and in GeoPandas, a chi-square test, a ratio normalized for effort, and the same comparison by the classes of a global habitat map, on Costa Rica's records. Everything runs in your own PostgreSQL with PostGIS and the Python environment, on the lab tables (the setup page).

Concept

Point patterns, s2, to Hot and cold spots, s5, asked where records cluster; this module asks whether some zones hold more records than they should, and what "should" can mean.

Worked example

Costa Rica's ecoregions as zones: area, records (those GBIF codes CR that lie inside the outline, g3), records per 1,000 km², the record-to-area ratio, and the jaguar's effort-normalized ratio:

Checked 2026-11-02 on PostgreSQL 16.15, PostGIS 3.4.2, GEOS 3.12.1, GDAL 3.8.4; Python 3.12.3 with requirements.txt sha256 d5bbe28024cc.

s6-zones · PostgreSQL with PostGIS 

WITH country AS (
  SELECT geom FROM country_shapes WHERE country_code = 'CR'
), zones AS (
  SELECT e.eco_id, e.name, ST_Intersection(e.geom, country.geom) AS geom
  FROM ecoregions AS e, country
  WHERE ST_Intersects(e.geom, country.geom)
), records AS (
  SELECT o.id, o.species_id, o.geom
  FROM occurrences AS o, country
  WHERE o.country_code = 'CR' AND ST_Contains(country.geom, o.geom)
), per_zone AS (
  SELECT z.eco_id, z.name, ST_Area(z.geom::geography) / 1e6 AS km2,
         count(r.id) AS records,
         count(r.id) FILTER (WHERE r.species_id = (SELECT id FROM species WHERE scientific_name = 'Panthera onca')) AS jaguar
  FROM zones AS z
  LEFT JOIN records AS r ON ST_Contains(z.geom, r.geom)
  GROUP BY z.eco_id, z.name, z.geom
)
SELECT eco_id, name, round(km2::numeric) AS km2, records, jaguar,
       round((1000 * records / km2)::numeric, 1) AS per_1000_km2,
       round(((records / sum(records) OVER ()) / (km2 / sum(km2) OVER ()))::numeric, 2) AS record_ratio,
       round(((jaguar::numeric / nullif(records, 0)) / (sum(jaguar) OVER () / sum(records) OVER ()))::numeric, 2) AS jaguar_ratio
FROM per_zone
ORDER BY km2 DESC, eco_id;
-- result: s6-zones · 7 rows · sha256 72c259664d1d
eco_id | name | km2 | records | jaguar | per_1000_km2 | record_ratio | jaguar_ratio
470 | Isthmian-Atlantic moist forests | 17288 | 215 | 0 | 12.4 | 1.23 | 0.00
471 | Isthmian-Pacific moist forests | 9012 | 96 | 0 | 10.7 | 1.05 | 0.00
458 | Costa Rican seasonal moist forests | 8400 | 107 | 0 | 12.7 | 1.26 | 0.00
506 | Talamancan montane forests | 8311 | 13 | 0 | 1.6 | 0.15 | 0.00
527 | Central American dry forests | 6062 | 62 | 61 | 10.2 | 1.01 | 7.29
… 2 more rows

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

Read it one CTE at a time. zones cuts each ecoregion that touches the outline to the outline (ST_Intersection, Distance and overlay, g4). records keeps the records coded CR that lie inside the outline. per_zone joins records to zones with ST_Contains and counts them, with LEFT JOIN and count(r.id) so a zone with no record gets 0, not no row; the jaguar's records are counted with FILTER. ST_Area on geography measures on the WGS 84 ellipsoid, in m². The last query divides each zone's share by the column totals, sum(…) OVER (); nullif keeps a zone with no record from dividing by zero.

591 of Costa Rica's 619 records lie inside its outline (g3 found the rest). By area alone the records are uneven: the Talamancan montane forests, 8,311 km² of mountains, hold 13 records (ratio 0.15), the three lowland moist forests 1.05 to 1.26. The jaguar's column looks different again: 61 of the 62 records in the Central American dry forests and all 7 in the Southern Mesoamerican Pacific mangroves are jaguars, 7.29 and 7.41 times its share of the records in zones, and it has none in the other five.

Lab

  1. The same zones in GeoPandas: an overlay of the ecoregions with the outline, then a spatial join of the records to the pieces, with areas in an equal-area projection (EPSG:6933, the global EASE-Grid 2.0), and a chi-square test of records against area. from_wkb reads the shapes PostGIS sends as well-known binary (Loading and validating shapes, g2):

    s6-geopandas · Python

    Runs after: py-connect

    import geopandas as gpd
import shapely
from scipy.stats import chisquare

with connect() as conn:
    outline = conn.execute("SELECT ST_AsBinary(geom) FROM country_shapes WHERE country_code = 'CR'").fetchone()[0]
    eco = conn.execute(
        """SELECT e.eco_id, e.name, ST_AsBinary(e.geom)
           FROM ecoregions AS e JOIN country_shapes AS c ON c.country_code = 'CR' AND ST_Intersects(e.geom, c.geom)
           ORDER BY e.eco_id"""
    ).fetchall()
    recs = conn.execute("SELECT id, longitude, latitude FROM occurrences WHERE country_code = 'CR' ORDER BY id").fetchall()

country = gpd.GeoDataFrame(geometry=[shapely.from_wkb(bytes(outline))], crs=4326)
zones = gpd.GeoDataFrame({"eco_id": [e[0] for e in eco], "name": [e[1] for e in eco]},
                         geometry=shapely.from_wkb([bytes(e[2]) for e in eco]), crs=4326)
pieces = gpd.overlay(zones, country, how="intersection", keep_geom_type=True)
points = gpd.GeoDataFrame({"id": [r[0] for r in recs]},
                          geometry=gpd.points_from_xy([r[1] for r in recs], [r[2] for r in recs]), crs=4326)
inside = points[points.within(country.geometry.iloc[0])]
counts = gpd.sjoin(inside, pieces, predicate="within").groupby("eco_id").size()

pieces["km2"] = pieces.to_crs(6933).area / 1e6
pieces["records"] = pieces["eco_id"].map(counts).fillna(0).astype(int)
pieces = pieces.sort_values(["km2", "eco_id"], ascending=[False, True])
print(f"records coded CR: {len(points)}  inside the outline: {len(inside)}  in a zone: {pieces['records'].sum()}")
print("eco_id | km2 | records")
for z in pieces.itertuples():
    print(f"{z.eco_id} | {z.km2:.0f} | {z.records}")
expected = pieces["records"].sum() * pieces["km2"] / pieces["km2"].sum()
test = chisquare(pieces["records"], expected)
print(f"chi-square: {test.statistic:.1f}  df: {len(pieces) - 1}  p: {test.pvalue:.2g}")
    -- result: s6-geopandas · python · sha256 1c182b76d465
records coded CR: 619  inside the outline: 591  in a zone: 504
eco_id | km2 | records
470 | 17288 | 215
471 | 9012 | 96
458 | 8400 | 107
506 | 8311 | 13
527 | 6061 | 62
617 | 636 | 7
613 | 198 | 4
chi-square: 77.4  df: 6  p: 1.2e-14

    In s6.py and s6.ipynb, with the blocks it runs after.

    The course's test checks that every zone s6-zones shows has the same records here, and an area within 1 km²: geography's ellipsoid and EPSG:6933's projection differ by less than that (6,062 against 6,061 km² for the dry forests). "In a zone" is 504, but 3 records lie where two pieces overlap (g4) and are counted in both, so 501 records are in a zone and 90 in the gaps the simplified ecoregions leave inside the outline; both engines count a record once per piece that contains it. The chi-square test rejects "records follow area": 77.4 on 6 degrees of freedom, p 1.2 × 10⁻¹⁴.

  2. Records by habitat class against the classes' estimated areas, for the whole country. Each record takes its class from occurrence_habitats; each class's area is country_habitats' estimate:

    s6-habitat · PostgreSQL with PostGIS 

    WITH records AS (
  SELECT oh.habitat_code, count(*) AS records
  FROM occurrences AS o
  JOIN occurrence_habitats AS oh ON oh.occurrence_id = o.id
  WHERE o.country_code = 'CR'
  GROUP BY oh.habitat_code
)
SELECT ch.habitat_code, hc.name, ch.est_area_km2, coalesce(r.records, 0) AS records,
       round((coalesce(r.records, 0) / sum(coalesce(r.records, 0)) OVER ())
             / (ch.est_area_km2 / sum(ch.est_area_km2) OVER ()), 2) AS record_ratio
FROM country_habitats AS ch
LEFT JOIN records AS r ON r.habitat_code = ch.habitat_code
LEFT JOIN habitat_classes AS hc ON hc.code = ch.habitat_code
WHERE ch.country_code = 'CR'
ORDER BY ch.est_area_km2 DESC, ch.habitat_code;
    -- result: s6-habitat · 34 rows · sha256 3494dd1c6e04
habitat_code | name | est_area_km2 | records | record_ratio
1401 | Arable Land | 16108.04 | 236 | 1.21
109 | Forest – Subtropical/tropical moist montane | 13095.13 | 31 | 0.20
106 | Forest – Subtropical/tropical moist lowland | 7108.44 | 128 | 1.49
201 | Savanna - Dry | 3528.48 | 35 | 0.82
105 | Forest – Subtropical/tropical dry | 2236.93 | 42 | 1.55
… 29 more rows

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

    Arable land, 31 % of the estimated area, holds 236 of the 619 records (ratio 1.21); 216 of them are red-eyed tree frogs, recorded in 1 km cells that are mostly farmland: the cell's class, not the pond the frog sat in. Moist montane forest, a quarter of the area, holds 31 (0.20), and plantations 2 on 2,107 km² (0.08). The highest ratios belong to small classes at the coast: Marine Neritic, the sea over the continental shelf, has 59 records on 514 km² (9.50), and its coral reefs 7 on 29 km² (20.01). The 1:50m outline takes in some coastal water, and a 1 km cell that is mostly sea gets a marine class even when the record was made on the shore; 54 of the 59 lie inside the outline, up to 3.3 km from its edge, and 5 outside it. Such a ratio measures where the map's coastline, the outline and the records' coordinates disagree, not where frogs, sloths and jaguars live. The block also counts the 28 records outside the outline against areas inside it: the exercise fixes that.

Result

Costa Rica's 591 records inside the outline are spread unevenly over its ecoregions: the montane forests hold about a seventh of what their area predicts and the lowland forests a little more, and the chi-square test rejects records in proportion to area (77.4, 6 degrees of freedom, p 1.2 × 10⁻¹⁴). PostGIS and GeoPandas count the same records in every zone, with areas within 1 km². The jaguar's records sit almost only in the dry forests and the Pacific mangroves, 7.29 and 7.41 times its share. By habitat class, arable land holds 1.21 times its share of the records and the coastal marine classes 9.50 and 20.01. All of these are ratios of recording: they say where people recorded, and where one collar dataset's animals walked, not which forest a jaguar prefers.

Exercise

s6-habitat counts every record coded CR, including the 28 outside the outline, against areas inside it. Compare like with like: keep only the records inside the outline, as s6-zones does, and normalize the jaguar for effort by habitat class, as s6-zones does by ecoregion. For each class with a record: its records, the jaguar's records, the number of datasets the jaguar's records come from, and the jaguar's effort-normalized ratio. In which classes is the jaguar recorded, and from how many datasets?

Show answer

s6-ex · PostgreSQL with PostGIS 

WITH country AS (
  SELECT geom FROM country_shapes WHERE country_code = 'CR'
), records AS (
  SELECT oh.habitat_code, o.dataset_key,
         o.species_id = (SELECT id FROM species WHERE scientific_name = 'Panthera onca') AS jaguar
  FROM occurrences AS o
  JOIN country ON ST_Contains(country.geom, o.geom)
  JOIN occurrence_habitats AS oh ON oh.occurrence_id = o.id
  WHERE o.country_code = 'CR'
), per_class AS (
  SELECT habitat_code, count(*) AS records,
         count(*) FILTER (WHERE jaguar) AS jaguar,
         count(DISTINCT dataset_key) FILTER (WHERE jaguar) AS jaguar_datasets
  FROM records
  GROUP BY habitat_code
)
SELECT p.habitat_code, hc.name, p.records, p.jaguar, p.jaguar_datasets,
       round((p.jaguar::numeric / p.records) / (sum(p.jaguar) OVER () / sum(p.records) OVER ()), 2) AS jaguar_ratio
FROM per_class AS p
JOIN habitat_classes AS hc ON hc.code = p.habitat_code
ORDER BY p.jaguar DESC, p.records DESC, p.habitat_code;
-- result: s6-ex · 14 rows · sha256 fb2d0c58861d
habitat_code | name | records | jaguar | jaguar_datasets | jaguar_ratio
105 | Forest – Subtropical/tropical dry | 42 | 42 | 1 | 6.29
201 | Savanna - Dry | 35 | 35 | 1 | 6.29
900 | Marine Neritic | 54 | 12 | 1 | 1.40
1401 | Arable Land | 236 | 5 | 1 | 0.13
106 | Forest – Subtropical/tropical moist lowland | 108 | 0 | 0 | 0.00
… 9 more rows

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

Check

The outline drops 28 records: Marine Neritic goes from 59 to 54, its coral reefs from 7 to 4, and moist lowland forest from 128 to 108. The jaguar's 94 records are all inside the outline and fall in four classes. Every record in dry forest (42) and in dry savanna (35) is a jaguar, 6.29 times its share of 94 in 591; it has 12 of the 54 in Marine Neritic (1.40) and 5 of the 236 on arable land (0.13), and none anywhere else. jaguar_datasets is 1 in each: all 94 come from the Jaguar Movement Database, GPS collars worn by jaguars tracked in the dry north-west (s1). Normalizing for effort removes the other species' recording, not the collar dataset's: a ratio of one species from one dataset says where the collared animals walked, in a habitat map's classes at 1 km, not which habitat jaguars prefer. count(DISTINCT dataset_key) FILTER (WHERE jaguar) counts the datasets of one kind of row inside a group, as s5's exercise counted rows.

Dialects

Lab schema: the PostGIS tables

critwire Learn lab database, PostGIS schema learn_curriculum

Every column, in words · The diagram alone