16b-1/2: terrain polygons as constraints (GeoPackage, GeoJSON, GML; CORINE maps)

ROADMAP.md

@@ -45,7 +45,9 @@ increment that most needs a picture to check against
 | — | The serial phase: profile refine's serial insert-and-flip phase, then parallelise what the profile blames. Scaling tops out at about 2.0-2.2×. Profiled 2026-09-27: serial part about a third of 1-thread refine, mostly Lawson legalisation; the scan stops near 5× from load imbalance | profiled; designed as increment 21 (Ola's rulings 2026-09-27: L1 determinism, one path, at most 2 % more triangles). **21a shipped with branch `increment21a-quick-wins`'s PR**: dynamic scan blocks, active merge, reused flip stack; mesh bit-identical; on AC refine -10 to -15 % at 8 threads, ceiling 2.1x -> 2.3x (`docs/benchmarks/2026-09-27/21a-acceptance.md`). **21b shipped with branch `increment21b-lattice-incircle`'s PR**: an int64 lattice incircle answers 99.97-100 % of refine's incircle tests; mesh bit-identical; on battery refine -15 % at 8 threads, ceiling 2.3x -> 2.5x (`docs/benchmarks/2026-09-27/21b-acceptance.md`). **21c measured** (`docs/benchmarks/2026-09-27/21c/README.md`): option C costs 6.5-10 % more triangles; A1 about 1 %; A0 is bit-identical (on these two inputs, not proven), and with evaluate-once is modelled at only 7-11 % faster refine at 8 threads and slower at 4. **21d is deferred** (Ola, 2026-09-27: "Review and push 21c, then basin-work"), behind the work in "Order of work" below; domain decomposition parallelises scan and split together | `docs/increments/21-parallel-refine.md`, `docs/benchmarks/2026-09-27/serial-profile/README.md` |
 | — | 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 with branch `increment16b0-verifier-sweep`'s PR**: 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 next; before 20c | `docs/increments/16b-terrain-polygons.md`, `docs/increments/16-domain-polygon.md` (R6) |
+| 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) |
+| 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 |
 | — | `raster/`: grid-to-world geometry and bilinear sampling | shipped (`7785fea`), **no record** | none — predates the protocol |

docs/benchmarks/2026-09-28/16b0-acceptance/README.md

@@ -267,3 +267,7 @@ angle, in degrees; "worst" is the smallest in the mesh, in degrees.
   repo root: `.venv/bin/python docs/benchmarks/2026-09-28/16b0-acceptance/corine/scripts/drv.py
   --pkg <tree>/build-bench/pkg mesh --side 48 --tol 10 --features A --min-angle 25`
   (the committed `drv.py` also looks for `clc_*.py` beside itself).
+- `clc_mesh.py:62` and `drv.py:120` monkeypatch `cli._domain_chains`, which
+  16b-1/2 replaced with `chains.start_chains` (`5079da8`), so `drv.py mesh
+  --features A` and `clc_mesh.py` reproduce only with `<tree>` checked out at
+  16b-0's `18be867`.

docs/benchmarks/2026-09-28/16b12-acceptance.md

