16c: a land-cover class per triangle, natural colours

ROADMAP.md

@@ -46,7 +46,7 @@ increment that most needs a picture to check against
 | — | Release hardening: measure libc++'s `_LIBCPP_HARDENING_MODE_FAST` (and libstdc++'s assertions on the GCC leg), which bounds-check `std::vector`, `span` and the like in Release, on the 1 m benchmark; switch it on if the cost is small. Today an out-of-range index in shipped code the tests miss is undefined behaviour and crashes Python (Ola, 2026-09-27: "I'm surprised we don't have proper memory control"; CI's ASan/UBSan/TSan cover what the tests reach) | to measure (Ola, 2026-09-27); 21d, which it was placed after, is deferred behind the basin work | none yet |
 | 15 | A DEM in several tiles, and the domain in its own CRS: 15a Norway, many tiles in one projected CRS (Ola's 254 DTM10 UTM33 tiles; only the selected tiles must share a lattice, since 8 of them sit half a cell off); 15b the domain polygon reprojected into the DEM's CRS; 15c and 15d the São Francisco basin (geographic DEMs meshed in the DEM's own lattice frame, not resampled; window decoding and memory) | **15a shipped with branch `increment15-dem-mosaic`'s PR**: `--dem DIR` or several files, `--bbox`, tiles selected and stitched on one lattice, Q2-Q5 refusals; overlaps that disagree (real DTM10 tiles exported on different dates, up to 52 m) are split down the middle and reported per seam (Ola's Q1 revised, 2026-09-28). The design's acceptance box, 9 tiles and 10,051² nodes, meshes: 11.05 M triangles in 31 s on battery. **15b built with branch `increment15b-domain-crs`'s PR**: `--domain` in its own CRS (e.g. EPSG:4326) over a UTM33 mosaic, reprojected vertex by vertex. Then 16b (Ola, 2026-09-27: Norway first) | `docs/increments/15-dem-mosaic.md` |
 | 16b | Interior polygons and polylines as constraints ("terrain polygons": lakes, land cover, roads, rivers): `--features PATH`, a GeoJSON `FeatureCollection`, each feature naming a vocabulary property; closed or open `Breakline`s, crossings noded, off-node vertices with bilinear z. **Its working example is real data** (Ola, 2026-09-27): CORINE Land Cover 2018 over the benchmark tile `7908_3_10m_z33.tif`. The source is Ola's local copy, `rasputin_data/corine_sql/.../U2018_CLC2018_V2020_20u1.gpkg` (8.2 GB, EPSG:3035, a sibling of this repository, which also holds the 254-tile DTM10 archive for gap 6). It reads without GDAL: sqlite3 over its R-tree, the GeoPackage blob header stripped, `shapely.wkb`, then `pyproj` to 25833, so CRS stops in Python as before. Probed 2026-09-27 in 0.2 s: 60 polygons in 8 classes (heath, bare rock, sparse vegetation, bogs, intertidal flats, water, sea, urban), 11 068 vertices clipped to the tile, median segment 54 m against 10 m cells. The EEA's public ArcGIS service (`image.discomap.eea.europa.eu`, `Corine/CLC2018_WM`) returns the same 11 068 clipped vertices and is the route for anyone without the file. What it forces on 16b's design: neighbouring polygons share their boundaries, so each shared edge arrives twice; the polygons run past the domain and must be clipped; the extract is committed as a fixture with the Copernicus attribution. Placed before 20c because 20c may split constraint segments and should be designed and measured on inputs that have interior ones | **designed** (`@architect`, 2026-09-28): two PRs, 16b-0 (the noder's verifier indexed: it is quadratic, 25 s of a 30 s CORINE run) and 16b-1/2 (GeoJSON and GeoPackage features, class maps, clip as linework, the noder merges shared edges); Q1-Q6 ruled by Ola 2026-09-28 (Q6: the legacy CORINE GML stays and is read). **16b-0 shipped (#107, merged 2026-09-28)**: the verifier's pair search by sort and sweep, always on; `node` on a 48 km CORINE square 24.8 s → 0.067 s. **16b-1/2 built on branch `increment16b12-features`** (2026-09-28): GeoPackage, GeoJSON and the legacy GML (strict reader, fixture repaired) as `--features`; CORINE class maps; the pre-clip keeps whole edges and widens its margin per edge for long ones; the noder merges shared edges; net +642 production lines. Before 20c | `docs/increments/16b-terrain-polygons.md`, `docs/increments/16-domain-polygon.md` (R6) |
-| 16c | A land-cover label per triangle: which input polygon (CORINE class) each triangle lies in, by a flood fill over the unconstrained edges. The bits on an edge say what kind of line it is; which class lies on each side is 16c's | planned: ruled by Ola 2026-09-28 (16b's Q2) as a separate increment, next after 16b | `docs/increments/16b-terrain-polygons.md` (Q2, R6) |
+| 16c | A land-cover label per triangle: which input polygon (CORINE class) each triangle lies in, by a flood fill over the unconstrained edges. The bits on an edge say what kind of line it is; which class lies on each side is 16c's | **built** (2026-09-29, overnight; net 249 lines, one PR; Bygdin at 10 m: every triangle labelled, class shares equal to CORINE's clipped areas at two decimals). Designed by `@architect`, 2026-09-29: components across unconstrained edges, one point-in-polygon test per component, cell array `land_cover_code`, a natural-colour ParaView preset from `rasputin palette corine`; defaults D1-D7 for Ola to confirm; ~230 production lines. Ruled by Ola 2026-09-28 (16b's Q2) as a separate increment | `docs/increments/16c-landcover-labels.md`, `docs/increments/16b-terrain-polygons.md` (Q2, R6) |
 | 16d | A record of the input polylines in the mesh. Today a constraint edge keeps only its type bits; the source feature, its attributes (a river's name or order, a road's class) and which edges belong to one polyline are lost, and where edges from two features merge only the union of their bits survives. Proposed: a feature table (one row per input feature: source ID and attributes) and, per constraint edge, the list of features it came from (a list, since a merged edge belongs to several). The noder already tracks each piece's input chain; the open part is the output: a per-edge list in `.vtk`/`.ply`, or a side file | proposed by Ola 2026-09-28; to be designed by `@architect` once Ola has said what it is for; after 16b | - |
 | 20c | Soft quality criterion: a penalty that each Steiner node or constraint split must pay for in angle gained, instead of 20's hard 25°; applied at the start and during DEM refinement; may split constraint segments when that improves the mesh. Ola's rulings on 20's C1-C3 | to design after 16b (`@architect` measures cost against 20 first) | `docs/increments/20-start-quality.md` (Ola's rulings) |
 | — | Auto-catchment: the watershed upstream of a coordinate, computed from the DEM and handed to `--domain`, so a catchment no longer has to be supplied as a file (Ola, 2026-09-27: "not far into the future"). The textbook route is depression handling (Priority-Flood, Barnes, Lehman and Mulla 2014), D8 flow directions (O'Callaghan and Mark 1984) and accumulation, the pour point snapped to the strongest flow nearby, the upstream cells traced and their outline turned into a polygon; the literature check is `@architect`'s. Open for its design: whether it runs in the C++ core (a 10 m tile is 25 M cells); how a stair-stepped cell outline becomes a domain polygon, which meets input coarsening; and that a real catchment crosses tile edges, so it needs gap 6 (a DEM in several tiles) first. Legacy has nothing on it (`grep -rliE "watershed|flow.?acc|flow.?dir|pour.?point|catchment" legacy` returns no files) | to design; after gap 6, which it needs; placed after 20c, can move ahead of it on Ola's word | none yet |

docs/benchmarks/2026-09-29/bygdin-landcover.md

@@ -0,0 +1,28 @@
+# Increment 16c acceptance on Bygdin: summary (@perf, 2026-09-29)
+
+**Verdict: ACCEPTED** on the design's three pass criteria. AC power, Apple M1
+Max, commit `f7d5f14`. Method, tables and scripts:
+`2026-09-29/bygdin-landcover/README.md`.
+
+Reduced Bygdin catchment (from 22), DTM10, CORINE 2018 Norway extract with
+`--features-map corine`:
+
+| criterion | 10 m | 1 m | pass |
+|---|---|---|---|
+| `vtkPolyDataReader` reads it; `land_cover_code` on every cell, 0 on every line | 92,053 cells | 1,159,508 cells | yes |
+| seven codes, no code 0, each share within 0.01 pp of CORINE clipped | largest gap 0.000001 pp | the same | yes |
+| design's centroid oracle, triangles checked / disagreements | 83,171 / 0 | 1,145,434 / 0 | yes |
+| I1: unconstrained edges with differing codes | 0 of 116,398 | 0 of 1,705,117 | yes |
+
+Shares: 333 42.34 %, 332 21.32 %, 322 17.00 %, 512 16.46 %, 335 2.44 %,
+412 0.33 %, 142 0.10 %, equal to 22's table.
+
+Time, not gated (median of 3): the `land cover` phase takes 0.071 s at 10 m
+(83,171 triangles, total 2.08 s) and 1.296 s at 1 m (1,145,434 triangles,
+total 3.89 s). At 1 m that is more than `refine` (0.559 s), which is the
+condition the design's R2 sets for moving `regions` into the core.
+
+No bench or thread sweep: 16c's commits change no C++ (design, Acceptance).
+Pictures: `bygdin_landcover_oblique.png` (from the south-east, 2× vertical)
+and `bygdin_landcover_top.png`. ParaView preset: `corine_natural.json`. The
+manual ParaView import (design item 6) is left for Ola.

docs/benchmarks/2026-09-29/bygdin-landcover/README.md

@@ -0,0 +1,142 @@
+# Increment 16c acceptance: Bygdin in natural colours (@perf, 2026-09-29)
+
+The summary and the verdict are in `../bygdin-landcover.md`. This file has the
+method, the tables and where every number comes from. The design's acceptance
+section (`docs/increments/16c-landcover-labels.md`, "Acceptance") names the
+directory `16c-bygdin/`; this run was asked to use `bygdin-landcover/`.
+
+## What was run
+
+- **Commit**: `f7d5f14` on `increment16c-landcover-labels` (`logs/commit.txt`).
+  The 16c commits (`196147e` to `f7d5f14`) change nothing under `include/`,
+  `src/`, `tests/cpp/`, `tools/` or `CMakeLists.txt`
+  (`git diff --stat 196147e~1..f7d5f14 -- include src tests/cpp CMakeLists.txt tools`
+  is empty), so the bench and thread sweep are not required (design,
+  Acceptance).
+- **Build**: `cmake --build build-pyext -j --target _core` (exit 0), copied into
+  `.venv`; the installed and built `_core` sha256 are equal
+  (`logs/commit.txt`). Python is the editable `src_python/` tree.
+- **Machine**: Apple M1 Max, 10 cores, 32 GB. `mesh` used 10 threads.
+- **Power**: AC throughout. `pmset -g batt` at the start, after the timed
+  runs, and at the end (`logs/pmset_*.txt`): each reads "AC Power, 80 %, AC
+  attached".
+- **Inputs**:
+  - domain `bygdin_reduced_t20.geojson`, committed here, copied with
+    `git show increment22-autocatchment:docs/benchmarks/2026-09-29/bygdin/bygdin_reduced_t20.geojson`
+    (sha256 `8083fb5b...`, the same as 22's);
+  - DEM `../rasputin_data/DTM10_UTM33_20260925` (10 m, EPSG:25833);
+  - features `../rasputin_data/corine2018_dtm10_utm33.gpkg`, layer
+    `corine2018`, `--features-map corine`.
+
+Scripts, all in this directory:
+
+- `run.sh`: three timed `--binary` runs at `--tolerance 10` and at
+  `--tolerance 1`, then one `--ascii` run of each for the quality check, then
+  `rasputin palette corine --out corine_natural.json`. Every command runs under
+  `/usr/bin/time -l`; logs and `--stats` reports go to `logs/`.
+- `analyse.py`: every table below (`logs/analysis.txt` is its output). It reads
+  the meshes with VTK's `vtkPolyDataReader`. The CORINE reference is read
+  straight from the GeoPackage with `sqlite3` and shapely, not through
+  `tin_engine.feature_input`, so it does not share code with what it checks.
+  The quality check calls `tools/bench.py`'s `quality()` and
+  `read_vtk_ascii()` unchanged.
+- `render.py`: the two pictures.
+
+Meshes stay out of the repository. To regenerate: from the repository root
+with `.venv` active, `SCRATCH=<dir> bash docs/benchmarks/2026-09-29/bygdin-landcover/run.sh`,
+then `python .../analyse.py <dir>` and
+`python .../render.py <dir>/lc_t10_run1.vtk docs/benchmarks/2026-09-29/bygdin-landcover`.
+
+## 1. The file opens in VTK with `land_cover_code` on every cell
+
+| run | cells | lines | triangles | `land_cover_code` values | nonzero on lines | `FieldData` `land_cover_codes` |
+|---|---|---|---|---|---|---|
+| 10 m | 92,053 | 8,882 | 83,171 | 92,053 | 0 | present |
+| 1 m | 1,159,508 | 14,074 | 1,145,434 | 1,159,508 | 0 | present |
+
+The field text reads `CORINE Land Cover level-3 code, attribute Code_18, map
+corine; 0 = in no polygon, and every constraint line`. The active `SCALARS`
+stays `feature_mask` (13, ruling 5). The ASCII run's codes equal the binary
+run's, cell for cell, at both tolerances.
+
+## 2. Class shares against CORINE clipped to the catchment
+
+79 CORINE polygons meet the domain; clipped to it they cover 304.909550 km²,
+the domain's area. The mesh's plan area is 304.909551 km² at both tolerances.
+Triangle area (x, y) per code:
+
+| code | class | mesh km² | mesh share | CORINE clipped share | difference (pp) | 22's table |
+|---|---|---|---|---|---|---|
+| 333 | sparsely vegetated areas | 129.093253 | 42.33821 % | 42.33821 % | -0.000001 | 42.34 % |
+| 332 | bare rocks | 65.011569 | 21.32159 % | 21.32159 % | +0.000000 | 21.32 % |
+| 322 | moors and heathland | 51.849126 | 17.00476 % | 17.00476 % | +0.000000 | 17.00 % |
+| 512 | water bodies | 50.191119 | 16.46099 % | 16.46099 % | +0.000000 | 16.46 % |
+| 335 | glaciers and perpetual snow | 7.452202 | 2.44407 % | 2.44407 % | -0.000000 | 2.44 % |
+| 412 | peat bogs | 1.020748 | 0.33477 % | 0.33477 % | +0.000000 | 0.33 % |
+| 142 | sport and leisure facilities | 0.291534 | 0.09561 % | 0.09561 % | -0.000000 | 0.10 % |
+
+The 10 m and 1 m rows agree to the printed digits; the table is both. The
+same seven codes, no triangle with code 0, and the largest difference is
+0.000001 percentage points (limit 0.01).
+
+**The check can fail.** Relabelling the 142 region as 333 on the 10 m mesh
+moves 333's share to 42.43383 %, 0.096 pp off, and the oracle below then
+disagrees on 115 triangles (a one-off probe, not committed).
+
+## 3. I1 (spread) and I2 (oracle)
+
+| run | interior unconstrained edges | I1: edges with differing codes | I2: triangles checked (r > 3 mm) | I2: disagreements |
+|---|---|---|---|---|
+| 10 m | 116,398 | 0 | 83,171 of 83,171 | 0 |
+| 1 m | 1,705,117 | 0 | 1,145,434 of 1,145,434 | 0 |
+
+The oracle is the design's (R1): each triangle's centroid against the CORINE
+polygons (read independently, unclipped), smallest area wins; margin
+`2 × 1 mm`, so the exact set is `r > 3 mm`, which is every triangle here.
+`--stats` stderr, every run: `land cover: 124 regions, 0 outside every
+polygon, 0 in more than one, 0 thinner than the snap`.
+
+## 4. Time (recorded, not gated)
+
+Median of three binary runs; phase times from each run's `--stats` report,
+wall and RSS from `/usr/bin/time -l`.
+
+| tolerance | triangles | components | union-find rounds | land cover (s) | refine (s) | features clip (s) | total (s) | wall (s) | max RSS |
+|---|---|---|---|---|---|---|---|---|---|
+| 10 m | 83,171 | 124 | 5 | 0.071 (0.075, 0.071, 0.070) | 0.047 | 1.460 | 2.075 | 2.34 | 719 MB |
+| 1 m | 1,145,434 | 124 | 7 | 1.296 (1.303, 1.296, 1.263) | 0.559 | 1.447 | 3.887 | 4.36 | 859 MB |
+
+"Union-find rounds" counts the hook rounds of `landcover.regions`' outer loop:
+`analyse.py` runs a copy of that loop with a counter and checks that its
+result equals `regions()` (it does, at both tolerances).
+
+At 1 m the `land cover` phase (1.30 s) costs more than `refine` (0.56 s).
+The design's R2 names this as the condition for moving `regions` into the
+core. Where inside the phase the time goes was not profiled.
+
+## 5. Quality (the ASCII runs, `tools/bench.py`'s `quality()`)
+
+| tolerance | worst angle | max degree | max error | within tolerance | Delaunay checked / ambiguous / violations |
+|---|---|---|---|---|---|
+| 10 m | 0.00408° | 16 | 9.99973 m | yes | 116,398 / 1,095 / 0 |
+| 1 m | 0.00408° | 29 | 1.0 m | yes | 1,705,117 / 162,591 / 0 |
+
+The 10 m figures equal 22's "reduced + CORINE" row (83,171 triangles, 0.00408°,
+degree 16). 22 has no 1 m run with CORINE features; its 1 m run without them
+had max degree 18. Labelling does not touch the mesh, so these are the
+refine's figures on this input, not 16c's.
+
+## 6. The pictures and the ParaView preset
+
+- `bygdin_landcover_oblique.png`: the 10 m mesh from the south-east, z
+  exaggerated 2×; `bygdin_landcover_top.png`: map view. `render.py` draws the
+  triangles only (the constraint lines, code 0, are left out), coloured through
+  an indexed `vtkLookupTable` built from `tin_engine.palettes.CORINE_NATURAL`,
+  with a legend of the seven codes present.
+- `corine_natural.json`: written by `rasputin palette corine --out`.
+  **In ParaView**: colour by `land_cover_code` (cells), then Colour Map
+  Editor → Choose Preset → Import `corine_natural.json`, select "rasputin
+  CORINE natural", Apply (tick *Interpret Values As Categories* if it is not
+  on).
+- Not done here (design item 6, manual, for Ola): the ParaView import itself,
+  its version, and whether the preset switches to categorical by itself.