critwire · Learn · Learn's data

Rasters g7

Prerequisites: Grids (g6)

Read a raster in SQL: the value under a point with ST_Value, the cells inside a polygon with ST_Clip and ST_ValueCount, and areas from cell counts, with the error of the method that made the published habitat summaries. Everything runs in your own PostgreSQL with PostGIS, on the lab tables (the setup page).

Concept

How to read a raster in SQL, and how to turn its cells into an area you can defend:

The habitat map's cells narrow east–west away from the equator, and an outline counts each 100 m or 1 km cell whole or not at all, by its centre

Data

habitat_clips, the 100 m map over all of Costa Rica (cr) and Tasmania (tas), stored as in-database rasters in tiles of 256 × 256 cells, one row per tile, keyed (clip, tile) (the lab data). Each clip is its area's box, Tasmania's with Bruny and Maria Islands, and only the tiles holding some of the area are stored. Each clip's box and tile grid, the tiles kept and left out, and their cells, from the tiles' georeferences:

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

g7-tiles · PostgreSQL with PostGIS 

WITH box AS (
  SELECT clip, ST_Extent(ST_Envelope(rast)) AS ext, max(ST_ScaleX(rast)) AS d
  FROM habitat_clips
  GROUP BY clip
), grid AS (
  SELECT clip, round((ST_XMax(ext) - ST_XMin(ext)) / d)::integer AS width,
         round((ST_YMax(ext) - ST_YMin(ext)) / d)::integer AS height
  FROM box
)
SELECT g.clip, g.width || ' × ' || g.height AS box_cells,
       ceil(g.width / 256.0) || ' × ' || ceil(g.height / 256.0) AS box_tiles,
       count(*) AS tiles,
       ceil(g.width / 256.0) * ceil(g.height / 256.0) - count(*) AS left_out,
       sum(ST_Width(h.rast) * ST_Height(h.rast)) AS cells,
       sum(ST_Count(h.rast)) AS class_cells,
       ST_SameAlignment(h.rast) AS aligned
FROM habitat_clips AS h JOIN grid AS g USING (clip)
GROUP BY g.clip, g.width, g.height
ORDER BY g.clip;
-- result: g7-tiles · 2 rows · sha256 1514a2c9d499
clip | box_cells | box_tiles | tiles | left_out | cells | class_cells | aligned
cr | 3730 × 3474 | 15 × 14 | 117 | 93 | 7555072 | 7554437 | true
tas | 4131 × 3318 | 17 × 13 | 173 | 48 | 11206656 | 11206656 | true

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

Land lies in the outermost tiles on every side, so the kept tiles reach every edge of the box, and their extent is the box: 3,730 × 3,474 cells in 15 × 14 tiles for cr, 4,131 × 3,318 in 17 × 13 for tas. Tiles on the east and south edges are narrower or shorter, not padded. The 93 tiles left out of cr and the 48 left out of tas hold none of that area, only sea and, for cr, Nicaragua and Panama. The kept tiles hold 7,555,072 and 11,206,656 cells; 635 of cr's are nodata in the map itself, which ST_Count leaves out. ST_SameAlignment, used here as an aggregate, is true when every tile lies on one grid: the same cell size and skew, and corners whole cells apart. That is what lets ST_Union mosaic tiles without resampling a cell.