@@ -0,0 +1,76 @@
+# Increment 16b-1/2 acceptance (@perf, 2026-09-28/29): summary
+
+**Verdict (2026-09-28), superseded by the 2026-09-29 addendum (ACCEPTED): NOT ACCEPTED as designed, on one measure: the candidate read from
+Ola's European GeoPackage.** It takes **17.5-19.9 s**, against the
+design's admitted ~1 s (R5). Everything else passes. All runs were on
+**battery** and are compared against battery. The branch is `5f3a522`; the
+base is `origin/master` `14f5fe3`. Method, tables and raw data are in
+`16b12-acceptance/README.md`.
+
+- **The defect, with its cause measured.** `io/geopackage.py`
+  `query_features` joins the feature table to the R-tree. SQLite plans that
+  as a scan of all 2.4 M feature rows, with one R-tree probe per row
+  (`EXPLAIN QUERY PLAN`: `SCAN t`). The query takes 17.5-18.7 s. The
+  R-tree scan alone, which is what the design timed, takes 0.71 s. The same
+  345 rows fetched by primary key through an `IN` subquery take 0.85 s. On
+  top of that, `features clip` takes 8.4-8.6 s. It is thought to be the
+  reprojection of the 14.7 M candidate vertices, not profiled. 13.8 M of
+  those vertices are in rows that only the widening admits, and all of them
+  are dropped. With the Norway file (25833), read and clip take 0.9 s and
+  1.3 s. Reported, not fixed.
+- **The quarter circle (M5), committed extract, `--features-map corine`:**
+  M5 is reproduced.
+  - Start vertices: 5 627, as in M5.
+  - Triangles: 441 855 at 1 m (M5 441 863) and 54 840 at 10 m (equal).
+  - Tolerance held, with 0 nodes uncovered.
+  - `vtkPolyDataReader` reads the file with R10's fields.
+  - 8 242 of 9 345 constraint edges carry `land_cover` and 3 905 carry
+    `water`.
+  - `node` 0.013 s, against M5's 0.33 s.
+- **Ola's case:** a 287 km² ring in EPSG:4326 on the 4-tile corner of
+  `6602_1`/`6602_2`/`6603_3`/`6603_4` in DTM10_UTM33_20260925, run with
+  both GeoPackages.
+  - **The two routes give bit-identical meshes:** 1 235 226 triangles at
+    1 m and 63 020 at 10 m.
+  - stderr: `87 features kept, 258 (Europe) / 175 (Norway) dropped outside,
+    28 clipped, 0 empty skipped`, then `11316 input vertices, 5642 noded
+    vertices`. There is no scan line.
+  - At 1 m the whole run takes 30.3 s with the Europe file, 6.3 s with the
+    Norway file, and 3.4 s without features.
+- **The 1 m benchmark and thread sweep:** 4 pairs, no change.
+  - The `_core` sha256 and both mesh sha256s are identical in all 8 runs.
+  - Pooled over 42 cells, the change is a median of +0.27 % (range −3.3 %
+    to +3.4 %), with no cell above +5 %.
+  - bench.py flagged 12 quarter cells in pair 3. The same build moves up to
+    +13.6 % between its own runs, and the fourth pair, run as a tie-break,
+    was ACCEPTED.
+  - Ceiling 2.5-2.8x, as the base.
+- **CORINE baseline, M4's 48 km square, through the real CLI** (Europe file):
+  - It holds against 16b-0: 6 753 084 triangles at 1 m (16b-0: 6 752 517),
+    and 480 941 at 10 m (480 961), which is still **3.61x** the
+    featureless mesh. That is 20c's input.
+  - `node` is 0.15 s, because the CLI feeds the rings unmerged (16b-0's
+    layout B: 0.14 s).
+  - Refine takes 5.9 s at 1 m. Peak RSS is 3.1 GB, which includes writing
+    the file and `--stats`; 16b-0's scripts did neither.
+  - Feature input is 26 s of 48 s: the same defect.
+- **Also found:** `tools/bench.py:605` fails when `--label` has a `/` and
+  `--mesh-dir` is given, because the mesh's parent directory is never
+  created. The workaround was to create it first.
+- **Not measured:** a profile of `features clip`, and any AC run.
+
+## Addendum, 2026-09-29: re-timed after `d58d693` (battery, 64 → 63 %)
+
+`features read` from the Europe file, before → after `d58d693`'s
+primary-key fetch, median of 3 runs:
+
+- Catchment at 1 m: 17.90 → 1.77 s.
+- Catchment at 10 m: 19.08 → 0.98 s.
+- 48 km square at 1 m: 17.72 → 1.01 s.
+
+The totals fall from 30.3 to 14.0 s, from 28.2 to 9.9 s and from 48.0 to
+30.6 s. The three meshes are byte-identical to the earlier runs.
+**Verdict: the read defect is fixed; 16b-1/2 is ACCEPTED.**
+`features clip` (8.2-8.5 s) is a known open item, pending Ola's decision,
+and is not counted against the verdict. Details are in the README's
+addendum.

docs/benchmarks/2026-09-28/16b12-acceptance/16b12-r2/README.md

