Increment 36: GeoPackage output (--out X.gpkg) over sqlite3, QGIS-ready

ROADMAP.md

@@ -56,7 +56,7 @@ increment that most needs a picture to check against
 | 30d | Growing the catchment outline, made faster, with every mesh byte-identical: the same mitred buffer of the domain computed once instead of twice on the reprojected path; the CORINE search region grown from the domain's convex hull instead of taking the hull of the grown domain; and a long staircase outline (SMHI's raster-traced ones) grown in pieces of 1,000 edges and united, behind a gate that keeps smooth outlines on GEOS's single buffer. Measured by `@perf` on Lagan with the European CORINE GeoPackage: 103.4 s on master, about 17 s with all three (`docs/benchmarks/2026-10-07/30d-buffer/README.md`); from the performance review's finding F1 (PR #208). Ola, 2026-10-07: grow the outline in pieces; the gate is an open question, default yes | **merged (#215, 2026-10-07)** (branch `worktree-buffer-speed-2`; 29 net production lines; Python only, no refine or mesh code; Lagan 103.38 s → 16.97 s, meshes byte-identical; push waits on Ola) | `docs/increments/30d-outline-buffer-speed.md` |
 | 31 | Binary output by default: `rasputin mesh` writes binary `.vtk` (and `.ply`) unless `--ascii` is given, because the text writer formats one number at a time and costs minutes at basin scale (the performance review's finding F6 and question 3, PR #208). Ola, 2026-10-07: "change 3 to default binary". The flag pair `--binary/--ascii` stays; the help says binary is the default; README updated. Python only; no `@perf` run (not refine or mesh code) | **merged (#212, 2026-10-07)** (designed 2026-10-07; red `e3fe391a`/`9008811c`, green `29bc5ad3`; 4 net production lines against about 1, the overrun is `ruff format` wrapping the option); Ola ruled `.ply` binary too, as built; next the push | `docs/increments/31-binary-default.md` |
 | 32 | Land-cover borders simplified within a band, each class keeping its area exactly (Ola, 2026-10-08: "we need to implement is a CORINE simplification, similar to what we do on the outline of the auto-catchments"): the land-cover coverage cut into borders between junctions, each border reduced by the area-preserving segment collapse of the catchment outline (`reduce_ring`), with junctions and the domain outline fixed, a crossing test and a side test across all borders, and a band no border moves past (two-sided, proved by anchors on the source border). `--features-tolerance` becomes that band, default 50 m (half CORINE's positional accuracy), replacing GEOS's coverage simplifier; start angle 15° when land cover is simplified (band above 0), 25° otherwise. Prototype on the German fused catchment: 5.7 times fewer triangles at 50 m vertical tolerance | **merged (#222, 2026-10-11)** (green `82c0a7d4`, 519 net production lines, one PR; mutation round `30f939b6`, one mutant recorded as not killed); Ola ruled all four questions (2026-10-08); code review round 1 asked for changes (comments, docstrings and `project_structure.md` the change made false), fixed by `4d4f0411` and the docs commit after it; code review rounds 1-3 done (round 3 approved); `@perf`'s short timing check found the runs about 70 % slower (the outline rule on borders clipped to the catchment) and the worst angle far lower (the simplifier placing points millimetres from junctions); fix built (section 15: outline rule first, a 1 m clearance in the simplifier; green `05c8f73a`, 573 counted lines in all; mutation round M7 to M13 killed); its quick check `b192a2c6`: Numedalslagen −11 %, Lagan −1.3 %, outline rule and moved area as before, triangles −42 %; questions 5 and 6 decided on default while Ola was away, reversible; fix code review round 1 asked for prose fixes and the test file's red-step staging removed; fix code review round 3 approved; next the push on Ola's yes | `docs/increments/32-landcover-simplify.md` |
-| — | GeoPackage output, so QGIS opens a mesh as it opens any vector layer: `rasputin mesh --out X.gpkg`, chosen by the suffix of `--out` as `.vtk` and `.ply` are (Ola, 2026-10-05: "Is the vtk_to_phkg.py the way to go, or as a direct output from rasputin? I think the latter, and it should go into the roadmap and eventually into the codebase."). Proposed content: a `triangles` layer, one PolygonZ per triangle in the run's CRS (z at each of the three corners, as the stopgap writes), with the attributes the run already knows: land-cover code and class name (the code system's names, as `rasputin palette` has them), mean z, plan area; a layer of the constraint polygons (the `--features` polygons clipped to the domain, with their code); optionally a QGIS style for land cover stored in the file (the `layer_styles` table), so the triangles open coloured. Land-cover colour: when the mesh carries CORINE classes, each triangle and each constraint polygon also gets its class's colour as a `#rrggbb` text column, so a reader can colour by the column without the stored style (in QGIS, a data-defined fill colour). The Graz file of 2026-10-05 carries `clc_color` on its `triangles` and `corine` layers; its stored styles categorise on `clc_code` and do not read the column. The colours are those of `rasputin palette` (`src_python/tin_engine/palettes.py`), the same 26 the Graz file carries: natural colours, as Ola chose over the official CLC legend in 16c's ruling 2 (that legend paints peat bogs blue and the sea almost white). **One attribute is not in the mesh yet**: the ID of the source polygon each triangle lies in (a sub-catchment cell, a parcel). Today only a class code labels triangles, one code system per mesh (16c, 16e's D1); the stopgap finds the cell by a centroid lookup outside rasputin. It needs a per-triangle array from a polygon source's ID property: its own small step, done just before this writer, because the `.vtk` and `.ply` outputs gain the attribute too (Ola, 2026-10-05, accepting the default: "defaults on the cell ID"). Written with the standard library's `sqlite3` beside the existing reader `src_python/tin_engine/io/geopackage.py`, which can read it back in the tests; no GDAL, OGR or Fiona (CLAUDE.md §2). Python only: the C++ core is unchanged and never sees the file (the I/O boundary). Prior art and working spec: Ola's stopgap `vtk_to_gpkg.py`, outside the repository (`rasputin_data/lands/scripts/`), which converts a mesh `.vtk` with `sqlite3`, NumPy, Shapely and PyProj, no GDAL, and its first use, the Mur above Graz (2 451 734 triangles, 10 m DEM, tolerance 10 m, 2026-10-05): **the GeoPackage was 677.9 MB against 110.6 MB for the binary `.vtk`**; the whole stopgap run took 39 s and 4.1 GB peak, and that includes the centroid cell lookup and the CORINE clipping, so it is not the cost of the GeoPackage encoding alone: the design measures the writer on its own. A GeoPackage row holds the triangle's whole ring (four corners, the first repeated) and a bounding box, and the file adds a spatial index, where `.vtk` stores each vertex once. Questions for the design: file size (keep the spatial index or not; at most an opt-in variant of 2D polygons with mean z as an attribute, proposed default: full PolygonZ; writing in chunks to bound memory); the CRS entry for a CRS with no EPSG code; where the run record (increment 25) goes; whether the constraint lines get a layer of their own | proposed by Ola 2026-10-05 ("yes, roadmap row now, build after PR D"); not designed. Built after the Python audit's PR D (branch `audit-encoders`: the file encoders move out of `cli.py` into `io/`, and the mesh writers get one shared check; `docs/increments/python-audit.md` on branch `worktree-python-audit`, `8241e14`), so the new writer joins that shared check instead of a third copy | - |
+| 36 | GeoPackage output, so QGIS opens a mesh as it opens any vector layer: `rasputin mesh --out X.gpkg`, chosen by the suffix of `--out` as `.vtk` and `.ply` are (Ola, 2026-10-05: "Is the vtk_to_phkg.py the way to go, or as a direct output from rasputin? I think the latter, and it should go into the roadmap and eventually into the codebase."). Proposed content: a `triangles` layer, one PolygonZ per triangle in the run's CRS (z at each of the three corners, as the stopgap writes), with the attributes the run already knows: land-cover code and class name (the code system's names, as `rasputin palette` has them), mean z, plan area; a layer of the constraint polygons (the `--features` polygons clipped to the domain, with their code); optionally a QGIS style for land cover stored in the file (the `layer_styles` table), so the triangles open coloured. Land-cover colour: when the mesh carries CORINE classes, each triangle and each constraint polygon also gets its class's colour as a `#rrggbb` text column, so a reader can colour by the column without the stored style (in QGIS, a data-defined fill colour). The Graz file of 2026-10-05 carries `clc_color` on its `triangles` and `corine` layers; its stored styles categorise on `clc_code` and do not read the column. The colours are those of `rasputin palette` (`src_python/tin_engine/palettes.py`), the same 26 the Graz file carries: natural colours, as Ola chose over the official CLC legend in 16c's ruling 2 (that legend paints peat bogs blue and the sea almost white). **One attribute is not in the mesh yet**: the ID of the source polygon each triangle lies in (a sub-catchment cell, a parcel). Today only a class code labels triangles, one code system per mesh (16c, 16e's D1); the stopgap finds the cell by a centroid lookup outside rasputin. It needs a per-triangle array from a polygon source's ID property: its own small step, done just before this writer, because the `.vtk` and `.ply` outputs gain the attribute too (Ola, 2026-10-05, accepting the default: "defaults on the cell ID"). Written with the standard library's `sqlite3` beside the existing reader `src_python/tin_engine/io/geopackage.py`, which can read it back in the tests; no GDAL, OGR or Fiona (CLAUDE.md §2). Python only: the C++ core is unchanged and never sees the file (the I/O boundary). Prior art and working spec: Ola's stopgap `vtk_to_gpkg.py`, outside the repository (`rasputin_data/lands/scripts/`), which converts a mesh `.vtk` with `sqlite3`, NumPy, Shapely and PyProj, no GDAL, and its first use, the Mur above Graz (2 451 734 triangles, 10 m DEM, tolerance 10 m, 2026-10-05): **the GeoPackage was 677.9 MB against 110.6 MB for the binary `.vtk`**; the whole stopgap run took 39 s and 4.1 GB peak, and that includes the centroid cell lookup and the CORINE clipping, so it is not the cost of the GeoPackage encoding alone: the design measures the writer on its own. A GeoPackage row holds the triangle's whole ring (four corners, the first repeated) and a bounding box, and the file adds a spatial index, where `.vtk` stores each vertex once. Questions for the design: file size (keep the spatial index or not; at most an opt-in variant of 2D polygons with mean z as an attribute, proposed default: full PolygonZ; writing in chunks to bound memory); the CRS entry for a CRS with no EPSG code; where the run record (increment 25) goes; whether the constraint lines get a layer of their own | **built** (branch `worktree-gpkg-out`: red `87f41c3d`, green `c59246e5`, fix-up `7085dddf`, 377 net production lines, mutation round `dddce8e8` ten of ten killed; `@perf`'s acceptance `57c25074` ACCEPTED: Numedalslagen 3.93 s / 203 MB, Romsdalen 19.4 s / 810 MB, every mesh byte-identical; code review round 1's one blocker, an empty land-cover layer crashing the write, fixed by `55cc3c3e` + `18442b27`, 378 net lines; code review round 2 APPROVED; next the push on Ola's yes). Designed by `@architect`, 2026-10-11, on `81a551b3`; proposed by Ola 2026-10-05, "yes, roadmap row now, build after PR D"; PR D merged as #205, so the writer joins the encoders under `io/`). Measured on this design's probe: the quick check's Numedalslagen mesh (653,039 triangles) written in 4.0 s to 203 MB against 0.28 s and 29 MB for the binary `.vtk`, the stopgap 5.1 s; the Romsdalen slope mesh (3.1 million) in 21.5 s to 922 MB at 2.1 GB peak, the stopgap 28.7 s at 4.1 GB; the R-tree index is 60 % of the time; the same run writes the same bytes. Layers `triangles` (PolygonZ with code, name, colour, z_min, z_max, z_mean, area_m2, slope_deg), `lines`, `land_cover`, `land_cover_classes`, a `run` table with the run record, and the stored QGIS style; about 305 net lines, one PR; questions 1-5 answered with their defaults while Ola is away | `docs/increments/36-geopackage-output.md` |
 | — | A missing station number refused in plain words. In a reference-polygon file (`rasputin station-catchments --reference`), a feature whose `station` is null becomes the text "None" (`src_python/tin_engine/io/station_set.py@ed12512:89`). No station can match it, since a station number must look like `2.11.0` (`src_python/tin_engine/hydrography.py@ed12512:40`), so the polygon is dropped without a word and the station it belonged to is written with no comparison (`src_python/tin_engine/catchment_batch.py@ed12512:220`); two such polygons are refused as "station None appears twice". In a stations file a null number is refused, but in Pydantic's validation text rather than a sentence naming the file (`src_python/tin_engine/io/station_set.py@ed12512:79`). NVE's files have a number on every station, so only hand-made or third-party files are hit. The fix: both readers refuse a null `station`, naming the file; one or two lines. Plainer GeoJSON error messages in general are not planned until other people feed rasputin their own files | approved by Ola 2026-10-06; not designed (a small fix) | - |
 | — | A GeoTIFF whose coordinate-system GeoKeys hold garbage refused in plain words. The reader turns each code key into a whole number with `int()`, and a key stored as a double or as text can hold NaN, infinity or letters. Probed 2026-10-06 on hand-made files: NaN or text in `ProjectedCSTypeGeoKey` (3072) stops `rasputin mesh --dem` with Python's conversion message (for NaN, "Invalid value for --dem: cannot convert float NaN to integer"), which names neither the key nor the file, and infinity ends in a Python traceback (`OverflowError`) (`src_python/tin_engine/io/geotiff.py@ed12512:171`). The same holds for every other code key the reader reads: 1025, 4099, 3076, 3074, 3075, 2048 and 2054 (`src_python/tin_engine/io/geotiff.py@ed12512:310`, `src_python/tin_engine/io/geotiff.py@ed12512:176`, `src_python/tin_engine/io/geotiff.py@ed12512:339`, `src_python/tin_engine/io/geotiff.py@ed12512:352`, `src_python/tin_engine/io/geotiff.py@ed12512:362`, `src_python/tin_engine/io/geotiff.py@ed12512:499`, `src_python/tin_engine/io/geotiff.py@ed12512:423`). Worse, a fractional value is cut to a whole number first: 2.5 in 4099 is refused as 2, and 32633.5 in 3072 is read as EPSG:32633 without a word (`src_python/tin_engine/io/geotiff.py@ed12512:321`); 1.5 in 1025, 9001.7 in 4099 and 9001.5 in 3076 are likewise accepted as 1 and 9001 (probed by code review round 1 on this row's branch, exit 0). The parameter keys are already checked (`_one_number`, `src_python/tin_engine/io/geotiff.py@ed12512:454`). The fix: one check for every code key, a whole number or a refusal naming the key and its value. No real file is known to do this | approved by Ola 2026-10-06; not designed (a small fix) | - |
 | — | tifffile's warning about the NoData tag silenced. On a float32 DEM whose `GDAL_NODATA` (42113) is the lowest float32 written in 17 digits (`-3.4028234663852886e+38`), tifffile prints "parsing GDAL_NODATA tag raised ValueError: ... is not castable to float32" to stderr when the file is opened (`src_python/tin_engine/io/geotiff.py@ed12512:146`). It is noise: rasputin reads the tag itself and accepts the value, since it is that float32 exactly (`src_python/tin_engine/io/geotiff.py@ed12512:535`). The fix: silence that one tifffile message while the reader opens the file, leaving every other tifffile warning visible | approved by Ola 2026-10-06; not designed (a small PR) | - |