Most blocks below study two windows of 836 × 836 cells cut from them: cr, longitude −84.20…−83.45 and latitude 9.85…10.60, from San José to the Caribbean lowlands; tas, 147.25…148.00 and −43.10…−42.35, south-east Tasmania with Hobart. Each window is the 0.75° square, on a grid of 0.05° steps, holding the most of the sample's records for its clip. Also habitat_classes, country_habitats (method 3, every country), occurrence_habitats (each record's 1 km class), country_shapes and occurrences. Each windowed block names its window in win and builds it as one raster in w: ST_Clip cuts each stored tile to the window and ST_Union mosaics the pieces (ArcGIS Pro's Clip Raster, then Mosaic), and MATERIALIZED makes PostgreSQL build it once per query. ST_Clip's result covers every cell the polygon touches, so win's edges sit a fifth of a cell inside the window's outer cell edges: an edge on a cell boundary, or just past it, would add a column or row of nodata. The rasters are in the database, so no raster setting is needed: every block runs on a server's defaults.

Worked example

Each window's georeference from ST_MetaData, one cell per row built from it, and the cells' size and area on the ground against the map's 100 m:

g7-latitude · PostgreSQL with PostGIS 

WITH win (clip, env) AS (
  VALUES ('cr', ST_MakeEnvelope(-84.2007, 9.8493, -83.4501, 10.5999, 4326)),
         ('tas', ST_MakeEnvelope(147.2493, -43.1001, 147.9999, -42.3495, 4326))
), w AS MATERIALIZED (
  SELECT h.clip, ST_Union(ST_Clip(h.rast, win.env)) AS rast
  FROM habitat_clips AS h JOIN win ON win.clip = h.clip AND ST_Intersects(h.rast, win.env)
  GROUP BY h.clip
), meta AS (
  SELECT w.clip, m.*
  FROM w CROSS JOIN LATERAL ST_MetaData(w.rast) AS m
), cell_rows AS (
  SELECT clip, width, r.row = height / 2 AS is_mid,
         ST_MakeEnvelope(upperleftx, upperlefty + r.row * scaley,
                         upperleftx + scalex, upperlefty + (r.row - 1) * scaley, srid) AS cell
  FROM meta CROSS JOIN generate_series(1, height) AS r (row)
), mid AS (
  SELECT clip, ST_XMin(cell) AS x0, ST_XMax(cell) AS x1, ST_YMin(cell) AS y0, ST_YMax(cell) AS y1
  FROM cell_rows WHERE is_mid
)
SELECT c.clip,
       round(min(ST_YMin(c.cell))::numeric, 2) AS lat_south, round(max(ST_YMax(c.cell))::numeric, 2) AS lat_north,
       round(ST_Distance(ST_Point(m.x0, m.y0, 4326)::geography, ST_Point(m.x1, m.y0, 4326)::geography)::numeric, 1) AS mid_width_m,
       round(ST_Distance(ST_Point(m.x0, m.y0, 4326)::geography, ST_Point(m.x0, m.y1, 4326)::geography)::numeric, 1) AS mid_height_m,
       round(count(*) * max(c.width) * 0.01, 1) AS map_km2,
       round((sum(ST_Area(c.cell::geography)) * max(c.width) / 1e6)::numeric, 1) AS ground_km2,
       round((sum(ST_Area(c.cell::geography)) / (count(*) * 1e4))::numeric, 4) AS ground_over_map,
       round(cos(radians((m.y0 + m.y1) / 2))::numeric, 4) AS cos_mid_lat
FROM cell_rows AS c JOIN mid AS m USING (clip)
GROUP BY c.clip, m.x0, m.x1, m.y0, m.y1
ORDER BY c.clip;
-- result: g7-latitude · 2 rows · sha256 e063bd6e1c5b
clip | lat_south | lat_north | mid_width_m | mid_height_m | map_km2 | ground_km2 | ground_over_map | cos_mid_lat
cr | 9.85 | 10.60 | 98.4 | 99.4 | 6989.0 | 6834.8 | 0.9779 | 0.9841
tas | -43.10 | -42.35 | 73.6 | 99.8 | 6989.0 | 5131.4 | 0.7342 | 0.7346

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

Both windows hold 698,896 cells: 6,989.0 km² if every cell were 100 m by 100 m. On the ground the Costa Rica window covers 6,834.8 km² and the Tasmania window 5,131.4 km², a quarter less. A cell in the middle row is 98.4 m wide and 99.4 m tall at 10.2° N, and 73.6 m by 99.8 m at 42.7° S: the width follows cos φ (0.9841 and 0.7346). The height changes a little too, because the Earth is flattened: a degree of latitude is about 110.6 km near the equator and 111.1 km at 43°. So Costa Rica's ratio, 0.9779, sits a little below its cos φ. Within a window a cell's area changes from the north edge to the south by 0.2 % (9,768 to 9,790 m²) in Costa Rica and 1.2 % (7,298 to 7,386 m²) in Tasmania.

Lab

  1. The class under each record in the Costa Rica window: win and w build the window, as in every windowed block here, ST_Intersects keeps the records on it, and ST_Value reads the cell under each:

    g7-records · PostgreSQL with PostGIS 

    WITH win (clip, env) AS (
  VALUES ('cr', ST_MakeEnvelope(-84.2007, 9.8493, -83.4501, 10.5999, 4326))
), w AS MATERIALIZED (
  SELECT h.clip, ST_Union(ST_Clip(h.rast, win.env)) AS rast
  FROM habitat_clips AS h JOIN win ON win.clip = h.clip AND ST_Intersects(h.rast, win.env)
  GROUP BY h.clip
)
SELECT v.habitat_code, c.name, count(*) AS records
FROM w
JOIN occurrences AS o ON ST_Intersects(w.rast, o.geom)
CROSS JOIN LATERAL (SELECT ST_Value(w.rast, 1, o.geom)::integer AS habitat_code) AS v
LEFT JOIN habitat_classes AS c ON c.code = v.habitat_code
GROUP BY v.habitat_code, c.name
ORDER BY records DESC, v.habitat_code NULLS LAST;
    -- result: g7-records · 12 rows · sha256 e9484b7cc942
habitat_code | name | records
109 | Forest – Subtropical/tropical moist montane | 87
106 | Forest – Subtropical/tropical moist lowland | 30
307 | Shrubland – Subtropical/tropical high altitude | 23
507 | Wetlands (inland) – Permanent freshwater marshes/pools (unde… | 21
1401 | Arable Land | 13
… 7 more rows

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

    All 202 records in the window have a class: the window has no nodata cells. 87 are on moist montane forest, 30 on moist lowland forest, 23 on high-altitude shrubland, and 21 on small freshwater marshes, a class with 0.1 % of the window's land (step 3). Arable land, over a third of the cells, holds 13. Records gather where people go to look, so these counts describe the recording as much as the animals; Zonal statistics, s6, compares records with the classes' areas.

  2. The cells inside Costa Rica's outline, by class, as ST_ValueCount gives them:

    g7-valuecount · PostgreSQL with PostGIS 

    WITH win (clip, env) AS (
  VALUES ('cr', ST_MakeEnvelope(-84.2007, 9.8493, -83.4501, 10.5999, 4326))
), w AS MATERIALIZED (
  SELECT h.clip, ST_Union(ST_Clip(h.rast, win.env)) AS rast
  FROM habitat_clips AS h JOIN win ON win.clip = h.clip AND ST_Intersects(h.rast, win.env)
  GROUP BY h.clip
)
SELECT vc.value::integer AS habitat_code, vc.count AS cells
FROM w
JOIN country_shapes AS cs ON cs.country_code = 'CR'
CROSS JOIN LATERAL ST_ValueCount(ST_Clip(w.rast, cs.geom), 1) AS vc
ORDER BY cells DESC, habitat_code;
    -- result: g7-valuecount · 20 rows · sha256 f1cb0c12721a
habitat_code | cells
1401 | 253016
109 | 178008
106 | 102946
306 | 36961
1405 | 33473
… 15 more rows

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

    693,893 of the window's 698,896 cells have their centre inside the outline. The other 5,003 are the sea in its north-east corner, all of them marine neritic (900). Another 6,329 marine neritic cells lie inside the outline, along the coast: a 1:50m outline and a 100 m map don't draw the same coastline (Spatial joins, g3; Spatial indexes, g5). The 20 classes run from arable land, 253,016 cells, to tidepools, 13.

  3. The same cells as estimated areas. Every cell in a row has the same ground area, so count the cells by row, measure one cell per row on the ellipsoid, and multiply. MATERIALIZED makes the planner measure each row once, not once per class and row (g3's exercise):

    g7-estimate · PostgreSQL with PostGIS 

    WITH win (clip, env) AS (
  VALUES ('cr', ST_MakeEnvelope(-84.2007, 9.8493, -83.4501, 10.5999, 4326))
), w AS MATERIALIZED (
  SELECT h.clip, ST_Union(ST_Clip(h.rast, win.env)) AS rast
  FROM habitat_clips AS h JOIN win ON win.clip = h.clip AND ST_Intersects(h.rast, win.env)
  GROUP BY h.clip
), cells AS (
  SELECT p.val::integer AS habitat_code, p.y AS row
  FROM w
  JOIN country_shapes AS cs ON cs.country_code = 'CR'
  CROSS JOIN LATERAL ST_PixelAsCentroids(ST_Clip(w.rast, cs.geom), 1) AS p
), by_row AS (
  SELECT habitat_code, row, count(*) AS cells
  FROM cells
  GROUP BY habitat_code, row
), row_area AS MATERIALIZED (
  SELECT r.row, ST_Area(ST_PixelAsPolygon(w.rast, 1, r.row)::geography) / 1e6 AS cell_km2
  FROM w, (SELECT DISTINCT row FROM by_row) AS r
), est AS (
  SELECT b.habitat_code, sum(b.cells) AS cells, sum(b.cells * ra.cell_km2) AS est_km2
  FROM by_row AS b JOIN row_area AS ra USING (row)
  GROUP BY b.habitat_code
)
SELECT e.habitat_code, c.name, e.cells,
       round(e.cells * 0.01, 2) AS map_km2,
       round(e.est_km2::numeric, 2) AS est_km2,
       round((100 * e.est_km2 / sum(e.est_km2) OVER ())::numeric, 1) AS pct
FROM est AS e
LEFT JOIN habitat_classes AS c ON c.code = e.habitat_code
ORDER BY e.cells DESC, e.habitat_code;
    -- result: g7-estimate · 20 rows · sha256 9941c849a9ba
habitat_code | name | cells | map_km2 | est_km2 | pct
1401 | Arable Land | 253016 | 2530.16 | 2474.07 | 36.5
109 | Forest – Subtropical/tropical moist montane | 178008 | 1780.08 | 1741.50 | 25.7
106 | Forest – Subtropical/tropical moist lowland | 102946 | 1029.46 | 1006.14 | 14.8
306 | Shrubland – Subtropical/tropical moist | 36961 | 369.61 | 361.29 | 5.3
1405 | Urban Areas | 33473 | 334.73 | 327.59 | 4.8
… 15 more rows

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

    The window's land inside the outline is an estimated 6,785.90 km², from 836 area calls instead of 693,893. Each class's map area, cells × 0.01 km², is 2.3 % above its estimate: at 10° N a cell covers 0.978 of an equator cell (the worked example). Arable land is 36.5 % of the area, moist montane forest 25.7 % and moist lowland forest 14.8 %; the marsh class with 21 records, 0.1 %.

  4. A record's 100 m class against its 1 km class, for every record on both clips. occurrence_habitats gives each record the class with the largest fraction in its 1 km cell. Build that cell from the 1 km grid's numbers (upper-left corner −180.0099155291168, 90.00297126602864; cells 0.009915667822884942° on a side), and find the most common 100 m class inside it. A 1 km cell can straddle up to four tiles, so ST_Clip cuts each tile it touches and ST_Union mosaics the pieces, as w does for a window. A point on the edge between two tiles meets both, and its cell is the first of the tile to the east or south (a cell holds its west and north edges), so DISTINCT ON keeps each record once, on its highest-numbered tile:

    g7-scale · PostgreSQL with PostGIS 

    WITH grid (x0, y0, d) AS (
  VALUES (-180.0099155291168, 90.00297126602864, 0.009915667822884942)
), recs AS (
  SELECT DISTINCT ON (o.id) h.clip, ST_Value(h.rast, 1, o.geom)::integer AS class_100m, oh.habitat_code AS class_1km,
         ST_MakeEnvelope(g.x0 + k.col * g.d, g.y0 - (k.row + 1) * g.d, g.x0 + (k.col + 1) * g.d, g.y0 - k.row * g.d, 4326) AS cell_1km
  FROM occurrences AS o
  JOIN habitat_clips AS h ON ST_Intersects(h.rast, o.geom)
  JOIN occurrence_habitats AS oh ON oh.occurrence_id = o.id
  CROSS JOIN grid AS g
  CROSS JOIN LATERAL (SELECT floor((ST_X(o.geom) - g.x0) / g.d) AS col, floor((g.y0 - ST_Y(o.geom)) / g.d) AS row) AS k
  ORDER BY o.id, h.tile DESC
), majority AS (
  SELECT r.clip, r.class_100m, r.class_1km,
         (SELECT vc.value::integer FROM ST_ValueCount(cell.rast, 1) AS vc
          ORDER BY vc.count DESC, vc.value LIMIT 1) AS most_100m
  FROM recs AS r
  CROSS JOIN LATERAL (SELECT ST_Union(ST_Clip(h.rast, r.cell_1km)) AS rast
                      FROM habitat_clips AS h
                      WHERE h.clip = r.clip AND ST_Intersects(h.rast, r.cell_1km)) AS cell
)
SELECT clip, count(*) AS records,
       count(*) FILTER (WHERE class_100m = class_1km) AS same_at_record,
       count(*) FILTER (WHERE most_100m = class_1km) AS same_as_majority,
       count(*) FILTER (WHERE class_100m = most_100m) AS record_in_majority
FROM majority
GROUP BY clip
ORDER BY clip;
    -- result: g7-scale · 2 rows · sha256 dda12e4d18e0
clip | records | same_at_record | same_as_majority | record_in_majority
cr | 625 | 340 | 617 | 344
tas | 1991 | 1204 | 1921 | 1225

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

    The cr clip holds 625 records: 618 in Costa Rica, and 4 in Nicaragua and 3 in Panama on tiles that also hold some of Costa Rica. The record's own 100 m cell has its 1 km class for 340 of them (54 %) and 1,204 of the 1,991 in Tasmania (60 %). Yet the 1 km class is the most common 100 m class in its cell for 617 and 1,921: the two maps agree; most records just aren't on the majority. Only 344 of the Costa Rican records are; 78 sit on montane forest in 1 km cells that are mostly arable land. Of the 78 records whose 1 km class isn't the 100 m majority, 46 are in cells where two classes tie, which the block breaks toward the lower code (20 Tasmanian records share one cell with 48 cells each of temperate forest and pastureland), and 31 in cells where the majority leads by 9 cells or fewer of 121 to 132. The last is on the Tasmania clip's north edge, in a 1 km cell 14 % on the clip. A 1 km class describes the landscape around a record, a 100 m class the patch it is on; neither is the species' habitat.

  5. How good are the published summaries? Method 3 summed over each country, against the outline's area on the ellipsoid, in bands of error:

    g7-countries · PostgreSQL with PostGIS 

    WITH est AS (
  SELECT country_code, sum(est_area_km2) AS est_km2 FROM country_habitats GROUP BY country_code
), cmp AS (
  SELECT cs.country_code, ST_Area(cs.geom::geography) / 1e6 AS outline_km2, coalesce(e.est_km2, 0) AS est_km2
  FROM country_shapes AS cs LEFT JOIN est AS e USING (country_code)
), banded AS (
  SELECT *, abs(est_km2 - outline_km2) / outline_km2 AS err,
         CASE WHEN abs(est_km2 - outline_km2) / outline_km2 < 0.01 THEN 'under 1 %'
              WHEN abs(est_km2 - outline_km2) / outline_km2 < 0.05 THEN '1 to 5 %'
              WHEN abs(est_km2 - outline_km2) / outline_km2 < 0.2 THEN '5 to 20 %'
              ELSE '20 % or more' END AS error
  FROM cmp
)
SELECT coalesce(error, 'all') AS error, count(*) AS countries,
       round(max(outline_km2)::numeric, 2) AS largest_km2,
       round(sum(outline_km2)::numeric) AS outline_km2, round(sum(est_km2)) AS est_km2,
       round(sum(est_km2) / sum(outline_km2)::numeric, 6) AS est_over_outline,
       CASE WHEN GROUPING(error) = 0 THEN
         string_agg(country_code || ' ' || to_char(100 * (est_km2 - outline_km2) / outline_km2, 'FMSG990.0'), ', ' ORDER BY err DESC)
           FILTER (WHERE err >= 0.05) END AS worst_pct
FROM banded
GROUP BY ROLLUP (error)
ORDER BY error IS NOT NULL, min(err);
    -- result: g7-countries · 5 rows · sha256 96b0c186ffd8
error | countries | largest_km2 | outline_km2 | est_km2 | est_over_outline | worst_pct
all | 247 | 16973548.42 | 146715521 | 146579614 | 0.999074 | NULL
under 1 % | 215 | 16973548.42 | 146711158 | 146575271 | 0.999074 | NULL
1 to 5 % | 21 | 430.33 | 3957 | 3951 | 0.998359 | NULL
5 to 20 % | 8 | 134.94 | 390 | 375 | 0.961755 | SX -17.3, MF +10.6, CC -9.8, MO -6.2, PN +5.9, BL -5.8, VG -…
20 % or more | 3 | 10.21 | 16 | 17 | 1.050974 | VA -100.0, TV +28.9, TK -28.6

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

    Over all 247 outlines the estimate is 0.999074 of the outlines' area (−0.09 %), the figure critwire's ingest records with the files. 215 countries are within 1 %, every outline over 431 km² among them. The error grows as the outline shrinks, since a small island holds few 1 km cells and each counts whole or not at all: Vatican City (0.71 km²) holds no cell centre and has no rows; Tuvalu's atolls, 10.21 km², get 13.16 km² (+28.9 %); Tokelau −28.6 %, Sint Maarten −17.3 %. Read a small country's habitat areas as rough.

  6. The published summary against the 100 m map, for a whole country. Method 1 over all of Costa Rica: clip each tile the outline touches to it, count the cells by class and row, and give each row one geodesic area, as in step 3. A cell's area depends only on its latitude, so the block builds one cell per row at longitude 0. ST_Clip crops a tile to the outline, which moves its upper-left corner, so each row's centre latitude comes from the clipped tile's own georeference (ST_MetaData). Then set it beside country_habitats, method 3, class by class, and against the outline's area:

    g7-country · PostgreSQL with PostGIS 

    WITH cr AS (
  SELECT geom AS outline, ST_Area(geom::geography) / 1e6 AS outline_km2 FROM country_shapes WHERE country_code = 'CR'
), clipped AS (
  SELECT ST_Clip(h.rast, cr.outline) AS rast
  FROM habitat_clips AS h JOIN cr ON ST_Intersects(h.rast, cr.outline)
  WHERE h.clip = 'cr'
), by_row AS (
  SELECT p.val::integer AS habitat_code, m.upperlefty + (p.y - 0.5) * m.scaley AS lat, m.scalex AS d, count(*) AS cells
  FROM clipped AS c
  CROSS JOIN LATERAL ST_MetaData(c.rast) AS m
  CROSS JOIN LATERAL ST_PixelAsCentroids(c.rast, 1) AS p
  GROUP BY habitat_code, lat, d
), row_area AS MATERIALIZED (
  SELECT lat, ST_Area(ST_MakeEnvelope(0, lat - d / 2, d, lat + d / 2, 4326)::geography) / 1e6 AS cell_km2
  FROM (SELECT DISTINCT lat, d FROM by_row) AS r
), method1 AS (
  SELECT b.habitat_code, sum(b.cells * ra.cell_km2) AS km2
  FROM by_row AS b JOIN row_area AS ra USING (lat)
  GROUP BY b.habitat_code
), method3 AS (
  SELECT habitat_code, est_area_km2 AS km2 FROM country_habitats WHERE country_code = 'CR'
)
SELECT habitat_code, c.name,
       round(sum(m1.km2)::numeric, 2) AS method1_km2,
       round(sum(m3.km2), 2) AS method3_km2,
       round(100 * (sum(m3.km2) / sum(m1.km2)::numeric - 1), 1) AS diff_pct,
       round((100 * sum(m1.km2) / max(cr.outline_km2))::numeric, 2) AS pct_of_outline
FROM method1 AS m1
FULL JOIN method3 AS m3 USING (habitat_code)
LEFT JOIN habitat_classes AS c ON c.code = habitat_code
CROSS JOIN cr
GROUP BY ROLLUP ((habitat_code, c.name))
ORDER BY habitat_code IS NOT NULL, sum(m1.km2) DESC NULLS LAST, habitat_code;
    -- result: g7-country · 35 rows · sha256 8104db915c8f
habitat_code | name | method1_km2 | method3_km2 | diff_pct | pct_of_outline
NULL | NULL | 51297.88 | 51199.23 | -0.2 | 99.99
1401 | Arable Land | 16128.59 | 16108.04 | -0.1 | 31.44
109 | Forest – Subtropical/tropical moist montane | 13115.84 | 13095.13 | -0.2 | 25.57
106 | Forest – Subtropical/tropical moist lowland | 7127.03 | 7108.44 | -0.3 | 13.89
201 | Savanna - Dry | 3533.35 | 3528.48 | -0.1 | 6.89
… 30 more rows

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

    Over all of Costa Rica the two methods differ by 0.2 %: 51,199.23 km² by method 3 against 51,297.88 km² by method 1, whose cells cover 99.99 % of the outline's area. Every class over 100 km² agrees within 1 % except two marine ones, and the eight largest, each over 1,900 km², within 0.7 %. The marine classes come out larger by method 3: coral reef (908) by 6.4 %, seagrass (909) 1.6 % and marine neritic (900) 1.3 %. A 1 km cell on the coast whose centre is inside the outline counts whole, sea included: the border cells of the Concept. Rivers and streams (501, 502) come out 4.8 and 3.9 % smaller. Classes this small, a median of 5 and 3 cells of 100 m in each 1 km cell they touch, carry the largest relative differences, since a few cells either way is a large share of them. For a country this size, the published summary stands for the 100 m map to within 1 % for its main classes; read the small and coastal ones as rough.

Result

A cell of v004's 100 m grid covers 9,768 to 9,790 m² in the Costa Rica window and 7,298 to 7,386 m² in Tasmania's; counting cells as 100 m squares overstates by 2.3 % and 36 %. Costa Rica's window holds 202 records and an estimated 6,785.90 km² of land inside the outline, 36.5 % of it arable; its records sit mostly on forest (87 montane, 30 lowland). The clips are 290 aligned tiles, the ones holding some of Costa Rica or Tasmania. Over both clips a record's 100 m class matches its 1 km class for 54 % (Costa Rica) and 60 % (Tasmania) of records, although the 1 km class is the 100 m majority around all but 78 of the 2,616, nearly all of them ties or near ties. The published summaries total 0.999074 of the outlines' area; 215 of 247 countries are within 1 %, the worst are all under 135 km². For all of Costa Rica, method 3 is 0.2 % below method 1 at 100 m, within 1 % for every class over 100 km² but marine neritic and seagrass.

Exercise

Method 1 at 100 m is the reference for method 3. Compute both for the Tasmania window inside AU's outline, by class, from the 100 m cells alone. Give each 100 m cell its 1 km cell from the grid's numbers (step 4); take a class's fraction of a 1 km cell as its share of the 100 m centres in it; keep only the 1 km cells wholly inside the window; and count a 1 km cell when its centre is inside the outline. Which classes differ most, and why?

Show answer

g7-ex · PostgreSQL with PostGIS 

WITH win (clip, env) AS (
  VALUES ('tas', ST_MakeEnvelope(147.2493, -43.1001, 147.9999, -42.3495, 4326))
), w AS MATERIALIZED (
  SELECT h.clip, ST_Union(ST_Clip(h.rast, win.env)) AS rast
  FROM habitat_clips AS h JOIN win ON win.clip = h.clip AND ST_Intersects(h.rast, win.env)
  GROUP BY h.clip
), grid (x0, y0, d) AS (
  VALUES (-180.0099155291168, 90.00297126602864, 0.009915667822884942)
), au AS (
  SELECT geom AS outline FROM country_shapes WHERE country_code = 'AU'
), px AS (
  SELECT p.val::integer AS habitat_code, p.y AS row_100m,
         floor((ST_X(p.geom) - g.x0) / g.d)::integer AS col_1km, floor((g.y0 - ST_Y(p.geom)) / g.d)::integer AS row_1km,
         ST_Intersects(au.outline, p.geom) AS inside
  FROM w, au, grid AS g, ST_PixelAsCentroids(w.rast, 1) AS p
), area_100m AS MATERIALIZED (
  SELECT r.row_100m, ST_Area(ST_PixelAsPolygon(w.rast, 1, r.row_100m)::geography) / 1e6 AS km2
  FROM w, (SELECT DISTINCT row_100m FROM px) AS r
), cells_1km AS MATERIALIZED (
  SELECT k.col_1km, k.row_1km, ST_Area(k.cell::geography) / 1e6 AS km2,
         ST_Intersects(au.outline, ST_Centroid(k.cell)) AS centre_inside,
         ST_Covers(ST_Envelope(w.rast), k.cell) AS whole
  FROM (SELECT DISTINCT col_1km, row_1km FROM px) AS c
  CROSS JOIN grid AS g
  CROSS JOIN LATERAL (SELECT c.col_1km, c.row_1km,
                             ST_MakeEnvelope(g.x0 + c.col_1km * g.d, g.y0 - (c.row_1km + 1) * g.d,
                                             g.x0 + (c.col_1km + 1) * g.d, g.y0 - c.row_1km * g.d, 4326) AS cell) AS k
  CROSS JOIN w
  CROSS JOIN au
), counted AS (
  SELECT px.*, count(*) OVER w AS n_1km, count(*) FILTER (WHERE px.inside) OVER w AS n_inside,
         k.km2 AS km2_1km, k.centre_inside, a.km2 AS km2_100m
  FROM px JOIN cells_1km AS k USING (col_1km, row_1km) JOIN area_100m AS a USING (row_100m)
  WHERE k.whole
  WINDOW w AS (PARTITION BY px.col_1km, px.row_1km)
)
SELECT habitat_code,
       round(sum(km2_100m) FILTER (WHERE inside)::numeric, 2) AS method1_km2,
       round(sum(km2_1km / n_1km) FILTER (WHERE centre_inside)::numeric, 2) AS method3_km2,
       round((100 * (sum(km2_1km / n_1km) FILTER (WHERE centre_inside) / sum(km2_100m) FILTER (WHERE inside) - 1))::numeric, 1) AS diff_pct,
       round((100 * coalesce(sum(km2_100m) FILTER (WHERE inside AND n_inside < n_1km), 0)
                  / sum(km2_100m) FILTER (WHERE inside))::numeric, 1) AS border_pct
FROM counted
GROUP BY ROLLUP (habitat_code)
HAVING sum(km2_100m) FILTER (WHERE inside) >= 1
ORDER BY habitat_code NULLS FIRST;
-- result: g7-ex · 13 rows · sha256 53600dc8b50d
habitat_code | method1_km2 | method3_km2 | diff_pct | border_pct
NULL | 3953.67 | 3944.68 | -0.2 | 2.7
104 | 1536.61 | 1536.27 | 0.0 | 0.9
304 | 739.01 | 738.71 | 0.0 | 2.1
305 | 10.53 | 10.57 | 0.4 | 0.0
404 | 57.52 | 57.78 | 0.4 | 5.6
… 8 more rows

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

Check

The answer: assign the cells, count each 1 km cell's centres with a window function, and sum both ways side by side. The share of each class's method 1 area that lies in 1 km cells the outline cuts tells you why:

Over the window's land the methods differ by 0.2 %: 3,944.68 km² by method 3 against 3,953.67 km² by method 1. The classes that differ most, by 1.1 to 2.8 %, are intertidal (1200), urban areas (1405), inland wetlands (500), marine neritic (900) and seagrass (909): 5 to 15 % of their area lies in 1 km cells the outline crosses, where method 3 takes the whole cell or none of it. Classes away from the coast agree within 0.4 %, and what's left there is the grids' misfit: a 1 km cell spans 11.04 cells of 100 m, so a whole one holds 121, 132 or 144 of their centres, and sharing its area among them is close, not exact. An earlier comparison over all of Trinidad and Tobago and of Åland, on 1:10m outlines, found the same pattern, larger: marine classes +12.5 % and +13.9 %, with three quarters of their area in border cells. The comparison judges the summary's method, not the map: both columns are estimates from the same classification.

Dialects

Lab schema: the PostGIS tables

critwire Learn lab database, PostGIS schema learn_curriculum

Every column, in words · The diagram alone