@@ -0,0 +1,79 @@
+# Benchmark run `16b12-acceptance/16b12-r2`
+
+Generated by `tools/bench.py`.
+
+## Verdict
+
+- ACCEPTED
+
+Time threshold: 5.0 %.
+
+## Method
+
+- Started 2026-09-28T23:43:01.153728+02:00; tree `5f3a52291ef7eebe7092f0b0f0c4da02ec77b879` (dirty); bench.py blob `77765b181714c31714da476c2335c799c6a84021`. Build: Release, AppleClang 21.0.0.21000101, `-O3 -DNDEBUG`, _core sha256 `190587e7ba0ea96a74a8621767d7d20257c3968126a710db1364e1195afb27d2`.
+- Apple M1 Max, 8 P + 2 E cores, 32 GiB, macOS 27.0, Python 3.14.7, numpy 2.5.3.
+- Power **battery** (79%), `pmset -g batt` before and after (in run.json).
+- DEM `/Users/skavhaug/projects/rasputin/tests/fixtures/dem_archive/7908_3_10m_z33.tif` (sha256 `aabd0cbc28471ce8593e4c278381811bec3fbf4e1194058b3b4888d00c3af575`), tolerance 1.0, extra mesh args ``; domains: `quarter` (/Users/skavhaug/projects/rasputin/docs/benchmarks/2026-09-26/quarter.geojson), `tile` (the whole tile).
+- One child per sample, `/Users/skavhaug/projects/rasputin/.venv/bin/python /Users/skavhaug/projects/rasputin/tools/bench.py _child --pkg /Users/skavhaug/projects/rasputin/build-bench/pkg --threads <N> -- mesh --dem /Users/skavhaug/projects/rasputin/tests/fixtures/dem_archive/7908_3_10m_z33.tif --tolerance 1 --out <path> --binary`, repeats interleaved over thread counts; t=0 is the CLI's default, other counts are forced into `refine`. Quality: one `--ascii` run per domain at t=0, kept out of the repository; rerun the child with `--ascii --out PATH` to regenerate it.
+
+## `quarter`
+
+| threads | n | median s | min s | max s |
+|---:|---:|---:|---:|---:|
+| 0 | 5 | 0.1661 | 0.1617 | 0.1678 |
+| 1 | 5 | 0.4117 | 0.4048 | 0.4321 |
+| 2 | 5 | 0.2673 | 0.2618 | 0.2724 |
+| 3 | 5 | 0.2205 | 0.2179 | 0.2261 |
+| 4 | 5 | 0.1993 | 0.1978 | 0.1995 |
+| 5 | 5 | 0.1854 | 0.1824 | 0.1998 |
+| 6 | 5 | 0.1738 | 0.1731 | 0.1996 |
+| 7 | 5 | 0.1666 | 0.1659 | 0.1700 |
+| 8 | 5 | 0.1620 | 0.1605 | 0.1637 |
+| 9 | 5 | 0.1657 | 0.1605 | 0.1746 |
+| 10 | 5 | 0.1641 | 0.1629 | 0.1661 |
+| 11 | 5 | 0.1645 | 0.1585 | 0.1651 |
+| 12 | 5 | 0.1637 | 0.1611 | 0.1772 |
+| 13 | 5 | 0.1646 | 0.1586 | 0.1679 |
+| 14 | 5 | 0.1648 | 0.1628 | 0.1755 |
+| 15 | 5 | 0.1658 | 0.1594 | 0.1703 |
+| 16 | 5 | 0.1649 | 0.1603 | 0.1678 |
+| 17 | 5 | 0.1650 | 0.1607 | 0.1654 |
+| 18 | 5 | 0.1634 | 0.1604 | 0.1645 |
+| 19 | 5 | 0.1654 | 0.1641 | 0.1692 |
+| 20 | 5 | 0.1653 | 0.1611 | 0.1672 |
+
+Ceiling: 2.49x at 20 threads over 1, best 2.54x at 8; 2026-09-26: about 2.2x, flat from about 7.
+
+Quality: worst angle 0.3955 deg, median 45.00, share under 1 deg 0.00003, max degree 18, within tolerance True, Delaunay 0 violations of 641791 edges (80136 decided exactly), mesh sha256 `ccebf96a86c6c5e244e4a0281919de4e866fcfe789b66024a290ac2af33771a1`.
+
+## `tile`
+
+| threads | n | median s | min s | max s |
+|---:|---:|---:|---:|---:|
+| 0 | 5 | 0.1917 | 0.1833 | 0.1924 |
+| 1 | 5 | 0.4731 | 0.4649 | 0.5000 |
+| 2 | 5 | 0.3069 | 0.3007 | 0.3165 |
+| 3 | 5 | 0.2551 | 0.2515 | 0.2794 |
+| 4 | 5 | 0.2321 | 0.2262 | 0.2779 |
+| 5 | 5 | 0.2116 | 0.2104 | 0.2126 |
+| 6 | 5 | 0.2053 | 0.1993 | 0.2127 |
+| 7 | 5 | 0.1918 | 0.1911 | 0.2056 |
+| 8 | 5 | 0.1913 | 0.1842 | 0.2008 |
+| 9 | 5 | 0.1842 | 0.1810 | 0.1870 |
+| 10 | 5 | 0.1864 | 0.1834 | 0.1905 |
+| 11 | 5 | 0.1919 | 0.1817 | 0.1947 |
+| 12 | 5 | 0.1849 | 0.1825 | 0.1908 |
+| 13 | 5 | 0.1862 | 0.1823 | 0.2137 |
+| 14 | 5 | 0.1863 | 0.1826 | 0.1925 |
+| 15 | 5 | 0.1863 | 0.1854 | 0.1923 |
+| 16 | 5 | 0.1850 | 0.1845 | 0.1908 |
+| 17 | 5 | 0.1881 | 0.1825 | 0.1998 |
+| 18 | 5 | 0.1873 | 0.1823 | 0.1945 |
+| 19 | 5 | 0.1867 | 0.1823 | 0.1900 |
+| 20 | 5 | 0.1857 | 0.1844 | 0.1913 |
+
+Ceiling: 2.55x at 20 threads over 1, best 2.57x at 9; 2026-09-26: about 2.2x, flat from about 7.
+
+Quality: worst angle 0.6296 deg, median 45.00, share under 1 deg 0.00000, max degree 74, within tolerance True, Delaunay 0 violations of 692056 edges (94469 decided exactly), mesh sha256 `11741a81adfa17b34a1ea56875ec9e791414d17625eb1c58305e7c0cf10e005c`.
+
+<!-- bench.py: generated above this line; hand-written prose below is kept -->