diff --git a/ROADMAP.md b/ROADMAP.md index 12348301..430874f5 100644 --- a/ROADMAP.md +++ b/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 | diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover.md b/docs/benchmarks/2026-09-29/bygdin-landcover.md new file mode 100644 index 00000000..96bfe765 --- /dev/null +++ b/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. diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/README.md b/docs/benchmarks/2026-09-29/bygdin-landcover/README.md new file mode 100644 index 00000000..71081328 --- /dev/null +++ b/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= bash docs/benchmarks/2026-09-29/bygdin-landcover/run.sh`, +then `python .../analyse.py ` and +`python .../render.py /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. diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/analyse.py b/docs/benchmarks/2026-09-29/bygdin-landcover/analyse.py new file mode 100644 index 00000000..35a03ae1 --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/analyse.py @@ -0,0 +1,212 @@ +"""Increment 16c acceptance on Bygdin (@perf, 2026-09-29): every table in +README.md comes from this script. Usage, from the repository root, .venv active: + + python docs/benchmarks/2026-09-29/bygdin-landcover/analyse.py + + holds run.sh's meshes. Prints to stdout; run.sh's logs are read from +./logs. The CORINE reference is read here straight from the GeoPackage with +sqlite3 and shapely, not through tin_engine.feature_input, so it is independent +of the code under test. +""" + +from __future__ import annotations + +import re +import sqlite3 +import statistics +import sys +from pathlib import Path + +import numpy as np +import shapely +from vtkmodules.util.numpy_support import vtk_to_numpy +from vtkmodules.vtkIOLegacy import vtkPolyDataReader + +HERE = Path(__file__).resolve().parent +ROOT = HERE.parents[3] +sys.path.insert(0, str(ROOT / "tools")) +import bench # noqa: E402 + +from tin_engine.landcover import regions # noqa: E402 + +GPKG = ROOT.parent / "rasputin_data" / "corine2018_dtm10_utm33.gpkg" +MARGIN = 2 * 1e-3 # cli: margin = 2 * snap spacing, default 1 mm +# 22's table (bygdin/README.md on increment22-autocatchment), percent. +TABLE_22 = {333: 42.34, 332: 21.32, 322: 17.00, 512: 16.46, 335: 2.44, 412: 0.33, 142: 0.10} + + +def gpkg_geometry(blob: bytes) -> shapely.Geometry: + """A GeoPackage binary geometry: 8-byte header, envelope, then WKB.""" + flags = blob[3] + envelope = {0: 0, 1: 32, 2: 48, 3: 48, 4: 64}[(flags >> 1) & 0b111] + return shapely.from_wkb(blob[8 + envelope :]) + + +def corine(domain: shapely.Geometry) -> list[tuple[shapely.Geometry, int]]: + """Every CORINE polygon whose envelope meets the domain's, unclipped.""" + x0, y0, x1, y1 = domain.bounds + con = sqlite3.connect(f"file:{GPKG}?mode=ro", uri=True) + rows = con.execute( + "select c.geom, c.code_18 from corine2018 c join rtree_corine2018_geom r on c.fid = r.id " + "where r.maxx >= ? and r.minx <= ? and r.maxy >= ? and r.miny <= ?", + (x0, x1, y0, y1), + ).fetchall() + out = [(gpkg_geometry(g), int(code)) for g, code in rows] + return [(g, c) for g, c in out if g.intersects(domain)] + + +def read(path: Path) -> dict: + r = vtkPolyDataReader() + r.SetFileName(str(path)) + r.ReadAllFieldsOn() + r.Update() + pd = r.GetOutput() + arr = pd.GetCellData().GetArray("land_cover_code") + n_lines, n_polys = pd.GetNumberOfLines(), pd.GetNumberOfPolys() + lines = vtk_to_numpy(pd.GetLines().GetConnectivityArray()).reshape(-1, 2) + tris = vtk_to_numpy(pd.GetPolys().GetConnectivityArray()).reshape(-1, 3) + fd = pd.GetFieldData().GetAbstractArray("land_cover_codes") + return dict( + cells=pd.GetNumberOfCells(), n_lines=n_lines, n_polys=n_polys, + codes=None if arr is None else vtk_to_numpy(arr).astype(np.int64), + points=vtk_to_numpy(pd.GetPoints().GetData()).astype(np.float64), + tris=tris.astype(np.int64), lines=lines.astype(np.int64), + text=None if fd is None else fd.GetValue(0), + scalars=pd.GetCellData().GetScalars().GetName() if pd.GetCellData().GetScalars() else None, + ) # fmt: skip + + +def rounds(tri: np.ndarray, edges: np.ndarray) -> tuple[int, np.ndarray]: + """landcover.regions' loop, counting its hook rounds (outer iterations + that hooked something). Checked equal to regions() by the caller.""" + n = int(max(tri.max(), edges.max())) + 1 + key = lambda p: np.minimum(p[:, 0], p[:, 1]) * n + np.maximum(p[:, 0], p[:, 1]) # noqa: E731 + sides = np.concatenate([tri[:, [0, 1]], tri[:, [1, 2]], tri[:, [2, 0]]]) + keys, owner = key(sides), np.tile(np.arange(len(tri)), 3) + free = ~np.isin(keys, key(edges)) + order = np.argsort(keys[free], kind="stable") + keys, owner = keys[free][order], owner[free][order] + pair = np.flatnonzero(keys[1:] == keys[:-1]) + u, v = owner[pair], owner[pair + 1] + parent, count = np.arange(len(tri)), 0 + while True: + ru, rv = parent[u], parent[v] + if np.array_equal(ru, rv): + return count, parent + count += 1 + np.minimum.at(parent, np.maximum(ru, rv), np.minimum(ru, rv)) + while not np.array_equal(parent, jumped := parent[parent]): + parent = jumped + + +def oracle(xy, tri, polys) -> tuple[np.ndarray, np.ndarray]: + """The design's oracle (R1): each triangle's centroid against the + polygons, smallest area wins, ties to the smaller code; and the mask of + triangles with inradius > 1.5 * margin, where it is exact.""" + a, b, c = xy[tri[:, 0]], xy[tri[:, 1]], xy[tri[:, 2]] + per = sum(np.hypot(*(q - p).T) for p, q in ((b, c), (c, a), (a, b))) + cross = (b[:, 0] - a[:, 0]) * (c[:, 1] - a[:, 1]) - (c[:, 0] - a[:, 0]) * (b[:, 1] - a[:, 1]) + r = np.abs(cross) / per + tree = shapely.STRtree([p for p, _ in polys]) + found, which = tree.query(shapely.points((a + b + c) / 3), predicate="intersects") + area = np.array([p.area for p, _ in polys]) + code = np.array([k for _, k in polys]) + order = np.lexsort((code[which], area[which], found)) + found, which = found[order], which[order] + head = np.unique(found, return_index=True)[1] + out = np.zeros(len(tri), dtype=np.int64) + out[found[head]] = code[which[head]] + return out, r > 1.5 * MARGIN + + +def stats_phases(path: Path) -> dict[str, float]: + text = path.read_text() + out = {m[1]: float(m[2]) for m in re.finditer(r"^\| ([^|]+?) \| ([0-9.]+) \| [0-9.]+ % \|$", text, re.M)} + total = re.search(r"^\| \*?\*?total\*?\*? \| \*?\*?([0-9.]+)", text, re.M) + if total: + out["total"] = float(total[1]) + return out + + +def wall(path: Path) -> tuple[float, int]: + text = path.read_text() + return float(re.search(r"([0-9.]+) real", text)[1]), int(re.search(r"(\d+)\s+maximum resident", text)[1]) + + +def main(scratch: Path) -> None: + domain = shapely.from_geojson((HERE / "bygdin_reduced_t20.geojson").read_text()) + domain = shapely.union_all([g for g in shapely.get_parts(domain)]) if domain.geom_type == "GeometryCollection" else domain + polys = corine(domain) + ref: dict[int, float] = {} + for p, k in polys: + ref[k] = ref.get(k, 0.0) + shapely.intersection(p, domain).area + ref_total = sum(ref.values()) + print(f"CORINE polygons meeting the domain: {len(polys)}; clipped area {ref_total / 1e6:.6f} km2, domain {domain.area / 1e6:.6f} km2") + + for tol in ("10", "1"): + tag = f"lc_t{tol}" + print(f"\n== {tol} m ==") + m = read(scratch / f"{tag}_run1.vtk") + codes, tri, lines, pts = m["codes"], m["tris"], m["lines"], m["points"] + print(f"vtkPolyDataReader: {m['cells']} cells = {m['n_lines']} lines + {m['n_polys']} triangles; active SCALARS {m['scalars']}") + assert codes is not None, "land_cover_code missing" + print(f"land_cover_code: {len(codes)} values (one per cell: {len(codes) == m['cells']}); " + f"nonzero on lines: {int((codes[: m['n_lines']] != 0).sum())}") # fmt: skip + print(f"FieldData land_cover_codes: {m['text']!r}") + tc = codes[m["n_lines"] :] + xy = pts[:, :2] + a, b, c = xy[tri[:, 0]], xy[tri[:, 1]], xy[tri[:, 2]] + area = 0.5 * np.abs((b[:, 0] - a[:, 0]) * (c[:, 1] - a[:, 1]) - (c[:, 0] - a[:, 0]) * (b[:, 1] - a[:, 1])) + total = area.sum() + print(f"mesh area {total / 1e6:.6f} km2; triangles with code 0: {int((tc == 0).sum())}") + print("| code | mesh km2 | mesh share % | CORINE clipped share % | diff (pp) | 22's table % |") + worst = 0.0 + for k in sorted(set(ref) | set(np.unique(tc).tolist()), key=lambda k: -ref.get(k, 0)): + ms = 100 * area[tc == k].sum() / total + rs = 100 * ref.get(k, 0.0) / ref_total + worst = max(worst, abs(ms - rs)) + print(f"| {k} | {area[tc == k].sum() / 1e6:.6f} | {ms:.5f} | {rs:.5f} | {ms - rs:+.6f} | {TABLE_22.get(k, '-')} |") + print(f"largest share difference: {worst:.6f} pp; codes equal: {set(np.unique(tc).tolist()) == set(ref)}") + # I1: across every interior non-constraint edge, the codes agree. + n = len(pts) + key = lambda p: np.minimum(p[:, 0], p[:, 1]) * n + np.maximum(p[:, 0], p[:, 1]) # noqa: E731 + sides = np.concatenate([tri[:, [0, 1]], tri[:, [1, 2]], tri[:, [2, 0]]]) + k3, own = key(sides), np.tile(np.arange(len(tri)), 3) + free = ~np.isin(k3, key(lines)) + o = np.argsort(k3[free], kind="stable") + ks, ow = k3[free][o], own[free][o] + pr = np.flatnonzero(ks[1:] == ks[:-1]) + print(f"I1: {len(pr)} interior unconstrained edges, {int((tc[ow[pr]] != tc[ow[pr + 1]]).sum())} with differing codes") + # I2: the centroid oracle. + want, exact = oracle(xy, tri, polys) + print(f"I2: oracle checks {int(exact.sum())} of {len(tri)} triangles (r > {1.5 * MARGIN} m); " + f"disagreements {int((want[exact] != tc[exact]).sum())}; " + f"disagreements on the {int((~exact).sum())} unchecked: {int((want[~exact] != tc[~exact]).sum())}") # fmt: skip + nr, parent = rounds(tri, lines) + same = np.array_equal(parent, regions(tri, lines)) + print(f"union-find: {nr} hook rounds; {len(np.unique(parent))} components; equal to landcover.regions: {same}") + # The ASCII run: same codes, and bench.quality. + asc = read(scratch / f"{tag}_ascii.vtk") + print(f"ASCII run codes equal to binary run1: {np.array_equal(asc['codes'], codes)}") + vm = bench.read_vtk_ascii(scratch / f"{tag}_ascii.vtk") + err = (HERE / "logs" / f"{tag}_ascii.err").read_text() + max_err = float(re.search(r"achieved max error ([0-9.e+-]+)", err)[1]) + q = bench.quality(vm.points, vm.triangles, vm.edges, float(tol), max_err) + print(f"quality: worst angle {q.worst_angle:.5f} deg, max degree {q.max_degree}, " + f"max error {max_err} (within: {q.within_tolerance}), Delaunay checked/ambiguous/violations " + f"{q.delaunay_checked}/{q.delaunay_ambiguous}/{q.delaunay_violations}") # fmt: skip + # Times over the three binary runs. + runs = [stats_phases(HERE / "logs" / f"{tag}_run{i}.stats.md") for i in (1, 2, 3)] + walls = [wall(HERE / "logs" / f"{tag}_run{i}.err") for i in (1, 2, 3)] + for ph in ("land cover", "refine", "features clip", "total"): + vals = [r.get(ph, float("nan")) for r in runs] + print(f"time {ph}: median {statistics.median(vals):.3f} s (runs {', '.join(f'{v:.3f}' for v in vals)})") + print(f"wall: median {statistics.median(w for w, _ in walls):.2f} s (runs {', '.join(f'{w:.2f}' for w, _ in walls)}); " + f"max RSS median {statistics.median(r for _, r in walls) / 1e6:.0f} MB") # fmt: skip + for i in (1, 2, 3): + line = [ln for ln in (HERE / "logs" / f"{tag}_run{i}.err").read_text().splitlines() if ln.startswith("land cover")] + print(f"run{i} stderr: {line[0]}") + + +if __name__ == "__main__": + main(Path(sys.argv[1])) diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/bygdin_landcover_oblique.png b/docs/benchmarks/2026-09-29/bygdin-landcover/bygdin_landcover_oblique.png new file mode 100644 index 00000000..fcb7defc Binary files /dev/null and b/docs/benchmarks/2026-09-29/bygdin-landcover/bygdin_landcover_oblique.png differ diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/bygdin_landcover_top.png b/docs/benchmarks/2026-09-29/bygdin-landcover/bygdin_landcover_top.png new file mode 100644 index 00000000..27ce3fe5 Binary files /dev/null and b/docs/benchmarks/2026-09-29/bygdin-landcover/bygdin_landcover_top.png differ diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/bygdin_reduced_t20.geojson b/docs/benchmarks/2026-09-29/bygdin-landcover/bygdin_reduced_t20.geojson new file mode 100644 index 00000000..5a63f92e --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/bygdin_reduced_t20.geojson @@ -0,0 +1 @@ +{"type": "FeatureCollection", "crs": {"type": "name", "properties": {"name": "EPSG:25833"}}, "features": [{"type": "Feature", "properties": {"seed": [8.5425, 61.3512], "seed_crs": "EPSG:4326", "nodes": 3049095, "fine_vertices": 17812, "fine_area_m2": 304909550.0, "reduced_vertices": 740, "reduced_area_m2": 304909549.9999997, "outline_tolerance_m": 20.0, "windows": [[1438, 3024], [2236, 3646], [3520, 4813]]}, "geometry": {"type": "Polygon", "coordinates": [[[137640.0, 6832215.0], [137597.74193548388, 6832243.172043011], [137458.29391628047, 6832099.4653071035], [137272.3691336265, 6832052.957247976], [137184.90341039226, 6831980.659617026], [137201.01864562993, 6831903.119034275], [136903.35714285713, 6831628.214285715], [136951.62316462147, 6831459.250051314], [136817.0754716981, 6831303.773584906], [136830.00038878224, 6831200.698025623], [136684.96189964417, 6831184.071259593], [136435.6189184878, 6831273.199872388], [136395.8257108954, 6831332.053447732], [135865.0, 6831270.0], [135885.40185471406, 6831247.897990727], [135635.0, 6830980.0], [135645.86030125848, 6830903.977891191], [135592.49409990516, 6830828.017502136], [135632.2171388938, 6830610.685336256], [135633.68724762482, 6830376.88081057], [135589.39452549722, 6830304.652792333], [135964.71533574187, 6829950.284664258], [136513.89924577868, 6829517.547606546], [136530.10616212134, 6829406.420685099], [136610.09717929043, 6829289.689134087], [136572.75051124743, 6829053.128834356], [136676.59931295138, 6828870.89274336], [136835.26730887045, 6828688.502266878], [136877.84676395185, 6828693.759262805], [137088.33333333334, 6828465.0], [137057.04474659674, 6828423.815620678], [137186.22954293294, 6828252.867802737], [137245.515878252, 6828247.818994536], [137496.96684188143, 6827874.340167733], [137586.42857142858, 6827642.857142857], [137591.7123287671, 6827533.287671233], [137650.25901788575, 6827515.823043272], [137722.1739130435, 6827398.043478261], [137818.72049870715, 6827340.082686023], [138004.6911740526, 6827112.268315339], [138129.00346474705, 6827096.180873251], [138191.4623733689, 6826868.500055129], [138135.0, 6826840.0], [138257.3819876784, 6826602.968908627], [138238.62261836408, 6826547.302507902], [138328.53135983145, 6826454.936038019], [138303.061546746, 6826412.9448269205], [138419.98539481877, 6826358.102763013], [138385.0, 6826290.0], [138514.84023017494, 6826180.874040825], [138576.58487743896, 6825839.170865189], [138718.05190812348, 6825860.803012173], [138802.4414514038, 6825760.811100701], [138781.7282051282, 6825686.420512821], [138844.19485053458, 6825643.713724635], [138863.115942029, 6825468.115942029], [138808.125, 6825413.125], [138893.70368496655, 6825265.600430721], [138929.1959398952, 6825247.993253361], [138960.0, 6825125.0], [139175.0, 6824910.0], [139116.94111216627, 6824813.235186944], [138949.77487567315, 6824698.088308348], [139098.8443611255, 6824566.148784354], [139065.0, 6824470.0], [139120.51022533225, 6824410.39306786], [138917.7332545976, 6824214.83718576], [138941.75736597166, 6824169.45830872], [139011.99137029288, 6824182.391474895], [139117.5, 6824107.5], [139035.84347231753, 6823966.130352153], [139261.66666666666, 6823673.333333333], [139288.7288135593, 6823534.406779661], [139345.0, 6823450.0], [139371.095371097, 6823292.340466289], [139498.8679245283, 6823156.509433962], [139550.6470393466, 6823014.014679101], [139505.71212772815, 6822934.286085522], [139649.8381673627, 6822703.225457756], [139763.43880354054, 6822680.531590193], [139845.0, 6822563.895642195], [139845.0, 6822510.0], [139933.17778456188, 6822421.822215438], [140025.35194260805, 6822426.207276433], [140062.3109441052, 6822509.62188821], [140168.80157170922, 6822471.640471512], [140210.0, 6822525.0], [140505.9208167372, 6822229.079183263], [140604.20249812232, 6822247.359009317], [140804.2085916462, 6821961.682867177], [140899.38604496184, 6822003.772089924], [140936.32, 6822077.64], [141007.86755289164, 6822086.714784027], [141157.0, 6821912.0], [141098.07046615635, 6821872.320620127], [141142.49400479617, 6821822.482014389], [141112.51491139995, 6821701.06200086], [141355.0, 6821440.0], [141336.42857142858, 6821318.571428572], [141365.0, 6821290.0], [141224.0416141236, 6821147.389659521], [141408.40847831414, 6820937.04239157], [141395.0, 6820870.0], [141499.44444444444, 6820765.555555556], [141465.45526353584, 6820730.716645124], [141593.33333333334, 6820571.666666667], [141569.6974925791, 6820536.908077322], [141659.26065162907, 6820424.260651629], [141651.84448836243, 6820324.296281796], [141772.8260869565, 6820162.173913044], [141816.25, 6820061.25], [142002.6073224465, 6820199.636999488], [142075.57326495965, 6820154.970894979], [142145.0, 6820180.0], [142126.74311926606, 6820253.302752294], [142205.56356523783, 6820277.275755774], [142217.27272727274, 6820382.2727272725], [142326.66666666666, 6820308.333333333], [142388.09997058453, 6820353.3151237555], [142415.0, 6820220.0], [142562.42090701312, 6820072.579092987], [142496.78952869328, 6819997.666816089], [142545.0, 6819670.0], [142500.23076923078, 6819625.230769231], [142605.5162615763, 6819519.49576577], [142622.45828279847, 6819373.558505002], [142764.9480968858, 6819290.017301038], [142802.37950475726, 6819326.935949897], [142957.33098629283, 6819326.905978922], [143108.52788881605, 6819207.685083574], [143320.7669052955, 6819176.611933246], [143329.91124260356, 6819067.061143984], [143440.0, 6818955.0], [143500.0, 6818955.0], [143668.3153049157, 6818845.289064038], [143728.84615384616, 6818876.923076923], [143832.3291834707, 6818814.509030335], [143810.82940272454, 6818728.670749027], [143766.25954198473, 6818729.122137405], [143798.94779576827, 6818644.411090414], [143724.88005997002, 6818546.64167916], [143728.33333333334, 6818373.333333333], [143665.70455182076, 6818215.276290152], [143951.67372638808, 6818022.302804808], [144000.3285880899, 6817943.895927811], [144332.65235579392, 6817745.0], [144365.87257886506, 6817745.0], [144498.42865702798, 6817640.664382635], [144528.09024704012, 6817525.1893014405], [144490.5572646479, 6817377.004492921], [144252.18583446817, 6817108.268572732], [144020.63138396858, 6817064.226122432], [144103.0, 6816882.0], [144078.0923694779, 6816670.361445783], [144277.27753076606, 6816595.668464988], [144557.5, 6816682.5], [144680.52583922824, 6816643.4437025385], [144647.251725348, 6816502.276289624], [144749.5221573919, 6816395.708043541], [144745.86872586873, 6816358.281853282], [144850.0, 6816345.0], [145040.46826041545, 6816171.163978183], [145023.57142857142, 6816067.142857143], [145102.99888205703, 6815921.453325881], [145065.0, 6815880.0], [145013.74513079796, 6815591.887975379], [145075.07522179233, 6815561.809962153], [145124.17391542502, 6815188.856079254], [145175.52631578947, 6815144.736842105], [145195.0, 6815040.0], [145155.0, 6815000.0], [145172.14285714287, 6814947.142857143], [145111.0, 6814886.0], [145139.67534438922, 6814816.195210808], [145468.04373820877, 6814477.181769896], [145615.5698981443, 6814257.181863997], [145604.56123577632, 6814219.313332136], [145914.0, 6813891.8], [146001.37096774194, 6813903.629032258], [146144.44444444444, 6813760.555555556], [146350.0, 6813645.0], [146335.0, 6813620.0], [146401.43835616438, 6813606.712328767], [146405.0, 6813460.0], [146366.15659509707, 6813430.677432868], [146527.69600301265, 6813171.358123098], [146503.70486860673, 6813136.148679506], [146611.05376965224, 6813106.046539101], [146847.60120585765, 6813085.9007122805], [146966.527700098, 6812958.03581416], [147052.18302184745, 6812997.871228581], [147169.1780821918, 6812955.273972603], [147196.0, 6812995.0], [147346.35462052366, 6812939.574463671], [147421.92532942898, 6812993.97510981], [147479.09752825845, 6812965.4824469825], [147558.09108487266, 6813045.788113164], [147519.6041081912, 6813077.697945904], [147608.4982822893, 6813168.276441781], [147667.1334431631, 6813151.573311367], [147790.91385169327, 6813284.549521231], [147842.47737909516, 6813293.602184087], [148220.3207003201, 6813636.607356718], [148361.08580106302, 6813519.130599848], [148460.0, 6813485.0], [148493.6742456259, 6813531.714563435], [148606.3060144191, 6813516.527023941], [148728.69185386883, 6813553.855230316], [148879.0, 6813714.0], [148892.34732105123, 6813769.057699337], [149049.7259037944, 6813832.720556439], [149340.6027820711, 6813827.364760432], [149443.44871181578, 6813949.081502975], [149783.62941755238, 6814088.053182289], [149787.3673469388, 6814149.102040816], [149913.8162243507, 6814447.021466103], [150027.98398457805, 6814525.496019012], [150141.04797171257, 6814506.503744745], [150194.78260869565, 6814540.217391305], [150435.0, 6814530.0], [150360.75, 6814604.25], [150416.25, 6814698.75], [150445.8695652174, 6814670.869565218], [150597.77777777778, 6814822.777777778], [150799.40637330865, 6814738.635584262], [150994.22455816425, 6814787.81637211], [151012.97354325702, 6814963.590357457], [150951.63579225936, 6815083.9067151835], [150962.72727272726, 6815115.0], [151030.48751697718, 6815110.134990859], [151102.07909634485, 6815245.279159011], [151271.65451519756, 6815239.984295697], [151308.7682386734, 6815311.079564392], [151460.57106393375, 6815222.896297737], [151563.125, 6815341.875], [151684.30056524486, 6815324.256142117], [151834.57169459964, 6815543.575418995], [151893.27946552495, 6815707.847712303], [151880.50264550265, 6815858.783068783], [151990.6010042372, 6815986.497164915], [152060.12504789885, 6816209.749904202], [152035.0, 6816290.0], [152084.0909090909, 6816311.515151516], [152092.80805687205, 6816360.912322275], [152299.2829360342, 6816553.41516797], [152422.85714285713, 6816443.571428572], [152477.49710791605, 6816486.557593786], [152521.42857142858, 6816456.428571428], [152592.2393430768, 6816512.669524852], [152596.76479812863, 6816637.450847148], [152549.84938495586, 6816705.217555064], [152762.9001756771, 6816982.259353165], [152932.75378831074, 6816779.468768818], [153164.8737720111, 6816642.538832252], [153165.02191609266, 6816592.614276769], [153287.77621563748, 6816513.054280999], [153606.79801008, 6816526.356951625], [153739.00240476284, 6816450.161591488], [153688.694621263, 6816254.935961682], [153870.23303808243, 6816121.92212627], [153852.40340794573, 6816058.978151899], [154044.0, 6816041.0], [154143.77704498603, 6815941.222955014], [154418.8745345136, 6816050.63404243], [154483.23751478214, 6816023.997543892], [154682.17736722674, 6816036.966404154], [154750.1219986578, 6816065.8835090455], [154933.96544251748, 6816022.994869282], [155041.75461367719, 6816104.826975314], [155126.02943487625, 6816021.415705997], [155364.18056245425, 6815908.366837521], [155473.4518828452, 6815801.338912134], [155576.0787598512, 6815918.459955267], [155673.1800339771, 6815792.479119474], [155757.8260869565, 6815830.217391305], [155829.44544463628, 6815790.353100789], [155879.4713349174, 6815939.219827721], [156068.19036510444, 6815921.3572261715], [156182.09302325582, 6815835.930232558], [156202.26098923897, 6815726.969625716], [156260.625, 6815685.0], [156376.49083068728, 6815772.95243927], [156345.0, 6815850.0], [156415.0, 6815960.0], [156350.33945324935, 6816024.449926077], [156217.47260584292, 6816064.129555845], [156130.1851851852, 6816183.333333333], [156111.44192111958, 6816302.650209459], [156391.81677091605, 6816286.592340345], [156426.4556794017, 6816495.392917727], [156558.90990598194, 6816665.272419229], [156707.94926474956, 6816481.194127676], [157091.39304937274, 6816355.131812555], [157100.0, 6816425.0], [157305.47173560766, 6816362.353989709], [157390.1019995881, 6816381.039124564], [157273.75, 6816628.75], [157337.97297297296, 6816740.540540541], [157316.11111111112, 6816773.333333333], [157380.42857142858, 6816864.571428572], [157345.0, 6816900.0], [157474.2857142857, 6817029.285714285], [157510.0, 6816995.0], [157604.04255319148, 6817011.808510638], [157765.5383315261, 6816937.592302455], [157812.5, 6816973.0], [157925.68808132003, 6816903.948551952], [157895.0, 6816803.333333333], [157951.5116479476, 6816690.878120243], [157875.94978325407, 6816565.120443983], [157930.0, 6816505.0], [158094.7585828061, 6816526.831066811], [158441.32971239905, 6816305.462029215], [158484.7425422148, 6816089.556151916], [158755.6582352538, 6815868.595579216], [158717.7460571056, 6815769.164069157], [158766.58437423463, 6815731.4305657605], [159016.40058227236, 6815653.148119857], [158997.54284853733, 6815563.793145235], [159254.1424648113, 6815392.845778662], [159190.83333333334, 6815303.229166667], [159241.66666666666, 6815233.333333333], [159330.39325842698, 6815254.7191011235], [159444.41583315123, 6815179.457590648], [159518.08560806658, 6815001.26311421], [159610.0, 6814955.0], [159655.0, 6815000.0], [159669.50205278205, 6815115.172084419], [159871.07197186654, 6815149.299696093], [159952.22222222222, 6815107.222222222], [160069.10729420255, 6815179.692141028], [160134.23048735395, 6815139.763186853], [160204.5, 6815139.5], [160398.81895860942, 6815238.118430292], [160522.43506943982, 6815104.396783581], [160624.0604394343, 6815155.430899723], [160830.04368030006, 6814930.551773602], [161066.24030172775, 6814855.777698366], [161350.01393145724, 6814953.180551685], [161510.0, 6814915.0], [161598.26306397247, 6814975.643329668], [161780.19607843139, 6814952.647058823], [161840.0, 6814845.0], [161876.86315789475, 6814861.484210527], [161954.70612328488, 6814750.924636494], [162190.0934019556, 6814568.886670609], [162169.10988743644, 6814521.359502801], [162266.37049742692, 6814460.777893015], [162299.0909090909, 6814490.454545454], [162350.0, 6814465.0], [162403.04314329737, 6814538.998459168], [162378.6015243573, 6814574.734254772], [162522.26744186046, 6814651.220930233], [162516.74387723103, 6814616.8511016825], [162587.29036769632, 6814563.175514227], [162620.0, 6814415.0], [162751.02946044356, 6814479.826216484], [162840.0, 6814475.0], [162961.65952121164, 6814353.340478788], [163276.6738121217, 6814158.326187878], [163248.3261878783, 6814106.673812122], [163318.57142857142, 6814036.428571428], [163410.0, 6814085.0], [163484.88372093023, 6814020.813953488], [163609.2154664364, 6814065.748037321], [163892.62440494008, 6813782.779033474], [163741.62921348313, 6813622.471910113], [163552.71878099983, 6813603.240766249], [163471.222741433, 6813515.288161994], [163552.03058304757, 6813460.807038269], [163423.73463132296, 6813291.364642098], [163644.32819325296, 6813163.016771873], [163567.88869182774, 6812999.652704162], [163698.85066599218, 6812751.698001093], [163635.8378165902, 6812658.929392556], [163821.70311841829, 6812385.71047471], [163808.98071837446, 6812291.979156205], [163917.43243243243, 6812173.513513514], [163923.49914439756, 6812025.5785719445], [164084.99727816283, 6811926.5085076], [164173.68676790822, 6811790.963314203], [164263.90259560433, 6811485.964099594], [164428.94132231225, 6811306.58511767], [164441.15527007577, 6811201.8750435], [164543.2696481575, 6811027.505426047], [164525.81517928684, 6810915.243012267], [164687.76002290647, 6810745.419509138], [164648.83532978804, 6810678.5318045085], [164645.0, 6810520.0], [164790.76588878914, 6810427.424689603], [165041.0616470588, 6810353.613647059], [165255.01829893093, 6810409.439022539], [165355.3015424999, 6810379.725468889], [165807.77777777778, 6810492.777777778], [165862.34620659557, 6810453.302744165], [166164.88580153522, 6810652.314058491], [166357.33333333334, 6810733.222222222], [166504.79644088983, 6810702.478892937], [166738.2570838119, 6810924.066940101], [167062.79121971395, 6810825.960480761], [167499.8949211909, 6810833.021015761], [167622.72727272726, 6810912.2727272725], [167650.0, 6810885.0], [167833.04171334364, 6810991.066069858], [167898.69837386673, 6811139.245934667], [168105.2605042017, 6811334.478991597], [168024.4240491967, 6811467.551301144], [168059.70967741936, 6811663.548387097], [168036.46786122074, 6811744.73555327], [168048.9334216628, 6811997.866843326], [168105.70533241477, 6812299.294667585], [168046.04465184055, 6812358.955348159], [168134.59571527297, 6812469.595715273], [168120.12841631015, 6812599.02955548], [168305.32710050038, 6812698.043566311], [168440.4606585253, 6812940.626835678], [168455.73499161843, 6812932.546442499], [168521.4814814815, 6813123.518518519], [168216.6542235406, 6813511.0663074], [168345.10455040447, 6813676.958179838], [168406.41025641025, 6813652.435897436], [168454.5945945946, 6813719.594594595], [168498.42973970083, 6813714.518389612], [168563.73029364986, 6813866.044721538], [168520.9112981316, 6814040.718329472], [168533.679245283, 6814088.867924528], [168417.0148486697, 6814242.79852416], [168314.74741928652, 6814237.949877645], [168117.49437586972, 6814455.506079413], [168229.12993961514, 6814745.610354501], [168439.9903529381, 6814713.798655271], [168566.31719061558, 6814723.66341645], [168707.1066854023, 6814894.870021078], [168654.27544051432, 6814970.128122582], [168356.7192067656, 6815072.068766756], [168204.62177231206, 6815255.0], [168042.2355061998, 6815234.290297861], [167990.77338129497, 6815283.543165468], [167837.66177546134, 6815323.520864536], [167901.0237011057, 6815437.318948756], [167523.23344481905, 6815822.238529126], [167451.0630340981, 6815801.026141856], [167306.06762680024, 6815847.0788979335], [167190.0, 6815705.0], [167044.6100144439, 6815835.117958594], [167237.1842529126, 6816083.455123402], [167209.17775182932, 6816131.461973824], [167355.66460390794, 6816329.335396092], [167329.44444444444, 6816355.555555556], [167387.96205630354, 6816616.499388005], [167497.21773979685, 6816778.822526856], [167451.0272743894, 6816847.633378786], [167549.15078972108, 6817139.476878677], [167243.98866247002, 6817199.603574305], [167174.06875028307, 6817294.2011025585], [166996.3145335883, 6817359.908163894], [166971.0, 6817472.0], [166884.0909090909, 6817515.454545454], [166775.05378511976, 6817506.144659716], [166566.31215009288, 6817598.261576196], [166549.1625748903, 6817751.680748954], [166492.55927025742, 6817753.93929139], [166233.4210483845, 6818133.2116055], [166115.3551138529, 6818471.641239512], [166302.33333333334, 6818686.666666667], [166341.41249586048, 6818833.392177055], [166734.65753424657, 6819245.0], [166785.99877184516, 6819232.320050565], [167011.5880375926, 6819333.411962408], [166880.89496928419, 6819542.104949057], [166924.33247057255, 6819628.391934455], [167058.63821719208, 6819722.779751538], [167131.96713615023, 6819871.615023474], [167098.4764497793, 6819919.100137919], [167161.05042016806, 6820079.075630252], [167146.6634105652, 6820200.538060977], [167282.5, 6820387.5], [167248.75347801892, 6820436.769616026], [167256.65750346836, 6820618.65411646], [167206.68033090775, 6820686.246072778], [167267.62458181655, 6820786.998011628], [167215.5, 6820851.833333333], [167259.5918349431, 6820946.2756069675], [167216.78021978022, 6821043.472527472], [167354.9859947679, 6821224.385468366], [167282.85714285713, 6821420.238095238], [167235.0, 6821430.0], [167265.53796103233, 6821483.618292585], [166947.04765552984, 6821784.847523488], [166720.31987415714, 6821775.924246227], [166235.69674469688, 6822241.87931986], [165934.52915906822, 6822079.799817318], [165911.88442211057, 6821905.326633166], [165679.2638544093, 6822102.583799018], [165423.29185031544, 6822158.237403855], [165223.8708631178, 6822108.915472963], [165124.66666666666, 6822119.666666667], [165028.29268292684, 6822023.292682927], [164829.50572649855, 6822175.708886004], [164815.0, 6822247.133283693], [163998.30829096286, 6823129.58053719], [163792.85178620127, 6823026.627806563], [163805.71543846332, 6823001.367150615], [163419.6545146307, 6822556.426685418], [163341.53802613975, 6822599.542039724], [163180.67637670063, 6822569.993789102], [163078.72477212505, 6822603.932886762], [162872.4742136763, 6822509.121809682], [162653.07692307694, 6822468.076923077], [162592.44541484717, 6822515.305676856], [162393.33333333334, 6822501.666666667], [162248.53606027988, 6822609.230355221], [162144.63927978635, 6822637.4066103045], [161895.76822711964, 6822882.543396485], [161826.2174817898, 6823161.207075963], [161645.77714427325, 6823342.256978973], [161562.80898876404, 6823267.921348315], [161490.0, 6823285.0], [161266.79030192975, 6822970.667908133], [160971.81136296896, 6822883.205454997], [160959.26763371166, 6822774.73558198], [160806.80737806426, 6822591.108377964], [160412.43243243243, 6822545.810810811], [160232.1669561144, 6822688.572761424], [160162.72727272726, 6822677.337662337], [160116.66826003825, 6822705.004780115], [160145.0, 6822760.0], [159971.1759431634, 6822771.255124691], [159841.0655737705, 6822924.918032787], [159771.92534466315, 6822925.010838464], [159635.0, 6823080.0], [159650.36697247706, 6823161.758409786], [159531.98256241117, 6823235.420569002], [159420.0, 6823155.0], [159301.30434782608, 6823135.217391305], [159260.0, 6823205.0], [159084.67146380994, 6823195.956785547], [158756.73403076752, 6823385.953928209], [158718.37016574584, 6823475.46961326], [158601.9285083186, 6823496.844012202], [158573.4, 6823588.0], [158618.70986920333, 6823633.745541022], [158516.53957800812, 6823634.453374332], [158514.84711211777, 6823750.481313704], [158377.36809690003, 6823777.6319031], [158296.42857142858, 6823858.571428572], [158336.72951870546, 6824128.23415214], [158552.46797385416, 6824354.995921588], [158546.4730999146, 6824410.388556789], [158631.45851142038, 6824504.988622307], [158627.8440081232, 6824573.971657521], [158751.9642857143, 6824783.035714285], [158635.0, 6824900.0], [158610.74017743874, 6825026.858048174], [158528.06397840928, 6825124.598348944], [158584.43570265872, 6825242.722011213], [158517.10595515018, 6825372.10595515], [158679.34782608695, 6825534.347826087], [158626.98421610988, 6825650.749890815], [158247.49005137023, 6826057.537393148], [157917.51924390567, 6826093.847710612], [157786.67985941586, 6826008.639392729], [157378.7747430343, 6826299.240397112], [157280.0, 6826305.0], [156995.2896710623, 6826597.929300996], [156855.0, 6826887.580645162], [156855.0, 6826987.142857143], [156781.4229437834, 6827014.580683301], [156694.18284996282, 6826960.430362051], [156580.0, 6827035.0], [156427.22619874487, 6826978.528233891], [156322.20362574054, 6826978.105847667], [155696.28811047922, 6826444.645677502], [155698.97111913358, 6826325.956678701], [155648.69638047065, 6826220.796272789], [155678.87947168696, 6826174.505945391], [155347.17741730215, 6825877.605460952], [155324.0, 6825899.0], [155239.92119495236, 6825793.362526991], [155329.7510011255, 6825692.161659015], [155235.0, 6825520.0], [155273.1181127542, 6825469.420196538], [155244.836861259, 6825398.880280872], [155275.0, 6825185.998955788], [155144.2457153856, 6825030.026256759], [155146.80368488174, 6824888.237334569], [154966.50743680075, 6824752.349697339], [154940.0, 6824675.0], [154672.96470588236, 6824566.05882353], [154562.5179401709, 6824615.016627856], [154475.2427402175, 6824558.400979183], [154061.17947062736, 6824952.682466841], [153948.3932853717, 6824889.604316547], [153910.0, 6824925.0], [153633.8414655395, 6824675.296095573], [153471.88416499528, 6824573.500874762], [153426.8792042776, 6824646.075720117], [153309.49074264755, 6824577.633097253], [153247.0679408379, 6824447.49916059], [153120.6811731315, 6824334.545884579], [152990.0, 6824405.0], [152808.32406471862, 6824259.991821811], [152677.78021868292, 6824315.803335164], [152626.22609219502, 6824579.0548049705], [152384.7992433376, 6824779.837076458], [152278.09075573212, 6824846.237759655], [152213.0, 6824712.0], [152161.25, 6824706.25], [152091.3911917246, 6824766.539108512], [151782.18418519443, 6824750.416064381], [151575.3332200404, 6824670.751515386], [151482.09492582743, 6824535.961711694], [151237.8509637461, 6824588.767535829], [151048.6855839793, 6824570.609580332], [150772.67070690743, 6824635.934637659], [150647.1180714055, 6824694.774777073], [150517.2738001314, 6824882.720017374], [150490.57487511373, 6824992.36659076], [150368.1541865854, 6825058.422906707], [150143.10841580582, 6825308.745261659], [150059.24676187313, 6825284.784789107], [149914.90366404055, 6825339.163638309], [149866.64025356577, 6825380.546751189], [149803.0982650025, 6825359.3660883345], [149713.27636316314, 6825441.684162582], [149593.7657585323, 6825331.10442416], [149547.14285714287, 6825340.714285715], [149490.0, 6825295.0], [149274.90679094542, 6825546.364846871], [149305.0, 6825590.0], [149096.13636363635, 6825798.863636363], [149137.1672067933, 6825876.02289043], [148925.0, 6826050.0], [148928.31606217616, 6826129.481865285], [148692.09151447949, 6826397.939703955], [148622.70970824885, 6826338.784241047], [148560.55555555556, 6826383.333333333], [148590.88585017837, 6826421.236623067], [148546.1879647242, 6826482.32734796], [148627.16378508764, 6826627.871523432], [148499.26859107963, 6826755.395468691], [148682.25090291057, 6826952.749097089], [148666.00580270792, 6827060.038684719], [148705.0, 6827133.333333333], [148653.18181818182, 6827111.818181818], [148523.04537521815, 6827177.068062827], [148547.29508196723, 6827213.442622951], [148456.21325244984, 6827278.492767149], [148465.0, 6827340.0], [148325.0, 6827520.0], [148315.0, 6827586.0], [148135.7594936709, 6827449.24050633], [148000.12285927028, 6827414.87714073], [147905.08593415626, 6827465.951652011], [147724.3533032997, 6827485.0], [147546.80814456358, 6827709.264580327], [147618.99196179604, 6827783.092181444], [147447.1647981899, 6827912.16479819], [147353.0386189563, 6827864.118452933], [147237.3367831924, 6827982.47275187], [147137.62411810685, 6828023.250587929], [147043.8427676554, 6827940.959470064], [146646.2475308213, 6828237.569703922], [146510.6303528427, 6828282.464883227], [146480.02232370205, 6828326.315199159], [146243.46153846153, 6828433.269230769], [146166.36363636365, 6828385.909090909], [145892.8579279267, 6828449.293074313], [145698.8515359408, 6828273.418585753], [145603.37953091683, 6828313.336886994], [145625.0, 6828380.0], [145239.0185801349, 6828819.417609053], [145283.04689674792, 6828925.031260869], [145148.23529411765, 6829058.529411765], [145163.20202680613, 6829149.718862373], [144653.03129221033, 6829692.410361127], [144531.66666666666, 6829883.333333333], [144568.93858913428, 6829988.851458391], [144377.39666944387, 6830062.1077173175], [143979.65746991846, 6830084.060252395], [143870.0, 6830045.0], [143733.7129986326, 6830163.906018191], [143487.2764752713, 6830596.022720626], [143333.0561105112, 6830726.501611498], [143216.97162517396, 6830920.114947552], [142763.93285686092, 6831292.534366456], [142450.3096744802, 6831305.16582353], [142317.48989381708, 6831369.199303599], [142189.6591747816, 6831304.743607308], [142082.6455026455, 6831372.883597883], [142014.54293797343, 6831352.083016457], [141930.4563679374, 6831434.037702831], [141864.5673076923, 6831389.567307692], [141799.14123013322, 6831405.247070234], [141787.48680984325, 6831323.967080333], [141636.66666666666, 6831181.666666667], [141475.60008149847, 6831197.302349077], [141233.0548628429, 6830976.832917706], [141509.8184635631, 6830604.366053854], [141495.33898305084, 6830560.677966102], [141564.46000643945, 6830509.770019778], [141545.77362714784, 6830458.905365103], [141632.36842105264, 6830352.631578947], [141545.0, 6830250.0], [141623.8377763633, 6830200.622274657], [141686.23654181996, 6829929.228708523], [141649.55835962144, 6829859.681913775], [141715.0, 6829800.0], [141666.44749290444, 6829677.47398297], [141721.66666666666, 6829581.111111111], [141709.79804479572, 6829513.855587176], [141562.36967840735, 6829412.856355283], [141454.0950811691, 6829492.618725208], [141175.23841093152, 6829271.356306608], [141034.58521433573, 6829410.063663613], [140909.11902109734, 6829406.78201276], [140730.0, 6829515.0], [140600.0, 6829435.0], [140618.06886056298, 6829398.782824378], [140532.1081036239, 6829317.67283945], [140470.0, 6829355.0], [140373.125, 6829268.194444444], [140226.71247727206, 6829229.2699248735], [140045.26763457104, 6829090.899232198], [139907.9834137516, 6829228.536762], [139848.93116912674, 6829387.837838343], [139537.18753273966, 6829777.593713987], [139457.36522433392, 6830009.841604134], [139471.75458260515, 6830190.407476903], [139361.69032642775, 6830442.3604914965], [139385.76923076922, 6830663.846153846], [139297.1158099488, 6830979.101254613], [138957.26103517492, 6831621.836906897], [138742.14602379626, 6831882.475498062], [138744.0186756718, 6831907.097402352], [138398.62650113303, 6832267.215780782], [138303.2242100864, 6832205.139652085], [138029.34806321948, 6832279.209163358], [137863.5294117647, 6832232.647058823], [137825.43219565472, 6832281.692015447], [137640.0, 6832215.0]]]}}]} \ No newline at end of file diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/corine_natural.json b/docs/benchmarks/2026-09-29/bygdin-landcover/corine_natural.json new file mode 100644 index 00000000..38788392 --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/corine_natural.json @@ -0,0 +1,239 @@ +[ + { + "Name": "rasputin CORINE natural", + "Annotations": [ + "0", + "0 No polygon; constraint lines", + "111", + "111 Continuous urban fabric", + "112", + "112 Discontinuous urban fabric", + "121", + "121 Industrial or commercial units", + "122", + "122 Road and rail networks and associated land", + "123", + "123 Port areas", + "124", + "124 Airports", + "131", + "131 Mineral extraction sites", + "132", + "132 Dump sites", + "133", + "133 Construction sites", + "141", + "141 Green urban areas", + "142", + "142 Sport and leisure facilities", + "211", + "211 Non-irrigated arable land", + "212", + "212 Permanently irrigated land", + "213", + "213 Rice fields", + "221", + "221 Vineyards", + "222", + "222 Fruit trees and berry plantations", + "223", + "223 Olive groves", + "231", + "231 Pastures", + "241", + "241 Annual crops associated with permanent crops", + "242", + "242 Complex cultivation patterns", + "243", + "243 Land principally occupied by agriculture, with significant areas of natural vegetation", + "244", + "244 Agro-forestry areas", + "311", + "311 Broad-leaved forest", + "312", + "312 Coniferous forest", + "313", + "313 Mixed forest", + "321", + "321 Natural grasslands", + "322", + "322 Moors and heathland", + "323", + "323 Sclerophyllous vegetation", + "324", + "324 Transitional woodland-shrub", + "331", + "331 Beaches, dunes, sands", + "332", + "332 Bare rocks", + "333", + "333 Sparsely vegetated areas", + "334", + "334 Burnt areas", + "335", + "335 Glaciers and perpetual snow", + "411", + "411 Inland marshes", + "412", + "412 Peat bogs", + "421", + "421 Salt marshes", + "422", + "422 Salines", + "423", + "423 Intertidal flats", + "511", + "511 Water courses", + "512", + "512 Water bodies", + "521", + "521 Coastal lagoons", + "522", + "522 Estuaries", + "523", + "523 Sea and ocean" + ], + "IndexedColors": [ + 0.23529411764705882, + 0.23529411764705882, + 0.23529411764705882, + 0.5490196078431373, + 0.37254901960784315, + 0.35294117647058826, + 0.6588235294117647, + 0.5529411764705883, + 0.5254901960784314, + 0.5568627450980392, + 0.5411764705882353, + 0.5882352941176471, + 0.43137254901960786, + 0.43137254901960786, + 0.43137254901960786, + 0.49019607843137253, + 0.5294117647058824, + 0.5882352941176471, + 0.6588235294117647, + 0.6588235294117647, + 0.6588235294117647, + 0.7058823529411765, + 0.6078431372549019, + 0.47058823529411764, + 0.5411764705882353, + 0.49019607843137253, + 0.39215686274509803, + 0.7372549019607844, + 0.6823529411764706, + 0.596078431372549, + 0.5254901960784314, + 0.7215686274509804, + 0.43137254901960786, + 0.6392156862745098, + 0.8117647058823529, + 0.49411764705882355, + 0.9019607843137255, + 0.8352941176470589, + 0.5490196078431373, + 0.8509803921568627, + 0.8, + 0.43137254901960786, + 0.803921568627451, + 0.8470588235294118, + 0.6039215686274509, + 0.611764705882353, + 0.47843137254901963, + 0.26666666666666666, + 0.6627450980392157, + 0.6784313725490196, + 0.3686274509803922, + 0.5607843137254902, + 0.6039215686274509, + 0.3215686274509804, + 0.7058823529411765, + 0.8313725490196079, + 0.47843137254901963, + 0.8627450980392157, + 0.8117647058823529, + 0.5803921568627451, + 0.8235294117647058, + 0.7686274509803922, + 0.48627450980392156, + 0.7372549019607844, + 0.7686274509803922, + 0.49411764705882355, + 0.6627450980392157, + 0.7098039215686275, + 0.4549019607843137, + 0.30980392156862746, + 0.5607843137254902, + 0.24705882352941178, + 0.12156862745098039, + 0.35294117647058826, + 0.1803921568627451, + 0.20784313725490197, + 0.4549019607843137, + 0.2196078431372549, + 0.7607843137254902, + 0.8392156862745098, + 0.5411764705882353, + 0.6549019607843137, + 0.7686274509803922, + 0.4980392156862745, + 0.5411764705882353, + 0.5882352941176471, + 0.34509803921568627, + 0.49411764705882355, + 0.6509803921568628, + 0.35294117647058826, + 0.9137254901960784, + 0.8666666666666667, + 0.6980392156862745, + 0.5607843137254902, + 0.5607843137254902, + 0.5607843137254902, + 0.7254901960784313, + 0.7215686274509804, + 0.6274509803921569, + 0.29411764705882354, + 0.24705882352941178, + 0.22745098039215686, + 0.9333333333333333, + 0.9647058823529412, + 0.984313725490196, + 0.43137254901960786, + 0.5803921568627451, + 0.4392156862745098, + 0.5411764705882353, + 0.4, + 0.25882352941176473, + 0.4980392156862745, + 0.6235294117647059, + 0.5490196078431373, + 0.8470588235294118, + 0.8392156862745098, + 0.8, + 0.7019607843137254, + 0.7411764705882353, + 0.7019607843137254, + 0.2980392156862745, + 0.5568627450980392, + 0.7686274509803922, + 0.24313725490196078, + 0.4823529411764706, + 0.7137254901960784, + 0.3568627450980392, + 0.5764705882352941, + 0.7019607843137254, + 0.3176470588235294, + 0.5333333333333333, + 0.7058823529411765, + 0.16862745098039217, + 0.36470588235294116, + 0.5568627450980392 + ], + "NanColor": [ + 1.0, + 0.0, + 1.0 + ] + } +] diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/logs/analysis.txt b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/analysis.txt new file mode 100644 index 00000000..4fcecef6 --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/analysis.txt @@ -0,0 +1,57 @@ +CORINE polygons meeting the domain: 79; clipped area 304.909550 km2, domain 304.909550 km2 + +== 10 m == +vtkPolyDataReader: 92053 cells = 8882 lines + 83171 triangles; active SCALARS feature_mask +land_cover_code: 92053 values (one per cell: True); nonzero on lines: 0 +FieldData land_cover_codes: 'CORINE Land Cover level-3 code, attribute Code_18, map corine; 0 = in no polygon, and every constraint line' +mesh area 304.909551 km2; triangles with code 0: 0 +| code | mesh km2 | mesh share % | CORINE clipped share % | diff (pp) | 22's table % | +| 333 | 129.093253 | 42.33821 | 42.33821 | -0.000001 | 42.34 | +| 332 | 65.011569 | 21.32159 | 21.32159 | +0.000000 | 21.32 | +| 322 | 51.849126 | 17.00476 | 17.00476 | +0.000000 | 17.0 | +| 512 | 50.191119 | 16.46099 | 16.46099 | +0.000000 | 16.46 | +| 335 | 7.452202 | 2.44407 | 2.44407 | -0.000000 | 2.44 | +| 412 | 1.020748 | 0.33477 | 0.33477 | +0.000000 | 0.33 | +| 142 | 0.291534 | 0.09561 | 0.09561 | -0.000000 | 0.1 | +largest share difference: 0.000001 pp; codes equal: True +I1: 116398 interior unconstrained edges, 0 with differing codes +I2: oracle checks 83171 of 83171 triangles (r > 0.003 m); disagreements 0; disagreements on the 0 unchecked: 0 +union-find: 5 hook rounds; 124 components; equal to landcover.regions: True +ASCII run codes equal to binary run1: True +quality: worst angle 0.00408 deg, max degree 16, max error 9.999732907715497 (within: True), Delaunay checked/ambiguous/violations 116398/1095/0 +time land cover: median 0.071 s (runs 0.075, 0.071, 0.070) +time refine: median 0.047 s (runs 0.048, 0.047, 0.046) +time features clip: median 1.460 s (runs 1.498, 1.460, 1.419) +time total: median 2.075 s (runs 2.534, 2.075, 2.021) +wall: median 2.34 s (runs 3.10, 2.34, 2.29); max RSS median 719 MB +run1 stderr: land cover: 124 regions, 0 outside every polygon, 0 in more than one, 0 thinner than the snap +run2 stderr: land cover: 124 regions, 0 outside every polygon, 0 in more than one, 0 thinner than the snap +run3 stderr: land cover: 124 regions, 0 outside every polygon, 0 in more than one, 0 thinner than the snap + +== 1 m == +vtkPolyDataReader: 1159508 cells = 14074 lines + 1145434 triangles; active SCALARS feature_mask +land_cover_code: 1159508 values (one per cell: True); nonzero on lines: 0 +FieldData land_cover_codes: 'CORINE Land Cover level-3 code, attribute Code_18, map corine; 0 = in no polygon, and every constraint line' +mesh area 304.909551 km2; triangles with code 0: 0 +| code | mesh km2 | mesh share % | CORINE clipped share % | diff (pp) | 22's table % | +| 333 | 129.093253 | 42.33821 | 42.33821 | -0.000001 | 42.34 | +| 332 | 65.011569 | 21.32159 | 21.32159 | +0.000000 | 21.32 | +| 322 | 51.849126 | 17.00476 | 17.00476 | +0.000000 | 17.0 | +| 512 | 50.191119 | 16.46099 | 16.46099 | +0.000000 | 16.46 | +| 335 | 7.452202 | 2.44407 | 2.44407 | -0.000000 | 2.44 | +| 412 | 1.020748 | 0.33477 | 0.33477 | +0.000000 | 0.33 | +| 142 | 0.291534 | 0.09561 | 0.09561 | -0.000000 | 0.1 | +largest share difference: 0.000001 pp; codes equal: True +I1: 1705117 interior unconstrained edges, 0 with differing codes +I2: oracle checks 1145434 of 1145434 triangles (r > 0.003 m); disagreements 0; disagreements on the 0 unchecked: 0 +union-find: 7 hook rounds; 124 components; equal to landcover.regions: True +ASCII run codes equal to binary run1: True +quality: worst angle 0.00408 deg, max degree 29, max error 1.0 (within: True), Delaunay checked/ambiguous/violations 1705117/162591/0 +time land cover: median 1.296 s (runs 1.303, 1.296, 1.263) +time refine: median 0.559 s (runs 0.573, 0.559, 0.537) +time features clip: median 1.447 s (runs 1.447, 1.463, 1.436) +time total: median 3.887 s (runs 3.887, 3.952, 3.800) +wall: median 4.36 s (runs 4.36, 4.41, 4.26); max RSS median 859 MB +run1 stderr: land cover: 124 regions, 0 outside every polygon, 0 in more than one, 0 thinner than the snap +run2 stderr: land cover: 124 regions, 0 outside every polygon, 0 in more than one, 0 thinner than the snap +run3 stderr: land cover: 124 regions, 0 outside every polygon, 0 in more than one, 0 thinner than the snap diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/logs/commit.txt b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/commit.txt new file mode 100644 index 00000000..fa3accb1 --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/commit.txt @@ -0,0 +1,3 @@ +f7d5f14fe34e64cdfd77887394588ad87a0f8cef +190587e7ba0ea96a74a8621767d7d20257c3968126a710db1364e1195afb27d2 .venv/lib/python3.14/site-packages/tin_engine/_core.cpython-314-darwin.so +190587e7ba0ea96a74a8621767d7d20257c3968126a710db1364e1195afb27d2 build-pyext/_core.cpython-314-darwin.so diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/logs/exits.txt b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/exits.txt new file mode 100644 index 00000000..fcbe9f23 --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/exits.txt @@ -0,0 +1,9 @@ +lc_t10_run1 exit=0 +lc_t1_run1 exit=0 +lc_t10_run2 exit=0 +lc_t1_run2 exit=0 +lc_t10_run3 exit=0 +lc_t1_run3 exit=0 +lc_t10_ascii exit=0 +lc_t1_ascii exit=0 +palette exit=0 diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_ascii.err b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_ascii.err new file mode 100644 index 00000000..4ad06393 --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_ascii.err @@ -0,0 +1,23 @@ +79 features kept, 342 dropped outside, 35 clipped, 0 empty skipped +15846 input vertices, 8242 noded vertices +15989 start quality nodes inserted, 2442 start quality skips, 534 constraint feet, 0 feet refused, 17 rounds, 17879 points inserted, 39582 flips, 83171 triangles, achieved max error 9.999732907715497 m, 0 valid DEM nodes not covered, 15614 start triangles, 8242 start vertices off-node, 0 vertices without data dropped +mosaic of 2 tiles, 2195 x 3314 nodes +land cover: 124 regions, 0 outside every polygon, 0 in more than one, 0 thinner than the snap + 2.49 real 2.83 user 0.20 sys + 716095488 maximum resident set size + 0 average shared memory size + 0 average unshared data size + 0 average unshared stack size + 62781 page reclaims + 44 page faults + 0 swaps + 0 block input operations + 0 block output operations + 0 messages sent + 0 messages received + 0 signals received + 5 voluntary context switches + 541 involuntary context switches + 42995350613 instructions retired + 9676800505 cycles elapsed + 693666800 peak memory footprint diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_ascii.out b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_ascii.out new file mode 100644 index 00000000..c7635b28 --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_ascii.out @@ -0,0 +1,2 @@ +/private/tmp/claude-501/-Users-skavhaug-projects-rasputin/38487caf-f56e-46a3-b05b-867af1fb1619/scratchpad/lc/lc_t10_ascii.vtk +/Users/skavhaug/projects/rasputin/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_ascii.stats.md diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_ascii.stats.md b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_ascii.stats.md new file mode 100644 index 00000000..dfdda6f4 --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_ascii.stats.md @@ -0,0 +1,67 @@ +# rasputin mesh — statistics + +`rasputin mesh --dem ../rasputin_data/DTM10_UTM33_20260925 --domain docs/benchmarks/2026-09-29/bygdin-landcover/bygdin_reduced_t20.geojson --tolerance 10 --features ../rasputin_data/corine2018_dtm10_utm33.gpkg --features-layer corine2018 --features-map corine --ascii --out /private/tmp/claude-501/-Users-skavhaug-projects-rasputin/38487caf-f56e-46a3-b05b-867af1fb1619/scratchpad/lc/lc_t10_ascii.vtk --stats docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_ascii.stats.md` + +## Sizes + +| item | count | +|---|---| +| DEM nodes | 2195 × 3314 (10 m) | +| domain vertices | 740 (1 ring, 0 holes) | +| start vertices | 8242 | +| start triangles | 15614 | +| output vertices | 42110 | +| output triangles | 83171 | +| constraint edges | 8882 | +| vertices without data dropped | 0 | +| lc_t10_ascii.vtk | 5.0 MB | + +## Quality (plan view, x/y) + +| metric | median | < 1° | < 10° | worst | +|---|---|---|---|---| +| minimum angle | 36.87° | 0.07 % | 1.70 % | 0.00408° | + +| metric | median | p99 | max | ≥ 12 | ≥ 20 | +|---|---|---|---|---|---| +| vertex degree (triangles) | 6 | 9 | 16 | 31 | 0 | + +## Refinement + +| tolerance | achieved max error | rounds | inserted | carved | flips | uncovered | quality inserted | quality skipped | feet | +|---|---|---|---|---|---|---|---|---|---| +| 10 m | 9.999732907715497 m | 17 | 17879 | 0 | 39582 | 0 | 15989 | 2442 | 534 | + +## Timings + +Wall clock, `time.perf_counter_ns` (Python) and `std::chrono::steady_clock` +(inside `refine`), one run, no warm-up. Total is the `mesh` command body, from +argument checks to the last file written; interpreter start-up and imports are +not in it. Threads: 10 (hardware concurrency). + +| phase | seconds | share | +|---|---|---| +| domain read | 0.002 | 0.1 % | +| decode | 0.307 | 13.8 % | +| features read | 0.119 | 5.3 % | +| features clip | 1.436 | 64.2 % | +| start mesh: build | 0.000 | 0.0 % | +| start mesh: node | 0.017 | 0.8 % | +| start mesh: triangulate | 0.005 | 0.2 % | +| start mesh: constraint edges | 0.016 | 0.7 % | +| refine | 0.046 | 2.1 % | +| refine: legalise start | 0.000 | 0.0 % | +| refine: start quality | 0.015 | 0.7 % | +| refine: scan (parallel) | 0.013 | 0.6 % | +| refine: split + flip (serial) | 0.009 | 0.4 % | +| refine: setup + output | 0.008 | 0.4 % | +| trim | 0.003 | 0.1 % | +| land cover | 0.070 | 3.1 % | +| write: encode | 0.207 | 9.2 % | +| write: disk | 0.003 | 0.1 % | +| other | 0.004 | 0.2 % | +| **total** | **2.236** | **100 %** | + +Sub-rows sum to their parent and are not added to the total. + +Statistics computed in 0.015 s, not included above. diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_run1.err b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_run1.err new file mode 100644 index 00000000..93f78526 --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_run1.err @@ -0,0 +1,23 @@ +79 features kept, 342 dropped outside, 35 clipped, 0 empty skipped +15846 input vertices, 8242 noded vertices +15989 start quality nodes inserted, 2442 start quality skips, 534 constraint feet, 0 feet refused, 17 rounds, 17879 points inserted, 39582 flips, 83171 triangles, achieved max error 9.999732907715497 m, 0 valid DEM nodes not covered, 15614 start triangles, 8242 start vertices off-node, 0 vertices without data dropped +mosaic of 2 tiles, 2195 x 3314 nodes +land cover: 124 regions, 0 outside every polygon, 0 in more than one, 0 thinner than the snap + 3.10 real 2.77 user 0.36 sys + 768819200 maximum resident set size + 0 average shared memory size + 0 average unshared data size + 0 average unshared stack size + 64750 page reclaims + 1438 page faults + 0 swaps + 0 block input operations + 0 block output operations + 0 messages sent + 0 messages received + 0 signals received + 2023 voluntary context switches + 2727 involuntary context switches + 40098720267 instructions retired + 9746241268 cycles elapsed + 746423328 peak memory footprint diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_run1.out b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_run1.out new file mode 100644 index 00000000..d597813a --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_run1.out @@ -0,0 +1,2 @@ +/private/tmp/claude-501/-Users-skavhaug-projects-rasputin/38487caf-f56e-46a3-b05b-867af1fb1619/scratchpad/lc/lc_t10_run1.vtk +/Users/skavhaug/projects/rasputin/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_run1.stats.md diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_run1.stats.md b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_run1.stats.md new file mode 100644 index 00000000..a1dc3fac --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_run1.stats.md @@ -0,0 +1,67 @@ +# rasputin mesh — statistics + +`rasputin mesh --dem ../rasputin_data/DTM10_UTM33_20260925 --domain docs/benchmarks/2026-09-29/bygdin-landcover/bygdin_reduced_t20.geojson --tolerance 10 --features ../rasputin_data/corine2018_dtm10_utm33.gpkg --features-layer corine2018 --features-map corine --binary --out /private/tmp/claude-501/-Users-skavhaug-projects-rasputin/38487caf-f56e-46a3-b05b-867af1fb1619/scratchpad/lc/lc_t10_run1.vtk --stats docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_run1.stats.md` + +## Sizes + +| item | count | +|---|---| +| DEM nodes | 2195 × 3314 (10 m) | +| domain vertices | 740 (1 ring, 0 holes) | +| start vertices | 8242 | +| start triangles | 15614 | +| output vertices | 42110 | +| output triangles | 83171 | +| constraint edges | 8882 | +| vertices without data dropped | 0 | +| lc_t10_run1.vtk | 3.7 MB | + +## Quality (plan view, x/y) + +| metric | median | < 1° | < 10° | worst | +|---|---|---|---|---| +| minimum angle | 36.87° | 0.07 % | 1.70 % | 0.00408° | + +| metric | median | p99 | max | ≥ 12 | ≥ 20 | +|---|---|---|---|---|---| +| vertex degree (triangles) | 6 | 9 | 16 | 31 | 0 | + +## Refinement + +| tolerance | achieved max error | rounds | inserted | carved | flips | uncovered | quality inserted | quality skipped | feet | +|---|---|---|---|---|---|---|---|---|---| +| 10 m | 9.999732907715497 m | 17 | 17879 | 0 | 39582 | 0 | 15989 | 2442 | 534 | + +## Timings + +Wall clock, `time.perf_counter_ns` (Python) and `std::chrono::steady_clock` +(inside `refine`), one run, no warm-up. Total is the `mesh` command body, from +argument checks to the last file written; interpreter start-up and imports are +not in it. Threads: 10 (hardware concurrency). + +| phase | seconds | share | +|---|---|---| +| domain read | 0.011 | 0.4 % | +| decode | 0.621 | 24.5 % | +| features read | 0.227 | 9.0 % | +| features clip | 1.498 | 59.1 % | +| start mesh: build | 0.000 | 0.0 % | +| start mesh: node | 0.017 | 0.7 % | +| start mesh: triangulate | 0.005 | 0.2 % | +| start mesh: constraint edges | 0.016 | 0.6 % | +| refine | 0.048 | 1.9 % | +| refine: legalise start | 0.000 | 0.0 % | +| refine: start quality | 0.016 | 0.6 % | +| refine: scan (parallel) | 0.014 | 0.5 % | +| refine: split + flip (serial) | 0.010 | 0.4 % | +| refine: setup + output | 0.009 | 0.3 % | +| trim | 0.004 | 0.2 % | +| land cover | 0.075 | 3.0 % | +| write: encode | 0.002 | 0.1 % | +| write: disk | 0.001 | 0.0 % | +| other | 0.008 | 0.3 % | +| **total** | **2.534** | **100 %** | + +Sub-rows sum to their parent and are not added to the total. + +Statistics computed in 0.016 s, not included above. diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_run2.err b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_run2.err new file mode 100644 index 00000000..dd8772ea --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_run2.err @@ -0,0 +1,23 @@ +79 features kept, 342 dropped outside, 35 clipped, 0 empty skipped +15846 input vertices, 8242 noded vertices +15989 start quality nodes inserted, 2442 start quality skips, 534 constraint feet, 0 feet refused, 17 rounds, 17879 points inserted, 39582 flips, 83171 triangles, achieved max error 9.999732907715497 m, 0 valid DEM nodes not covered, 15614 start triangles, 8242 start vertices off-node, 0 vertices without data dropped +mosaic of 2 tiles, 2195 x 3314 nodes +land cover: 124 regions, 0 outside every polygon, 0 in more than one, 0 thinner than the snap + 2.34 real 2.66 user 0.21 sys + 719142912 maximum resident set size + 0 average shared memory size + 0 average unshared data size + 0 average unshared stack size + 64062 page reclaims + 109 page faults + 0 swaps + 0 block input operations + 0 block output operations + 0 messages sent + 0 messages received + 0 signals received + 51 voluntary context switches + 628 involuntary context switches + 39497552999 instructions retired + 9116605155 cycles elapsed + 696714224 peak memory footprint diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_run2.out b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_run2.out new file mode 100644 index 00000000..b676ba6c --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_run2.out @@ -0,0 +1,2 @@ +/private/tmp/claude-501/-Users-skavhaug-projects-rasputin/38487caf-f56e-46a3-b05b-867af1fb1619/scratchpad/lc/lc_t10_run2.vtk +/Users/skavhaug/projects/rasputin/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_run2.stats.md diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_run2.stats.md b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_run2.stats.md new file mode 100644 index 00000000..ca5be4fe --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_run2.stats.md @@ -0,0 +1,67 @@ +# rasputin mesh — statistics + +`rasputin mesh --dem ../rasputin_data/DTM10_UTM33_20260925 --domain docs/benchmarks/2026-09-29/bygdin-landcover/bygdin_reduced_t20.geojson --tolerance 10 --features ../rasputin_data/corine2018_dtm10_utm33.gpkg --features-layer corine2018 --features-map corine --binary --out /private/tmp/claude-501/-Users-skavhaug-projects-rasputin/38487caf-f56e-46a3-b05b-867af1fb1619/scratchpad/lc/lc_t10_run2.vtk --stats docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_run2.stats.md` + +## Sizes + +| item | count | +|---|---| +| DEM nodes | 2195 × 3314 (10 m) | +| domain vertices | 740 (1 ring, 0 holes) | +| start vertices | 8242 | +| start triangles | 15614 | +| output vertices | 42110 | +| output triangles | 83171 | +| constraint edges | 8882 | +| vertices without data dropped | 0 | +| lc_t10_run2.vtk | 3.7 MB | + +## Quality (plan view, x/y) + +| metric | median | < 1° | < 10° | worst | +|---|---|---|---|---| +| minimum angle | 36.87° | 0.07 % | 1.70 % | 0.00408° | + +| metric | median | p99 | max | ≥ 12 | ≥ 20 | +|---|---|---|---|---|---| +| vertex degree (triangles) | 6 | 9 | 16 | 31 | 0 | + +## Refinement + +| tolerance | achieved max error | rounds | inserted | carved | flips | uncovered | quality inserted | quality skipped | feet | +|---|---|---|---|---|---|---|---|---|---| +| 10 m | 9.999732907715497 m | 17 | 17879 | 0 | 39582 | 0 | 15989 | 2442 | 534 | + +## Timings + +Wall clock, `time.perf_counter_ns` (Python) and `std::chrono::steady_clock` +(inside `refine`), one run, no warm-up. Total is the `mesh` command body, from +argument checks to the last file written; interpreter start-up and imports are +not in it. Threads: 10 (hardware concurrency). + +| phase | seconds | share | +|---|---|---| +| domain read | 0.002 | 0.1 % | +| decode | 0.310 | 14.9 % | +| features read | 0.131 | 6.3 % | +| features clip | 1.460 | 70.4 % | +| start mesh: build | 0.000 | 0.0 % | +| start mesh: node | 0.017 | 0.8 % | +| start mesh: triangulate | 0.005 | 0.2 % | +| start mesh: constraint edges | 0.016 | 0.8 % | +| refine | 0.047 | 2.3 % | +| refine: legalise start | 0.000 | 0.0 % | +| refine: start quality | 0.015 | 0.7 % | +| refine: scan (parallel) | 0.015 | 0.7 % | +| refine: split + flip (serial) | 0.009 | 0.4 % | +| refine: setup + output | 0.008 | 0.4 % | +| trim | 0.003 | 0.1 % | +| land cover | 0.071 | 3.4 % | +| write: encode | 0.002 | 0.1 % | +| write: disk | 0.003 | 0.1 % | +| other | 0.007 | 0.3 % | +| **total** | **2.075** | **100 %** | + +Sub-rows sum to their parent and are not added to the total. + +Statistics computed in 0.015 s, not included above. diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_run3.err b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_run3.err new file mode 100644 index 00000000..e595180a --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_run3.err @@ -0,0 +1,23 @@ +79 features kept, 342 dropped outside, 35 clipped, 0 empty skipped +15846 input vertices, 8242 noded vertices +15989 start quality nodes inserted, 2442 start quality skips, 534 constraint feet, 0 feet refused, 17 rounds, 17879 points inserted, 39582 flips, 83171 triangles, achieved max error 9.999732907715497 m, 0 valid DEM nodes not covered, 15614 start triangles, 8242 start vertices off-node, 0 vertices without data dropped +mosaic of 2 tiles, 2195 x 3314 nodes +land cover: 124 regions, 0 outside every polygon, 0 in more than one, 0 thinner than the snap + 2.29 real 2.61 user 0.22 sys + 718979072 maximum resident set size + 0 average shared memory size + 0 average unshared data size + 0 average unshared stack size + 63059 page reclaims + 109 page faults + 0 swaps + 0 block input operations + 0 block output operations + 0 messages sent + 0 messages received + 0 signals received + 0 voluntary context switches + 607 involuntary context switches + 39539025182 instructions retired + 9016158015 cycles elapsed + 696550384 peak memory footprint diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_run3.out b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_run3.out new file mode 100644 index 00000000..7a094ceb --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_run3.out @@ -0,0 +1,2 @@ +/private/tmp/claude-501/-Users-skavhaug-projects-rasputin/38487caf-f56e-46a3-b05b-867af1fb1619/scratchpad/lc/lc_t10_run3.vtk +/Users/skavhaug/projects/rasputin/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_run3.stats.md diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_run3.stats.md b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_run3.stats.md new file mode 100644 index 00000000..72cd45c2 --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_run3.stats.md @@ -0,0 +1,67 @@ +# rasputin mesh — statistics + +`rasputin mesh --dem ../rasputin_data/DTM10_UTM33_20260925 --domain docs/benchmarks/2026-09-29/bygdin-landcover/bygdin_reduced_t20.geojson --tolerance 10 --features ../rasputin_data/corine2018_dtm10_utm33.gpkg --features-layer corine2018 --features-map corine --binary --out /private/tmp/claude-501/-Users-skavhaug-projects-rasputin/38487caf-f56e-46a3-b05b-867af1fb1619/scratchpad/lc/lc_t10_run3.vtk --stats docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t10_run3.stats.md` + +## Sizes + +| item | count | +|---|---| +| DEM nodes | 2195 × 3314 (10 m) | +| domain vertices | 740 (1 ring, 0 holes) | +| start vertices | 8242 | +| start triangles | 15614 | +| output vertices | 42110 | +| output triangles | 83171 | +| constraint edges | 8882 | +| vertices without data dropped | 0 | +| lc_t10_run3.vtk | 3.7 MB | + +## Quality (plan view, x/y) + +| metric | median | < 1° | < 10° | worst | +|---|---|---|---|---| +| minimum angle | 36.87° | 0.07 % | 1.70 % | 0.00408° | + +| metric | median | p99 | max | ≥ 12 | ≥ 20 | +|---|---|---|---|---|---| +| vertex degree (triangles) | 6 | 9 | 16 | 31 | 0 | + +## Refinement + +| tolerance | achieved max error | rounds | inserted | carved | flips | uncovered | quality inserted | quality skipped | feet | +|---|---|---|---|---|---|---|---|---|---| +| 10 m | 9.999732907715497 m | 17 | 17879 | 0 | 39582 | 0 | 15989 | 2442 | 534 | + +## Timings + +Wall clock, `time.perf_counter_ns` (Python) and `std::chrono::steady_clock` +(inside `refine`), one run, no warm-up. Total is the `mesh` command body, from +argument checks to the last file written; interpreter start-up and imports are +not in it. Threads: 10 (hardware concurrency). + +| phase | seconds | share | +|---|---|---| +| domain read | 0.002 | 0.1 % | +| decode | 0.313 | 15.5 % | +| features read | 0.122 | 6.1 % | +| features clip | 1.419 | 70.2 % | +| start mesh: build | 0.000 | 0.0 % | +| start mesh: node | 0.017 | 0.8 % | +| start mesh: triangulate | 0.005 | 0.3 % | +| start mesh: constraint edges | 0.016 | 0.8 % | +| refine | 0.046 | 2.3 % | +| refine: legalise start | 0.000 | 0.0 % | +| refine: start quality | 0.015 | 0.8 % | +| refine: scan (parallel) | 0.013 | 0.7 % | +| refine: split + flip (serial) | 0.009 | 0.4 % | +| refine: setup + output | 0.008 | 0.4 % | +| trim | 0.003 | 0.1 % | +| land cover | 0.070 | 3.5 % | +| write: encode | 0.002 | 0.1 % | +| write: disk | 0.001 | 0.0 % | +| other | 0.004 | 0.2 % | +| **total** | **2.021** | **100 %** | + +Sub-rows sum to their parent and are not added to the total. + +Statistics computed in 0.014 s, not included above. diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_ascii.err b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_ascii.err new file mode 100644 index 00000000..87c83302 --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_ascii.err @@ -0,0 +1,23 @@ +79 features kept, 342 dropped outside, 35 clipped, 0 empty skipped +15846 input vertices, 8242 noded vertices +15989 start quality nodes inserted, 2442 start quality skips, 5726 constraint feet, 0 feet refused, 24 rounds, 549527 points inserted, 1152221 flips, 1145434 triangles, achieved max error 1 m, 0 valid DEM nodes not covered, 15614 start triangles, 8242 start vertices off-node, 0 vertices without data dropped +mosaic of 2 tiles, 2195 x 3314 nodes +land cover: 124 regions, 0 outside every polygon, 0 in more than one, 0 thinner than the snap + 7.02 real 7.80 user 0.43 sys + 1045397504 maximum resident set size + 0 average shared memory size + 0 average unshared data size + 0 average unshared stack size + 174571 page reclaims + 44 page faults + 0 swaps + 0 block input operations + 0 block output operations + 0 messages sent + 0 messages received + 0 signals received + 21 voluntary context switches + 764 involuntary context switches + 99887680145 instructions retired + 26140403929 cycles elapsed + 866092352 peak memory footprint diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_ascii.out b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_ascii.out new file mode 100644 index 00000000..2766df19 --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_ascii.out @@ -0,0 +1,2 @@ +/private/tmp/claude-501/-Users-skavhaug-projects-rasputin/38487caf-f56e-46a3-b05b-867af1fb1619/scratchpad/lc/lc_t1_ascii.vtk +/Users/skavhaug/projects/rasputin/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_ascii.stats.md diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_ascii.stats.md b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_ascii.stats.md new file mode 100644 index 00000000..4b1c6039 --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_ascii.stats.md @@ -0,0 +1,67 @@ +# rasputin mesh — statistics + +`rasputin mesh --dem ../rasputin_data/DTM10_UTM33_20260925 --domain docs/benchmarks/2026-09-29/bygdin-landcover/bygdin_reduced_t20.geojson --tolerance 1 --features ../rasputin_data/corine2018_dtm10_utm33.gpkg --features-layer corine2018 --features-map corine --ascii --out /private/tmp/claude-501/-Users-skavhaug-projects-rasputin/38487caf-f56e-46a3-b05b-867af1fb1619/scratchpad/lc/lc_t1_ascii.vtk --stats docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_ascii.stats.md` + +## Sizes + +| item | count | +|---|---| +| DEM nodes | 2195 × 3314 (10 m) | +| domain vertices | 740 (1 ring, 0 holes) | +| start vertices | 8242 | +| start triangles | 15614 | +| output vertices | 573758 | +| output triangles | 1145434 | +| constraint edges | 14074 | +| vertices without data dropped | 0 | +| lc_t1_ascii.vtk | 69.0 MB | + +## Quality (plan view, x/y) + +| metric | median | < 1° | < 10° | worst | +|---|---|---|---|---| +| minimum angle | 45.00° | 0.10 % | 1.36 % | 0.00408° | + +| metric | median | p99 | max | ≥ 12 | ≥ 20 | +|---|---|---|---|---|---| +| vertex degree (triangles) | 6 | 9 | 29 | 1377 | 37 | + +## Refinement + +| tolerance | achieved max error | rounds | inserted | carved | flips | uncovered | quality inserted | quality skipped | feet | +|---|---|---|---|---|---|---|---|---|---| +| 1 m | 1 m | 24 | 549527 | 0 | 1152221 | 0 | 15989 | 2442 | 5726 | + +## Timings + +Wall clock, `time.perf_counter_ns` (Python) and `std::chrono::steady_clock` +(inside `refine`), one run, no warm-up. Total is the `mesh` command body, from +argument checks to the last file written; interpreter start-up and imports are +not in it. Threads: 10 (hardware concurrency). + +| phase | seconds | share | +|---|---|---| +| domain read | 0.002 | 0.0 % | +| decode | 0.310 | 4.7 % | +| features read | 0.118 | 1.8 % | +| features clip | 1.418 | 21.6 % | +| start mesh: build | 0.000 | 0.0 % | +| start mesh: node | 0.017 | 0.3 % | +| start mesh: triangulate | 0.005 | 0.1 % | +| start mesh: constraint edges | 0.016 | 0.2 % | +| refine | 0.536 | 8.2 % | +| refine: legalise start | 0.000 | 0.0 % | +| refine: start quality | 0.015 | 0.2 % | +| refine: scan (parallel) | 0.094 | 1.4 % | +| refine: split + flip (serial) | 0.377 | 5.7 % | +| refine: setup + output | 0.049 | 0.7 % | +| trim | 0.042 | 0.6 % | +| land cover | 1.256 | 19.1 % | +| write: encode | 2.810 | 42.8 % | +| write: disk | 0.036 | 0.5 % | +| other | 0.005 | 0.1 % | +| **total** | **6.571** | **100 %** | + +Sub-rows sum to their parent and are not added to the total. + +Statistics computed in 0.201 s, not included above. diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_run1.err b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_run1.err new file mode 100644 index 00000000..49cfde19 --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_run1.err @@ -0,0 +1,23 @@ +79 features kept, 342 dropped outside, 35 clipped, 0 empty skipped +15846 input vertices, 8242 noded vertices +15989 start quality nodes inserted, 2442 start quality skips, 5726 constraint feet, 0 feet refused, 24 rounds, 549527 points inserted, 1152221 flips, 1145434 triangles, achieved max error 1 m, 0 valid DEM nodes not covered, 15614 start triangles, 8242 start vertices off-node, 0 vertices without data dropped +mosaic of 2 tiles, 2195 x 3314 nodes +land cover: 124 regions, 0 outside every polygon, 0 in more than one, 0 thinner than the snap + 4.36 real 5.25 user 0.31 sys + 855375872 maximum resident set size + 0 average shared memory size + 0 average unshared data size + 0 average unshared stack size + 95870 page reclaims + 44 page faults + 0 swaps + 0 block input operations + 0 block output operations + 0 messages sent + 0 messages received + 0 signals received + 21 voluntary context switches + 844 involuntary context switches + 54320270403 instructions retired + 17478869271 cycles elapsed + 696370160 peak memory footprint diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_run1.out b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_run1.out new file mode 100644 index 00000000..82e23d06 --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_run1.out @@ -0,0 +1,2 @@ +/private/tmp/claude-501/-Users-skavhaug-projects-rasputin/38487caf-f56e-46a3-b05b-867af1fb1619/scratchpad/lc/lc_t1_run1.vtk +/Users/skavhaug/projects/rasputin/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_run1.stats.md diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_run1.stats.md b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_run1.stats.md new file mode 100644 index 00000000..36aa4f43 --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_run1.stats.md @@ -0,0 +1,67 @@ +# rasputin mesh — statistics + +`rasputin mesh --dem ../rasputin_data/DTM10_UTM33_20260925 --domain docs/benchmarks/2026-09-29/bygdin-landcover/bygdin_reduced_t20.geojson --tolerance 1 --features ../rasputin_data/corine2018_dtm10_utm33.gpkg --features-layer corine2018 --features-map corine --binary --out /private/tmp/claude-501/-Users-skavhaug-projects-rasputin/38487caf-f56e-46a3-b05b-867af1fb1619/scratchpad/lc/lc_t1_run1.vtk --stats docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_run1.stats.md` + +## Sizes + +| item | count | +|---|---| +| DEM nodes | 2195 × 3314 (10 m) | +| domain vertices | 740 (1 ring, 0 holes) | +| start vertices | 8242 | +| start triangles | 15614 | +| output vertices | 573758 | +| output triangles | 1145434 | +| constraint edges | 14074 | +| vertices without data dropped | 0 | +| lc_t1_run1.vtk | 48.5 MB | + +## Quality (plan view, x/y) + +| metric | median | < 1° | < 10° | worst | +|---|---|---|---|---| +| minimum angle | 45.00° | 0.10 % | 1.36 % | 0.00408° | + +| metric | median | p99 | max | ≥ 12 | ≥ 20 | +|---|---|---|---|---|---| +| vertex degree (triangles) | 6 | 9 | 29 | 1377 | 37 | + +## Refinement + +| tolerance | achieved max error | rounds | inserted | carved | flips | uncovered | quality inserted | quality skipped | feet | +|---|---|---|---|---|---|---|---|---|---| +| 1 m | 1 m | 24 | 549527 | 0 | 1152221 | 0 | 15989 | 2442 | 5726 | + +## Timings + +Wall clock, `time.perf_counter_ns` (Python) and `std::chrono::steady_clock` +(inside `refine`), one run, no warm-up. Total is the `mesh` command body, from +argument checks to the last file written; interpreter start-up and imports are +not in it. Threads: 10 (hardware concurrency). + +| phase | seconds | share | +|---|---|---| +| domain read | 0.002 | 0.1 % | +| decode | 0.310 | 8.0 % | +| features read | 0.118 | 3.0 % | +| features clip | 1.447 | 37.2 % | +| start mesh: build | 0.000 | 0.0 % | +| start mesh: node | 0.017 | 0.4 % | +| start mesh: triangulate | 0.005 | 0.1 % | +| start mesh: constraint edges | 0.016 | 0.4 % | +| refine | 0.573 | 14.7 % | +| refine: legalise start | 0.000 | 0.0 % | +| refine: start quality | 0.016 | 0.4 % | +| refine: scan (parallel) | 0.097 | 2.5 % | +| refine: split + flip (serial) | 0.409 | 10.5 % | +| refine: setup + output | 0.050 | 1.3 % | +| trim | 0.045 | 1.2 % | +| land cover | 1.303 | 33.5 % | +| write: encode | 0.017 | 0.4 % | +| write: disk | 0.027 | 0.7 % | +| other | 0.005 | 0.1 % | +| **total** | **3.887** | **100 %** | + +Sub-rows sum to their parent and are not added to the total. + +Statistics computed in 0.202 s, not included above. diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_run2.err b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_run2.err new file mode 100644 index 00000000..e323a0f1 --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_run2.err @@ -0,0 +1,23 @@ +79 features kept, 342 dropped outside, 35 clipped, 0 empty skipped +15846 input vertices, 8242 noded vertices +15989 start quality nodes inserted, 2442 start quality skips, 5726 constraint feet, 0 feet refused, 24 rounds, 549527 points inserted, 1152221 flips, 1145434 triangles, achieved max error 1 m, 0 valid DEM nodes not covered, 15614 start triangles, 8242 start vertices off-node, 0 vertices without data dropped +mosaic of 2 tiles, 2195 x 3314 nodes +land cover: 124 regions, 0 outside every polygon, 0 in more than one, 0 thinner than the snap + 4.41 real 5.25 user 0.35 sys + 876838912 maximum resident set size + 0 average shared memory size + 0 average unshared data size + 0 average unshared stack size + 100890 page reclaims + 44 page faults + 0 swaps + 0 block input operations + 0 block output operations + 0 messages sent + 0 messages received + 0 signals received + 16 voluntary context switches + 1208 involuntary context switches + 54424424256 instructions retired + 17620119152 cycles elapsed + 699614192 peak memory footprint diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_run2.out b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_run2.out new file mode 100644 index 00000000..49173035 --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_run2.out @@ -0,0 +1,2 @@ +/private/tmp/claude-501/-Users-skavhaug-projects-rasputin/38487caf-f56e-46a3-b05b-867af1fb1619/scratchpad/lc/lc_t1_run2.vtk +/Users/skavhaug/projects/rasputin/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_run2.stats.md diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_run2.stats.md b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_run2.stats.md new file mode 100644 index 00000000..2c2e8f9e --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_run2.stats.md @@ -0,0 +1,67 @@ +# rasputin mesh — statistics + +`rasputin mesh --dem ../rasputin_data/DTM10_UTM33_20260925 --domain docs/benchmarks/2026-09-29/bygdin-landcover/bygdin_reduced_t20.geojson --tolerance 1 --features ../rasputin_data/corine2018_dtm10_utm33.gpkg --features-layer corine2018 --features-map corine --binary --out /private/tmp/claude-501/-Users-skavhaug-projects-rasputin/38487caf-f56e-46a3-b05b-867af1fb1619/scratchpad/lc/lc_t1_run2.vtk --stats docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_run2.stats.md` + +## Sizes + +| item | count | +|---|---| +| DEM nodes | 2195 × 3314 (10 m) | +| domain vertices | 740 (1 ring, 0 holes) | +| start vertices | 8242 | +| start triangles | 15614 | +| output vertices | 573758 | +| output triangles | 1145434 | +| constraint edges | 14074 | +| vertices without data dropped | 0 | +| lc_t1_run2.vtk | 48.5 MB | + +## Quality (plan view, x/y) + +| metric | median | < 1° | < 10° | worst | +|---|---|---|---|---| +| minimum angle | 45.00° | 0.10 % | 1.36 % | 0.00408° | + +| metric | median | p99 | max | ≥ 12 | ≥ 20 | +|---|---|---|---|---|---| +| vertex degree (triangles) | 6 | 9 | 29 | 1377 | 37 | + +## Refinement + +| tolerance | achieved max error | rounds | inserted | carved | flips | uncovered | quality inserted | quality skipped | feet | +|---|---|---|---|---|---|---|---|---|---| +| 1 m | 1 m | 24 | 549527 | 0 | 1152221 | 0 | 15989 | 2442 | 5726 | + +## Timings + +Wall clock, `time.perf_counter_ns` (Python) and `std::chrono::steady_clock` +(inside `refine`), one run, no warm-up. Total is the `mesh` command body, from +argument checks to the last file written; interpreter start-up and imports are +not in it. Threads: 10 (hardware concurrency). + +| phase | seconds | share | +|---|---|---| +| domain read | 0.002 | 0.1 % | +| decode | 0.360 | 9.1 % | +| features read | 0.139 | 3.5 % | +| features clip | 1.463 | 37.0 % | +| start mesh: build | 0.000 | 0.0 % | +| start mesh: node | 0.017 | 0.4 % | +| start mesh: triangulate | 0.005 | 0.1 % | +| start mesh: constraint edges | 0.016 | 0.4 % | +| refine | 0.559 | 14.1 % | +| refine: legalise start | 0.000 | 0.0 % | +| refine: start quality | 0.015 | 0.4 % | +| refine: scan (parallel) | 0.098 | 2.5 % | +| refine: split + flip (serial) | 0.393 | 9.9 % | +| refine: setup + output | 0.053 | 1.3 % | +| trim | 0.046 | 1.2 % | +| land cover | 1.296 | 32.8 % | +| write: encode | 0.017 | 0.4 % | +| write: disk | 0.025 | 0.6 % | +| other | 0.006 | 0.2 % | +| **total** | **3.952** | **100 %** | + +Sub-rows sum to their parent and are not added to the total. + +Statistics computed in 0.201 s, not included above. diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_run3.err b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_run3.err new file mode 100644 index 00000000..44e4aa9c --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_run3.err @@ -0,0 +1,23 @@ +79 features kept, 342 dropped outside, 35 clipped, 0 empty skipped +15846 input vertices, 8242 noded vertices +15989 start quality nodes inserted, 2442 start quality skips, 5726 constraint feet, 0 feet refused, 24 rounds, 549527 points inserted, 1152221 flips, 1145434 triangles, achieved max error 1 m, 0 valid DEM nodes not covered, 15614 start triangles, 8242 start vertices off-node, 0 vertices without data dropped +mosaic of 2 tiles, 2195 x 3314 nodes +land cover: 124 regions, 0 outside every polygon, 0 in more than one, 0 thinner than the snap + 4.26 real 5.19 user 0.28 sys + 859078656 maximum resident set size + 0 average shared memory size + 0 average unshared data size + 0 average unshared stack size + 99904 page reclaims + 44 page faults + 0 swaps + 0 block input operations + 0 block output operations + 0 messages sent + 0 messages received + 0 signals received + 17 voluntary context switches + 786 involuntary context switches + 54294758790 instructions retired + 17269081294 cycles elapsed + 695993328 peak memory footprint diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_run3.out b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_run3.out new file mode 100644 index 00000000..2a9cb9c0 --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_run3.out @@ -0,0 +1,2 @@ +/private/tmp/claude-501/-Users-skavhaug-projects-rasputin/38487caf-f56e-46a3-b05b-867af1fb1619/scratchpad/lc/lc_t1_run3.vtk +/Users/skavhaug/projects/rasputin/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_run3.stats.md diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_run3.stats.md b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_run3.stats.md new file mode 100644 index 00000000..f8fc38f3 --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_run3.stats.md @@ -0,0 +1,67 @@ +# rasputin mesh — statistics + +`rasputin mesh --dem ../rasputin_data/DTM10_UTM33_20260925 --domain docs/benchmarks/2026-09-29/bygdin-landcover/bygdin_reduced_t20.geojson --tolerance 1 --features ../rasputin_data/corine2018_dtm10_utm33.gpkg --features-layer corine2018 --features-map corine --binary --out /private/tmp/claude-501/-Users-skavhaug-projects-rasputin/38487caf-f56e-46a3-b05b-867af1fb1619/scratchpad/lc/lc_t1_run3.vtk --stats docs/benchmarks/2026-09-29/bygdin-landcover/logs/lc_t1_run3.stats.md` + +## Sizes + +| item | count | +|---|---| +| DEM nodes | 2195 × 3314 (10 m) | +| domain vertices | 740 (1 ring, 0 holes) | +| start vertices | 8242 | +| start triangles | 15614 | +| output vertices | 573758 | +| output triangles | 1145434 | +| constraint edges | 14074 | +| vertices without data dropped | 0 | +| lc_t1_run3.vtk | 48.5 MB | + +## Quality (plan view, x/y) + +| metric | median | < 1° | < 10° | worst | +|---|---|---|---|---| +| minimum angle | 45.00° | 0.10 % | 1.36 % | 0.00408° | + +| metric | median | p99 | max | ≥ 12 | ≥ 20 | +|---|---|---|---|---|---| +| vertex degree (triangles) | 6 | 9 | 29 | 1377 | 37 | + +## Refinement + +| tolerance | achieved max error | rounds | inserted | carved | flips | uncovered | quality inserted | quality skipped | feet | +|---|---|---|---|---|---|---|---|---|---| +| 1 m | 1 m | 24 | 549527 | 0 | 1152221 | 0 | 15989 | 2442 | 5726 | + +## Timings + +Wall clock, `time.perf_counter_ns` (Python) and `std::chrono::steady_clock` +(inside `refine`), one run, no warm-up. Total is the `mesh` command body, from +argument checks to the last file written; interpreter start-up and imports are +not in it. Threads: 10 (hardware concurrency). + +| phase | seconds | share | +|---|---|---| +| domain read | 0.002 | 0.1 % | +| decode | 0.313 | 8.2 % | +| features read | 0.119 | 3.1 % | +| features clip | 1.436 | 37.8 % | +| start mesh: build | 0.000 | 0.0 % | +| start mesh: node | 0.017 | 0.5 % | +| start mesh: triangulate | 0.005 | 0.1 % | +| start mesh: constraint edges | 0.017 | 0.4 % | +| refine | 0.537 | 14.1 % | +| refine: legalise start | 0.000 | 0.0 % | +| refine: start quality | 0.015 | 0.4 % | +| refine: scan (parallel) | 0.094 | 2.5 % | +| refine: split + flip (serial) | 0.378 | 10.0 % | +| refine: setup + output | 0.049 | 1.3 % | +| trim | 0.042 | 1.1 % | +| land cover | 1.263 | 33.2 % | +| write: encode | 0.015 | 0.4 % | +| write: disk | 0.027 | 0.7 % | +| other | 0.005 | 0.1 % | +| **total** | **3.800** | **100 %** | + +Sub-rows sum to their parent and are not added to the total. + +Statistics computed in 0.206 s, not included above. diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/logs/palette.out b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/palette.out new file mode 100644 index 00000000..706d3645 --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/palette.out @@ -0,0 +1 @@ +docs/benchmarks/2026-09-29/bygdin-landcover/corine_natural.json diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/logs/pmset_end.txt b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/pmset_end.txt new file mode 100644 index 00000000..a692306c --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/pmset_end.txt @@ -0,0 +1,2 @@ +Now drawing from 'AC Power' + -InternalBattery-0 (id=7929955) 80%; AC attached; not charging present: true diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/logs/pmset_start.txt b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/pmset_start.txt new file mode 100644 index 00000000..a692306c --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/pmset_start.txt @@ -0,0 +1,2 @@ +Now drawing from 'AC Power' + -InternalBattery-0 (id=7929955) 80%; AC attached; not charging present: true diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/logs/pmset_timed_end.txt b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/pmset_timed_end.txt new file mode 100644 index 00000000..a692306c --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/logs/pmset_timed_end.txt @@ -0,0 +1,2 @@ +Now drawing from 'AC Power' + -InternalBattery-0 (id=7929955) 80%; AC attached; not charging present: true diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/render.py b/docs/benchmarks/2026-09-29/bygdin-landcover/render.py new file mode 100644 index 00000000..8716e7ee --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/render.py @@ -0,0 +1,155 @@ +"""Render the 10 m Bygdin mesh in natural land-cover colours, offscreen +(@perf, 2026-09-29, increment 16c acceptance item 5). Usage, from the +repository root, .venv active (the `viewer` extra's vtk): + + python docs/benchmarks/2026-09-29/bygdin-landcover/render.py + +Writes bygdin_landcover_oblique.png (from the south-east, 2x vertical +exaggeration) and bygdin_landcover_top.png (map view). The colours come from +tin_engine.palettes.CORINE_NATURAL, the table `rasputin palette corine` writes. +Only the triangles are drawn; constraint lines (code 0) are left out, since on +a shaded 3D surface they fight the triangles for depth. +""" + +from __future__ import annotations + +import sys +from pathlib import Path + +import numpy as np +from vtkmodules.util.numpy_support import numpy_to_vtk, numpy_to_vtkIdTypeArray, vtk_to_numpy +from vtkmodules.vtkCommonCore import vtkLookupTable, vtkPoints +from vtkmodules.vtkCommonDataModel import vtkCellArray, vtkPolyData +from vtkmodules.vtkFiltersCore import vtkPolyDataNormals +from vtkmodules.vtkIOImage import vtkPNGWriter +from vtkmodules.vtkIOLegacy import vtkPolyDataReader +from vtkmodules.vtkRenderingAnnotation import vtkScalarBarActor +from vtkmodules.vtkRenderingCore import ( + vtkActor, vtkPolyDataMapper, vtkRenderer, vtkRenderWindow, vtkWindowToImageFilter, +) # fmt: skip +import vtkmodules.vtkRenderingOpenGL2 # noqa: F401 (registers the OpenGL backend) +import vtkmodules.vtkRenderingFreeType # noqa: F401 (text for the legend) + +from tin_engine.palettes import CORINE_NATURAL + +EXAGGERATION = 2.0 + + +def surface(path: Path) -> tuple[vtkPolyData, list[int]]: + """The triangles alone, centred on the origin, z exaggerated, with the + triangle part of `land_cover_code`; and the codes present.""" + r = vtkPolyDataReader() + r.SetFileName(str(path)) + r.Update() + pd = r.GetOutput() + n_lines = pd.GetNumberOfLines() + codes = vtk_to_numpy(pd.GetCellData().GetArray("land_cover_code"))[n_lines:].astype(np.int32) + pts = vtk_to_numpy(pd.GetPoints().GetData()).astype(np.float64) + pts[:, :2] -= pts[:, :2].mean(axis=0) # UTM offsets would cost float32 precision + pts[:, 2] = (pts[:, 2] - pts[:, 2].min()) * EXAGGERATION + tri = vtk_to_numpy(pd.GetPolys().GetConnectivityArray()).astype(np.int64).reshape(-1, 3) + out = vtkPolyData() + p = vtkPoints() + p.SetData(numpy_to_vtk(pts, deep=True)) + out.SetPoints(p) + cells = vtkCellArray() + cells.SetData(numpy_to_vtkIdTypeArray(np.arange(0, 3 * len(tri) + 1, 3), deep=True), + numpy_to_vtkIdTypeArray(tri.ravel(), deep=True)) # fmt: skip + out.SetPolys(cells) + arr = numpy_to_vtk(codes, deep=True) + arr.SetName("land_cover_code") + out.GetCellData().SetScalars(arr) + return out, sorted(set(codes.tolist())) + + +def lookup(present: list[int]) -> vtkLookupTable: + """Indexed mode: one annotated colour per code present.""" + lut = vtkLookupTable() + lut.IndexedLookupOn() + lut.SetNumberOfTableValues(len(present)) + for i, code in enumerate(present): + label, hexa = CORINE_NATURAL[code] + rgb = [int(hexa[j : j + 2], 16) / 255 for j in (1, 3, 5)] + lut.SetTableValue(i, *rgb, 1.0) + lut.SetAnnotation(code, f"{code} {label}") + lut.SetNanColor(1.0, 0.0, 1.0, 1.0) + return lut + + +def render(poly: vtkPolyData, lut: vtkLookupTable, out: Path, oblique: bool) -> None: + normals = vtkPolyDataNormals() + normals.SetInputData(poly) + normals.SplittingOff() + mapper = vtkPolyDataMapper() + mapper.SetInputConnection(normals.GetOutputPort()) + mapper.SetLookupTable(lut) + mapper.SetScalarModeToUseCellData() + mapper.UseLookupTableScalarRangeOn() + actor = vtkActor() + actor.SetMapper(mapper) + actor.GetProperty().SetAmbient(0.25 if oblique else 0.45) + actor.GetProperty().SetDiffuse(0.8 if oblique else 0.6) + bar = vtkScalarBarActor() + bar.SetLookupTable(lut) + bar.SetNumberOfLabels(0) + bar.SetMaximumWidthInPixels(60) + # Annotations are drawn left of the swatches; leave them room there. + bar.SetPosition(0.19, 0.52 if oblique else 0.04) + bar.SetWidth(0.21) + bar.SetHeight(0.44 if oblique else 0.40) + bar.GetAnnotationTextProperty().SetColor(0.1, 0.1, 0.1) + bar.GetAnnotationTextProperty().SetFontSize(14) + bar.GetAnnotationTextProperty().ShadowOff() + bar.SetTitle("") + ren = vtkRenderer() + ren.AddActor(actor) + ren.AddViewProp(bar) + if oblique: + ren.GradientBackgroundOn() + ren.SetBackground(0.93, 0.95, 0.97) + ren.SetBackground2(0.62, 0.74, 0.88) + else: + ren.SetBackground(1.0, 1.0, 1.0) + win = vtkRenderWindow() + win.SetOffScreenRendering(1) + win.SetSize(1400, 900) + win.SetMultiSamples(8) + win.AddRenderer(ren) + x0, x1, y0, y1, z0, z1 = poly.GetBounds() + cx, cy, cz = (x0 + x1) / 2, (y0 + y1) / 2, (z0 + z1) / 2 + cam = ren.GetActiveCamera() + cam.SetFocalPoint(cx + (0.05 * (x1 - x0) if oblique else 0), cy - (0.08 * (y1 - y0) if oblique else 0), cz) + if oblique: + # From the south-east, looking north-west, about 30 degrees down. + d = 1.45 * max(x1 - x0, y1 - y0) + cam.SetPosition(cx + 0.55 * d, cy - 0.75 * d, cz + 0.55 * d) + cam.SetViewUp(0, 0, 1) + cam.SetViewAngle(26) + else: + cam.SetFocalPoint(cx - 0.04 * (x1 - x0), cy, cz) + cam.SetPosition(cx - 0.04 * (x1 - x0), cy, cz + 10 * (y1 - y0)) + cam.SetViewUp(0, 1, 0) + cam.ParallelProjectionOn() + cam.SetParallelScale(0.55 * (y1 - y0)) + ren.ResetCameraClippingRange() + win.Render() + grab = vtkWindowToImageFilter() + grab.SetInput(win) + grab.Update() + png = vtkPNGWriter() + png.SetCompressionLevel(9) + png.SetFileName(str(out)) + png.SetInputConnection(grab.GetOutputPort()) + png.Write() + + +def main(mesh: Path, out_dir: Path) -> None: + poly, present = surface(mesh) + lut = lookup(present) + render(poly, lut, out_dir / "bygdin_landcover_oblique.png", oblique=True) + render(poly, lut, out_dir / "bygdin_landcover_top.png", oblique=False) + print(f"codes drawn: {present}") + + +if __name__ == "__main__": + main(Path(sys.argv[1]), Path(sys.argv[2])) diff --git a/docs/benchmarks/2026-09-29/bygdin-landcover/run.sh b/docs/benchmarks/2026-09-29/bygdin-landcover/run.sh new file mode 100644 index 00000000..f631c329 --- /dev/null +++ b/docs/benchmarks/2026-09-29/bygdin-landcover/run.sh @@ -0,0 +1,38 @@ +#!/bin/bash +# Increment 16c acceptance on Bygdin (@perf, 2026-09-29): the reduced +# catchment meshed with the CORINE extract, at 10 m and 1 m. Raw logs go to +# ./logs; meshes go to $SCRATCH and are not committed. +# Usage: SCRATCH=/some/dir bash run.sh (from the repository root, .venv active) +set -u +HERE=docs/benchmarks/2026-09-29/bygdin-landcover +LOG=$HERE/logs +DEM=../rasputin_data/DTM10_UTM33_20260925 +GPKG=../rasputin_data/corine2018_dtm10_utm33.gpkg +DOM=$HERE/bygdin_reduced_t20.geojson +mkdir -p "$LOG" "$SCRATCH" +: > "$LOG/exits.txt" + +mesh() { # $1 tolerance, $2 tag, $3 ascii|binary + /usr/bin/time -l rasputin mesh --dem $DEM --domain $DOM --tolerance "$1" \ + --features $GPKG --features-layer corine2018 --features-map corine \ + --"$3" --out "$SCRATCH/$2.vtk" --stats "$LOG/$2.stats.md" \ + > "$LOG/$2.out" 2> "$LOG/$2.err" + echo "$2 exit=$?" | tee -a "$LOG/exits.txt" +} + +git rev-parse HEAD > "$LOG/commit.txt" +shasum -a 256 .venv/lib/python3.*/site-packages/tin_engine/_core*.so \ + build-pyext/_core*.so >> "$LOG/commit.txt" +pmset -g batt > "$LOG/pmset_start.txt" +# Three timed binary runs at each tolerance; the median is reported. +for i in 1 2 3; do + mesh 10 "lc_t10_run$i" binary + mesh 1 "lc_t1_run$i" binary +done +pmset -g batt > "$LOG/pmset_timed_end.txt" +# One ASCII run each, for tools/bench.py's quality() (it reads ASCII only). +mesh 10 lc_t10_ascii ascii +mesh 1 lc_t1_ascii ascii +pmset -g batt > "$LOG/pmset_end.txt" +rasputin palette corine --out "$HERE/corine_natural.json" > "$LOG/palette.out" 2>&1 +echo "palette exit=$?" | tee -a "$LOG/exits.txt" diff --git a/docs/increments/16c-landcover-labels.md b/docs/increments/16c-landcover-labels.md new file mode 100644 index 00000000..f4afeb27 --- /dev/null +++ b/docs/increments/16c-landcover-labels.md @@ -0,0 +1,660 @@ +# Increment 16c — a land-cover class per triangle, in natural colours + +Status: **built (green `0487ed0`), under review.** Designed by `@architect` +2026-09-29, night. Written before +`@tester`, per `docs/increments/README.md` step 1, on branch +`increment16c-landcover-labels` off master `b4847d7` (16b-1/2 merged, 22 not). +Ola was asleep while this was written; every question a design would +otherwise put to Ola is decided below as a **Default (2026-09-29), for Ola to +confirm** (collected under "Defaults for Ola to confirm"). + +**Closes.** `ROADMAP.md`'s 16c row, 16b's Q2 (a) as ruled on 2026-09-28 +("a label per triangle by flood fill across unconstrained edges ..., written +as a cell field `land_cover` with the CLC code"), and Ola's request of +2026-09-29: *"For the example, I would like the land cover polygons colored in +a natural way for the vtk-vizualisation."* The example is increment 22's +Bygdin catchment meshed with `--features +--features-map corine`. After this increment, + +```sh +rasputin mesh --dem ../rasputin_data/DTM10_UTM33_20260925 \ + --domain bygdin_reduced_t20.geojson --tolerance 10 \ + --features ../rasputin_data/corine2018_dtm10_utm33.gpkg \ + --features-layer corine2018 --features-map corine --out bygdin.vtk +rasputin palette corine --out corine_natural.json +``` + +writes a `.vtk` whose triangles carry their CORINE code in the cell array +`land_cover_code`, and a ParaView colour preset that paints those codes in +natural colours (forest green, rock grey, heath light green, bog brown, water +blue, glacier white-blue). + +**Not closed.** "Water on exactly one side" as an edge property (16b's +Q2 note): the labels make it computable, but nothing here writes it. +Hydro-flattening lakes. A raster land-cover source (MapBiomas). Labels from a +`--features-map` without class codes (`property`; see R2). A colour table for +any code system other than CORINE. + +## Ola's direction this design is built on + +1. **16b's Q2 (a)**, ruled 2026-09-28: a separate increment, next after 16b; + a label per triangle by flood fill across unconstrained edges, the method + of Shewchuk's Triangle (`-A`), a cell field with the CLC code, about 80 + lines of pure Python over the trimmed mesh and the features. +2. **Natural colours, not the official CLC legend** (2026-09-29, via the + main session's brief): the official legend colours peat bogs blue + (`077-077-255`) and the sea almost white (`230-242-255`) + (`Legend/CLC_legend.csv` in Ola's download, read for this design). +3. **The I/O boundary** (`CLAUDE.md` §2): file decoding and CRS stay in + Python; nothing here goes near `_core`. +4. **Increment 13, ruling 5**: every cell array covers every cell; per-triangle + data gives lines a fill value. +5. **Lean** (2026-09-29): red, green and review in one night; no throwaway + implementations. + +## Prior art: legacy and literature + +### Literature + +- **Regional attributes on a constrained triangulation.** Shewchuk, + "Triangle: Engineering a 2D quality mesh generator and Delaunay + triangulator", *Applied Computational Geometry*, LNCS 1148, 203-222, 1996. + Triangle's `-A` takes a list of (seed point, attribute), finds the triangle + containing each seed, and spreads the attribute across every edge that is + not a segment; a triangle no seed reaches gets 0. Recalled (16b cited it the + same way), not reread tonight. + **What differs, and why.** R1 spreads first and looks up second: it finds + the connected components of triangles across non-constraint edges with no + seeds at all, then asks which input polygon each *component* lies in, using + one point per component. Triangle's order needs one seed per polygon piece, + chosen from the input, and that fails in three cases this project meets: + a polygon cut into two pieces by the domain's boundary (the quarter circle + does that), a polygon cut in two by a polyline running across it (a road), + and a seed that lands in the wrong triangle because the polygon is thinner + there than the noder's snap. Components-first has no seeds to place, so it + has none of the three. The spread itself (constant across non-constraint + edges, blocked by constraints, 0 where nothing applies) is Triangle's. +- **Connected components, vectorised.** Union-find with path compression: + Tarjan, "Efficiency of a good but not linear set union algorithm", *J. ACM* + 22(2):215-225, 1975. The array form R1 uses, hooking the larger root onto + the smaller and then pointer-jumping until every label is a root, is the + hook-and-shortcut scheme of Shiloach and Vishkin, "An O(log n) parallel + connectivity algorithm", *J. Algorithms* 3(1):57-67, 1982. Both recalled. + **What differs:** Shiloach-Vishkin hooks conditionally and bounds its + rounds by O(log n); R1 hooks every root to the smallest root it is joined + to and loops to a fixed point, which is correct (labels only decrease, and + hooking larger onto smaller cannot make a cycle) but claims no round bound. + The acceptance records the rounds. +- **Point in polygon** is GEOS's (shapely 2 `STRtree.query(..., + predicate="intersects")` over prepared polygons). Nothing is claimed about it + beyond what shapely documents. +- **Where the labels came from before.** The legacy labelled cell centres, + one Python point test per cell (see Legacy). The per-centre test is kept, + as the **test oracle**, not as production (R1). +- **The CORINE nomenclature.** Kosztra, Büttner, Hazeu and Arnold, *Updated + CLC illustrated nomenclature guidelines*, European Environment Agency, 2019: + the 44 level-3 classes and their codes. Recalled. The class names used in + R4 are the `LABEL3` column of `Legend/CLC_legend.csv` in Ola's download + (`../rasputin_data/corine_sql/u2018_clc2018_v2020_20u1_geoPackage/Legend/`), + read for this design. +- **Natural colour for land cover.** Patterson and Kelso, "Hal Shelton + revisited: designing and producing natural-color maps with satellite land + cover data", *Cartographic Perspectives* 47:28-55, 2004: land cover drawn in + the colours the ground has when seen from above (dark greens for conifers, + tans and greys for bare ground, blues for water, white for ice), instead of + a categorical legend. Recalled. R4's table follows that convention by hand; + it is taste, and is not claimed to reproduce their palette. +- **The ParaView preset file.** ParaView's colour-map presets are a JSON list + of objects; a categorical preset carries `"Name"`, `"Annotations"` (value, + label, value, label, ...) and `"IndexedColors"` (r, g, b in 0-1, one triple + per annotated value, in order), plus an optional `"NanColor"`. Recalled + from ParaView's own presets file, **not checked against a running ParaView + tonight** (none is installed). It is the acceptance's manual item, and the + fallback if it fails is in R4. + +**Novelty: none claimed.** Labelling a constrained triangulation's faces by +region is standard (Triangle, 1996). No search was needed because nothing new +is claimed; if a later increment claims a guarantee for region labels on +snapped input, it searches first (16b listed the queries). + +### Legacy + +```sh +$ grep -rlE 'cover_type|cover_color|def color\(' legacy/ | grep -v pycache | sort +legacy/rasputin/application.py +legacy/rasputin/globcov_repository.py +legacy/rasputin/gml_repository.py +legacy/rasputin/land_cover_repository.py +legacy/rasputin/tin_repository.py +legacy/rasputin/wfs_repository.py +legacy/tests/test_gml_repository.py +``` + +What matters, read directly: + +- `legacy/rasputin/gml_repository.py:184-226`, `land_cover`: every cell + centre tested against every polygon, one `shapely.Point` at a time (`# + TODO: Move to C++ for speed!`), raising on a centre in no polygon. + **Carried as the test oracle only** (R1, "the oracle"), with two changes: + the test is vectorised, and a centre in no polygon is 0, not an error, + because a general input need not cover the domain. +- `legacy/rasputin/gml_repository.py:69-117`, `LandCoverMetaInfo.color`: the + official CLC RGB triples. **Not carried**: Ola asked for natural colours + (direction 2). +- `legacy/rasputin/application.py:126-148`: writes the face fields + `cover_type` (the code) and `cover_color` (RGB). **Carried:** a per-face + code. **Not carried:** a per-face colour array (R4 says why) and the names + (`cover_type` collides with nothing, but R3's `land_cover_code` says what it + holds and sits next to 16b's `land_cover` bit). + +`@migration-expert` is not needed: nothing numeric is ported. + +## The blueprint: data flow and boundaries + +``` +cli.mesh (composition root) + │ + ├─ feature_input.open_features(...) 16b, extended [R5] + │ a feature from a map with class codes also keeps + │ code: int (its attribute value, e.g. Code_18 "322" -> 322) + │ polygon: Polygon | MultiPolygon | None, in the DEM's CRS + │ a coded polygon with no boundary inside the domain but covering it + │ is kept, with no lines + │ + ├─ _dem_mesh(...) -> Trimmed unchanged + │ vertices (N, 3), triangles (T, 3), edges (E, 2), edge_masks + │ + ├─ landcover.label_triangles( NEW, pure [R1, R2] + │ vertices, triangles, edges, + │ polygons=[(f.polygon, f.code) ...], margin=2 * snap spacing) + │ -> CoverLabels(codes (T,) int32, regions, outside, overlapped, thin) + │ regions(triangles, edges) -> (T,) component ids numpy only + │ one point per component -> STRtree query -> smallest-area rule + │ timed as the phase "land cover"; one stderr line + │ + ├─ io.vtk_legacy.write_vtk(..., triangle_codes=codes) [R3] + │ cell array land_cover_code: 0 on every LINES cell, the code on each + │ triangle; FieldData string land_cover_codes names the code system + ├─ io.ply.write_ply(..., face_codes=codes) [R3] + │ + └─ rasputin palette corine [--out FILE] NEW command [R4] + palettes.paraview_preset(CORINE_NATURAL) -> JSON +``` + +Boundaries: + +- `landcover.py` imports numpy and shapely, nothing first-party, never + `_core`, never a path. It takes arrays and `(geometry, code)` pairs and + returns arrays and counts. It is unit-testable on a hand-made mesh of three + triangles. +- `palettes.py` is data plus one function; it imports nothing first-party and + knows nothing about meshes. The CLI command and the acceptance's render + script both read it, so there is one colour table. +- The writers stay pure (`bytes` out) and gain one optional array each. +- No C++ changes, no bindings, no new dependency (numpy, shapely and, for + the picture only, the existing `viewer` extra's `vtk`). + +## Rulings + +### R1. The label: components first, one point per component + +**Production.** For a trimmed mesh with triangles `T`, constraint edges `E` +and coded polygons `P`: + +1. **Components** (`regions`). Each triangle's three edges as undirected keys + `min(a, b) * N + max(a, b)` (int64; N is the vertex count). A key in `E` + is blocked (`np.isin` on the same keys). The unblocked keys are sorted; + two equal neighbours in the sorted order are one interior edge and join + its two triangles. (A key appears at most twice: the trimmed mesh is a + manifold.) Components by the array union-find of the literature section: + `parent = arange(T)`; repeat { `ru, rv = parent[u], parent[v]`; + stop if `ru == rv` everywhere; `np.minimum.at(parent, max(ru, rv), + min(ru, rv))`; pointer-jump `parent = parent[parent]` until it is fixed }. + The component id of a triangle is its component's smallest triangle + index, so ids do not depend on the order edges are visited. +2. **One point per component.** Each triangle's incentre and inradius `r` + (incentre `(a·A + b·B + c·C)/(a + b + c)` with `a = |BC|` and so on, `r = + 2·area/(a + b + c)`, in x and y). Per component, the triangle with the + largest `r`, ties to the lowest triangle index; its incentre is the + component's point. +3. **Which polygon.** `STRtree([p for p, _ in P]).query(points, + predicate="intersects")`. A point in no polygon gives its component 0. A + point in several gives the code of the one with the **smallest area**, ties + to the smaller code (Default D2). +4. **Spread.** Every triangle gets its component's code. + +**Why one point is enough, and when it is not.** A component is bounded by +constraint edges, and the input polygon boundaries lie within `δ` of the +constraint edges made from them. `05-noder.md` (fix 8, the one-way Hausdorff +form) bounds the other direction: every output edge derived from an input +segment lies within `h/√2` of it per noding round, `k·h/√2` after `k` rounds, +for snap spacing `h`. The output chain runs, inside that tube, from near one +end of the input segment to near the other, so it crosses the perpendicular +through any point of the segment within the same distance: `δ <= k·h/√2` +(0.7 mm per round at the default `h = 1 mm`). `margin = 2h` covers two +rounds. Refinement's points on constraint segments (20b's feet) are rounded +to doubles, which adds nanometres. This bound says when the lookup is exact; +the tests rest on the oracle, not on it. +The incircle of a triangle lies inside the triangle, and no constraint edge +enters a triangle, so the incentre is at least `r` from every constraint edge +and at least `r - δ` from every input boundary. With `margin = 2h > δ`, a +component whose largest `r` exceeds `margin` has a point that is strictly on +one side of every input boundary, so the lookup is exact up to GEOS's own +arithmetic. A component whose largest `r` is at most `margin` is counted as +`thin`: it is still labelled by its point, but the stderr line says how many +there were. (If the noder took more than two rounds, `δ` can exceed +`margin`; the count is then optimistic, and the oracle is the check.) On CORINE such a component is a sliver the snap made between two +boundaries that nearly coincide. + +**Is the polygon the same object the constraints came from?** Yes, by R5: the +polygon tested is the feature's geometry moved into the DEM's CRS vertex by +vertex with the same transform that moved its chains, so its edges are the +straight lines the engine draws. The domain clip is linework on the chains +and is not repeated on the polygon; the triangles only exist inside the +domain anyway. + +**The oracle (tests only).** Each triangle's **centroid** is tested against +`P` with the same predicate and the same overlap rule, for every triangle with +`r > 1.5·margin`. The centroid is at least `2r/3` from the triangle's sides +(its distance to side `a` is `2·area/(3a)`, and `a <= perimeter/2`), so for +those triangles it is further than `margin > δ` from any input boundary and +the oracle's answer is exact. Production and oracle must agree on every such +triangle. This is the legacy's per-centre test; it borrows the producer's +*predicate* (a point against the input polygons) and none of its records (no +components, no chosen points), per the `computational-geometry` skill's rule. +It is O(T) point queries, which is why it is the oracle and not production. + +**Degeneracy and edge cases.** + +| case | outcome | +|---|---| +| overlapping polygons (not CORINE; any GeoJSON) | smallest area wins, ties to the smaller code; the component is counted in `overlapped` (D2) | +| two polygons sharing a boundary that the noder merged (every CORINE neighbour pair) | one constraint edge between them; its two sides are different components, each labelled by its own point | +| boundaries closer than the snap, merged or not | at most a sliver component, labelled by its point, counted `thin` | +| a polygon with holes | shapely's test excludes the hole; the hole's components get whatever polygon fills it, else 0 | +| a lake inside a forest, the forest with a hole (CORINE) | lake components 512, forest components 31x | +| a lake inside a forest, the forest without a hole (hand-made input) | lake components the lake's code, by the smallest-area rule | +| a polygon clipped by the domain into several pieces | each piece is its own component and is labelled | +| a polyline across a polygon (a road through a forest) | two components, both the forest's code | +| a polygon covering the whole domain, no boundary inside | kept with no lines (R5); the single component gets its code | +| a triangle outside the domain | none exist: the engine's mesh is the domain's | +| a triangle dropped by `trim` (no DEM data) | not in the file, not labelled; components are computed after trimming | +| a component in no polygon | 0, counted in `outside` | +| no coded polygons at all (`property` map, or no `--features`) | no labelling, no array, no field (R2) | +| a constraint edge with mask 0 (the domain boundary, an unclassified edge) | blocks the spread like any constraint (it is in `E`) | + +**Determinism.** The codes are a function of the triangle array, the +constraint edges and the polygons as a set: component ids are +smallest-index, the chosen point is largest-`r` then smallest-index, and the +overlap rule is order-free. 16b's I3 (rows in reverse order give a +bit-identical file) therefore extends to the new array with no new code; its +existing test covers it. + +### R2. Where it runs: Python, over the trimmed mesh's arrays + +In `src_python/tin_engine/landcover.py`, on numpy arrays, after `trim` and +before the writers. Not in the C++ core, because: + +- the class codes and polygons are attribute data read in Python, in a CRS + handled in Python; the core would need neither the polygons nor the codes, + only the component pass, and that is about 25 lines of numpy; +- a C++ pass would need bindings, a rebuild, a C++ suite and, touching + `include/terrain/mesh/`, `@perf`'s bench acceptance, which is not what a + night is for. + +**Cost at scale.** One sort of `3T` int64 keys, one `np.isin`, the union-find +rounds (each O(T)), and one point query per component, not per triangle. At +1 m on Bygdin (1.13 M triangles without features) that is a few million keys; +**not measured** (no prototype, per the retrospective agenda of 2026-09-29). +The acceptance records the phase time and the round count at 10 m and at 1 m. +If the phase costs more than the refine it follows at 1 m, moving `regions` +into the core is the next step, and the Python version stays as its oracle. + +**Which maps label.** Only a class map that declares a code system does +(`ClassMap.codes`, R5): `corine` and `clc18_kode` (and `corine-water`, whose +kept features are the water polygons only, so land is 0). The `property` map's +values are vocabulary names, not codes; it writes no labels (Default D3). + +### R3. What the files carry + +**`.vtk`** (`io/vtk_legacy.py`, `write_vtk(..., triangle_codes=None)`): + +- A cell array **`land_cover_code`**, `int`, one value per cell in the file's + cell order: **0 on every `LINES` cell**, then each triangle's code (13's + rulings 4 and 5: lines first, every array covers every cell, lines get the + fill value). It is written in the `FIELD` block beside the per-feature + arrays, not as the `SCALARS` block, which stays `feature_mask` (13, ruling + 5 is not reopened). The `FIELD` count includes it; with no feature arrays + the block is written for it alone. +- **Not named `land_cover`**: 16b's vocabulary bit 7 is `land_cover`, and + ruling 6 already writes a 0/1 cell array under that name. The writer refuses + `triangle_codes` with a vocabulary that names a property `land_cover_code`. +- A dataset string **`land_cover_codes`** in `FieldData`: the code system, + the attribute, the map, and what 0 means, e.g. + `CORINE Land Cover level-3 code, attribute Code_18, map corine; 0 = in no + polygon, and every constraint line`. `land_cover_codes` joins `RESERVED`. +- `triangle_codes` must have one entry per triangle and fit `int32`; else + `ValueError`. + +**`.ply`** (`io/ply.py`, `write_ply(..., face_codes=None)`): on the face file +only, a face property `int land_cover_code` after `vertex_indices`, and the +same text as a header comment `land_cover_codes ...`. Refused with `edges` +(the edge file has no faces). MDAL is expected to offer it as a face dataset +in QGIS; recalled, not measured, and not an acceptance item. + +**stderr**, one line after the `trim` line: +`land cover: regions, outside every polygon, in more than one, + thinner than the snap`. The phase `land cover` appears in `--stats`'s +phase table through the existing clock; no other report row (16b's ruling +that a table in a run report is output, not data). + +### R4. Natural colours: a ParaView preset, generated from one table + +**Where the table lives.** `src_python/tin_engine/palettes.py`: +`CORINE_NATURAL: dict[int, tuple[str, str]]`, code → (CLC `LABEL3`, `#rrggbb`), +all 44 CLC level-3 classes plus 0, and `paraview_preset(table, name) -> +list[dict[str, object]]`. It ships in the package, is type-checked, and the +tests pin it. Not a JSON data file in the package: `wheel.packages` would +carry it, but a Python table is one object that the CLI and the render +script both import, with no file-location lookup. + +**How a user gets it.** A new command, `rasputin palette corine [--out +FILE]`: the ParaView preset as JSON, to `FILE` or stdout. An unknown name is a +usage error listing the known ones. `mesh` does not write it beside the mesh +(Default D4): one more file per run, and a second name to keep from +colliding with `--out-edges` and `--stats`, for a table that never changes +between runs. + +**How it is applied in ParaView** (the manual acceptance item): + +1. Open `bygdin.vtk`; *Color By* `land_cover_code` (cell data). +2. Colour Map Editor → *Choose Preset* → *Import*, pick + `corine_natural.json`, select "rasputin CORINE natural", *Apply*. Imported + presets persist in ParaView's settings, so this is once per machine. +3. If the categorical mode is not switched on by the preset, tick + *Interpret Values As Categories*. +4. Optionally *Save current colour map as default for arrays named + `land_cover_code`* so every later file opens coloured. + +Constraint lines carry 0 and draw in the colour for 0, a dark grey, which +reads as a class boundary. A code the table lacks draws in `NanColor`, +magenta, so it is seen rather than blended in. + +**Rejected: a per-cell RGB array in the file** (the legacy's `cover_color`). +It needs no import, but it bakes one taste into every mesh file, adds three +bytes per cell, and a second viewer would want its own. It stays the fallback +if the preset fails the manual check. + +**The table.** Natural colours, chosen by hand after Patterson and Kelso; the +official CLC colour is not used for any class. Hex is the value `@tester` pins +by family (Tests), not by exact triple. + +| code | class (`LABEL3`) | colour | in the Norway extract | +|---|---|---|---| +| 0 | no polygon; constraint lines | `#3c3c3c` | (lines) | +| 111 | Continuous urban fabric | `#8c5f5a` | yes | +| 112 | Discontinuous urban fabric | `#a88d86` | yes | +| 121 | Industrial or commercial units | `#8e8a96` | yes | +| 122 | Road and rail networks and associated land | `#6e6e6e` | yes | +| 123 | Port areas | `#7d8796` | yes | +| 124 | Airports | `#a8a8a8` | yes | +| 131 | Mineral extraction sites | `#b49b78` | yes | +| 132 | Dump sites | `#8a7d64` | yes | +| 133 | Construction sites | `#bcae98` | yes | +| 141 | Green urban areas | `#86b86e` | yes | +| 142 | Sport and leisure facilities | `#a3cf7e` | yes | +| 211 | Non-irrigated arable land | `#e6d58c` | yes | +| 212 | Permanently irrigated land | `#d9cc6e` | | +| 213 | Rice fields | `#cdd89a` | | +| 221 | Vineyards | `#9c7a44` | | +| 222 | Fruit trees and berry plantations | `#a9ad5e` | yes | +| 223 | Olive groves | `#8f9a52` | | +| 231 | Pastures | `#b4d47a` | yes | +| 241 | Annual crops associated with permanent crops | `#dccf94` | | +| 242 | Complex cultivation patterns | `#d2c47c` | yes | +| 243 | Land principally occupied by agriculture, with significant areas of natural vegetation | `#bcc47e` | yes | +| 244 | Agro-forestry areas | `#a9b574` | | +| 311 | Broad-leaved forest | `#4f8f3f` | yes | +| 312 | Coniferous forest | `#1f5a2e` | yes | +| 313 | Mixed forest | `#357438` | yes | +| 321 | Natural grasslands | `#c2d68a` | yes | +| 322 | Moors and heathland | `#a7c47f` | yes | +| 323 | Sclerophyllous vegetation | `#8a9658` | | +| 324 | Transitional woodland-shrub | `#7ea65a` | yes | +| 331 | Beaches, dunes, sands | `#e9ddb2` | yes | +| 332 | Bare rocks | `#8f8f8f` | yes | +| 333 | Sparsely vegetated areas | `#b9b8a0` | yes | +| 334 | Burnt areas | `#4b3f3a` | yes | +| 335 | Glaciers and perpetual snow | `#eef6fb` | yes | +| 411 | Inland marshes | `#6e9470` | yes | +| 412 | Peat bogs | `#8a6642` | yes | +| 421 | Salt marshes | `#7f9f8c` | | +| 422 | Salines | `#d8d6cc` | | +| 423 | Intertidal flats | `#b3bdb3` | yes | +| 511 | Water courses | `#4c8ec4` | yes | +| 512 | Water bodies | `#3e7bb6` | yes | +| 521 | Coastal lagoons | `#5b93b3` | | +| 522 | Estuaries | `#5188b4` | yes | +| 523 | Sea and ocean | `#2b5d8e` | yes | + +The "in the Norway extract" column is from +`sqlite3 -readonly ../rasputin_data/corine2018_dtm10_utm33.gpkg "select +code_18, count(*) from corine2018 group by code_18"`, run for this design: 34 +codes, all present in the table. The unclassified codes 990, 995 and 999 do +not occur there and are left out; they draw in `NanColor`. + +### R5. What `feature_input` keeps for labelling + +- `ClassMap` gains `codes: str = ""`, the name of the code system. Non-empty + means the map's attribute values are integer class codes and labels are + written. `corine`, `corine-water` and `clc18_kode` set it to `CORINE Land + Cover level-3 code`; `property` leaves it empty. +- `TerrainFeature` gains `code: int | None = None` and `polygon: Polygon | + MultiPolygon | None = None`, both set only for a coded map: `code` is + `int(str(value))` of the feature's attribute, refused (`FeatureError`, + naming the feature and value) unless it is an integer in 1 .. 2³¹-1 — a + list value is refused too; `polygon` is the feature's polygonal parts, + moved into the DEM's CRS with the **same** transform that moved its chains + (on the moved-first path it is `moved` itself; on 16b's geographic + pre-clip path the whole polygon is moved once more). Lines have no polygon. +- **A coded polygon with no boundary in the domain** is today counted + `outside` and dropped. For a coded map it is kept with `lines=()` if it + covers the domain, tested by one point, `domain.polygon.point_on_surface()`: + no boundary crosses the domain's interior, so the interior is wholly inside + the polygon or wholly outside it. Without this, a catchment lying entirely + inside one large CORINE polygon would be labelled 0 everywhere. It counts as + kept, not clipped. +- Uncoded maps behave exactly as in 16b: no polygon is moved or kept, so no + cost is added to them. +- Codes are checked only on the features the map keeps: under + `corine-water`, an unlisted value is dropped before its code is read, as + 16b drops any unlisted value. + +### R6. The CLI + +- `mesh`: labelling runs whenever `--features` is given with a coded map and + the output is `.vtk` or `.ply`. No new flag (Default D1). +- `palette`: R4. +- The `mesh` help text gains one sentence naming `land_cover_code` and + pointing at `rasputin palette`. + +## Invariants + +- **I1 (spread).** For every interior edge that is not a constraint edge, its + two triangles have the same code. Checked on the arrays alone, no geometry. +- **I2 (oracle).** For every triangle with `r > 1.5·margin`, the code equals + the centroid oracle's (R1). +- **I3 (file shape).** `land_cover_code` has one value per cell, 0 on every + `LINES` cell, and the triangle part equals `CoverLabels.codes`. +- **I4 (determinism).** 16b's reversed-rows run gives a bit-identical file. +- **I5 (area).** On a fixture whose polygons are known, the triangle area + (x, y) per code equals the area of `polygon ∩ domain` per code within + `margin × (boundary length inside the domain)` plus rounding. + +## Defaults for Ola to confirm + +Each marked **Default (2026-09-29), for Ola to confirm**. + +- **D1.** Labels are always written for a coded map; no `--no-labels` flag. +- **D2.** Overlapping polygons: the smallest area wins, ties to the smaller + code. The alternative, the last feature in source order wins (painter's + rule), makes the result depend on row order, which 16b's I3 forbids. +- **D3.** The `property` map writes no labels. Giving its names integer ids + would need an id table in the file; nobody has asked for it. +- **D4.** `rasputin palette corine` writes the preset on request; `mesh` + does not write it beside the mesh. +- **D5.** Names: the cell array `land_cover_code`, the field + `land_cover_codes`, the command `palette`, the preset "rasputin CORINE + natural". 16b's Q2 said "a cell field `land_cover`"; that name is taken by + bit 7's 0/1 array. +- **D6.** The colours in R4's table. Taste, and Ola's to change; the tests + pin only the families (green forest, grey rock, light-green heath, brown + bog, blue water, white-blue glacier). +- **D7.** No suite is named invariant-critical, so no mutation round: every + labelling fixture is checked against the independent oracle (I2) and I1, + which is what a mutation round would test (Ola's lean-brief rule). + +## Tests for `@tester` (the red suite) + +Lean: no throwaway implementation, no mutation round (D7). Two new files and +additions to three. + +**`tests/python/test_landcover.py`** — pure, no `_core`, hand-made meshes: + +- `regions`: two triangles sharing an unconstrained edge are one component; + the same with the edge in `E` are two; a strip of 1 000 triangles with no + constraints is one component (the loop reaches a fixed point); ids are the + smallest triangle index; permuting `E`'s rows changes nothing. +- `label_triangles`: a square cut by a constrained diagonal with a polygon + on each side gives each side its code; no polygons gives all 0 and + `outside` = regions; two nested polygons give the inner one's code + (smallest area) and count `overlapped`; a sliver component (largest `r` <= + margin) is counted `thin`; `codes` is int32 of length T. +- The oracle as a test helper (`landcover_oracle(vertices, triangles, + polygons, margin)`, per R1), used by the CLI tests below. + +**`tests/python/test_cli_mesh_landcover.py`** — through `rasputin mesh` on the +synthetic `bumpy` DEM of `test_cli_mesh_features.py`, GeoJSON features with a +`Code_18` attribute and `--features-map corine`. For each fixture: I1, I2, +I3, and the expected codes named below. + +1. **Two squares side by side** (311 and 512) sharing an edge: the triangles + on either side of the shared edge carry 311 and 512; nothing is 0 inside + the squares, the rest of the domain is 0. +2. **A hole**: a forest (312) with a square hole and nothing in it: the hole + is 0. +3. **A lake inside a forest**: (a) the forest holed and the lake filling the + hole: 512 and 312; (b) the forest not holed: still 512 in the lake + (D2), and the stderr line says `in more than one` > 0. +4. **A polygon clipped by the domain into two pieces** (a band crossing a + concave notch of the domain): both pieces carry its code. +5. **A road across a polygon** (a `LineString` feature with a code, under + `corine`): both sides carry the polygon's code. +6. **A polygon covering the domain** with no boundary inside: every triangle + carries its code; the feature is counted kept. +7. **Boundaries 0.1 mm apart** (two squares whose shared side is offset + by less than the snap): both squares labelled correctly away from the + seam (I2), and the run succeeds. +8. **Refusals**: a `Code_18` of `"forest"` under `corine` is a usage error + naming the feature and value; `--features-map property` writes no + `land_cover_code` and no `land_cover_codes`. +9. **`.ply`**: the face file has `property int land_cover_code` and the + comment; the edge file has neither. +10. **VTK readback** (`importorskip("vtk")`, the `viewer` extra's CI step): + `vtkPolyDataReader` sees `land_cover_code` with `GetNumberOfCells()` + values, 0 on the line cells. + +**Additions to `test_cli_mesh_features.py::TestCommittedExtract`** (the +quarter circle on the committed tile and extract, at 10 m): I1, I2, I3 on the +real mesh; I5 per code against shapely's `intersection` of each extract +polygon (moved to EPSG:25833) with the quarter circle, within 0.01 % of the +domain's area; the reversed-rows test (I4) passes unchanged. + +**Additions to `test_io_vtk_legacy.py` and `test_io_ply.py`**: the array's +place and fill in ASCII and binary, the length refusal, the refusal of a +vocabulary naming `land_cover_code`, and the reserved field name. + +**`tests/python/test_palettes.py`**: every one of the 34 Norway codes (R4's +list) and 0 is in `CORINE_NATURAL`; the 44 CLC codes are all there; families: +311-313 have g > r and g > b; 322 has g > r and g > b and is lighter than +312; 332 is grey (r = g = b); 412 has r > g > b (brown); every 5xx has b > r +and b > g; 335 has every channel >= 0xe0 and b >= r; `paraview_preset` +gives one object with `Name`, `Annotations` of length 2 × (entries), +`IndexedColors` of length 3 × (entries) in [0, 1], in the same order; +`rasputin palette corine` writes that JSON to stdout and to `--out`; an +unknown name is a usage error. + +## Acceptance (`@perf`, AC power recorded) + +On this branch, with the Bygdin reduced catchment taken from 22's branch +(no merge needed; 16c does not depend on 22's code): + +```sh +git show increment22-autocatchment:docs/benchmarks/2026-09-29/bygdin/bygdin_reduced_t20.geojson \ + > $SCRATCH/bygdin_reduced_t20.geojson +rasputin mesh --dem ../rasputin_data/DTM10_UTM33_20260925 \ + --domain $SCRATCH/bygdin_reduced_t20.geojson --tolerance 10 \ + --features ../rasputin_data/corine2018_dtm10_utm33.gpkg \ + --features-layer corine2018 --features-map corine \ + --out $SCRATCH/bygdin_lc.vtk --stats bygdin_lc.stats.md +rasputin palette corine --out corine_natural.json +``` + +Evidence under `docs/benchmarks/2026-09-29/16c-bygdin/`: + +1. **It opens in VTK with the array**: `vtkPolyDataReader` reads the file; + `land_cover_code` is present with one value per cell and 0 on every line. +2. **Class shares against CORINE clipped to the catchment.** Triangle area + (x, y) per code, as shares of the mesh's area, against 22's table + (`bygdin/README.md` on 22's branch: 333 42.34 %, 332 21.32 %, 322 + 17.00 %, 512 16.46 %, 335 2.44 %, 412 0.33 %, 142 0.10 %). **Pass:** the + same seven codes, no 0 triangles, and every share within 0.01 percentage + points. The expected gap is `δ` times the boundary length, well under a + hectare; anything larger is a leak between regions. +3. **I1 and I2** on the Bygdin mesh (the oracle over all ~83 000 triangles). +4. **Time**: the `land cover` phase and the union-find rounds at 10 m, and + one run at `--tolerance 1` with the same features for the scale. Recorded, + not gated (R2). +5. **A picture**, `bygdin_landcover.png`: `render.py` in the evidence + directory, using the `viewer` extra's `vtk` offscreen (`vtkRenderWindow` + with `SetOffScreenRendering(1)`, `vtkWindowToImageFilter`, `vtkPNGWriter`; + this design checked that the chain writes a PNG with vtk 9.7.0 in the + repository's venv on this Mac), a map view from above with a + `vtkLookupTable` in indexed mode built from `palettes.CORINE_NATURAL`, and + a legend of the codes present. No new dependency. +6. **Manual, for Ola**: import the preset in ParaView and record the version + and whether it switched to categorical by itself (R4, steps 2-3). + +`@perf`'s bench and thread sweep are not required: nothing under +`include/terrain/refinement/` or `include/terrain/mesh/` or what drives them +changes; labelling runs after the mesh is finished. + +## LOC + +Production lines as `CLAUDE.md` §2 counts them (tests excluded): + +| file | estimate | as built (green `0487ed0`) | +|---|---|---| +| `landcover.py` (`regions`, incentres, `label_triangles`, `CoverLabels`) | 75 | 83 | +| `palettes.py` (45-entry table, `paraview_preset`) | 60 | 64 | +| `feature_input.py` (`codes`, `code`, `polygon`, coded refusals, covering polygon) | 30 | 22 | +| `io/vtk_legacy.py` (array, fill, field, refusals) | 15 | 19 | +| `io/ply.py` (face property, comment, refusal) | 12 | 16 | +| `cli.py` (label call and phase, stderr line, fields and comments, `palette` command, help) | 40 | 45 | +| **total** | **~230** | **net 249 (262 added, 13 removed)** | + +The estimate's worst case was about 320. Well under the 700-line ceiling; one +PR. + +**As-built writer API.** `write_vtk(..., triangle_codes=None, +land_cover_codes="")`: the codes array and the `land_cover_codes` field text +are separate arguments, the text written only with the array. +`write_ply(..., face_codes=None)`. + +## Not in scope + +- An edge property "water on exactly one side" (computable from the labels). +- Labels for the `property` map, or for any map without integer codes (D3). +- Colour tables other than CORINE's; a colour array in the mesh file (R4). +- Moving `regions` into the C++ core (R2 says when it would be). +- Changing which array is the `.vtk` file's active `SCALARS` (13, ruling 5). diff --git a/project_structure.md b/project_structure.md index 7f5e37ff..60f2355d 100644 --- a/project_structure.md +++ b/project_structure.md @@ -79,7 +79,9 @@ bindings/ src_python/tin_engine/ # public Python API (distribution name: rasputin) __init__.py # re-exports from tin_engine._core - cli.py # Typer entry point declared in pyproject + cli.py # Typer entry point declared in pyproject; + # `rasputin palette NAME [--out FILE]` writes a + # ParaView colour preset (16c) raster.py # the ONLY adapter from decoded data into _core grid_domain.py # DEM extent -> stride-subsampled nodes + outer ring; # pure numpy, never imports _core @@ -103,6 +105,13 @@ src_python/tin_engine/ # public Python API (distribution name: rasputin) # moved to the DEM's CRS and clipped to the domain # as linework -> FeatureSet (16b); opens GeoJSON # and .gml itself; never imports _core + landcover.py # regions, label_triangles: a land-cover code per + # triangle, components across unconstrained edges, + # one point-in-polygon test per component (16c); + # numpy and shapely, never imports _core + palettes.py # CORINE_NATURAL (code -> label, colour) and + # paraview_preset(); data, imports nothing + # first-party (16c) chains.py # start_chains: the domain's rings, then every # feature line, as (indices, role, mask) (16b); # never imports _core @@ -117,11 +126,18 @@ src_python/tin_engine/ # public Python API (distribution name: rasputin) # carries positions rather than feature names svg.py # (Scene, SvgStyle) -> str; the stylesheet lives here fixtures.py # the synthetic gallery, declarative; `rasputin draw - io/ # all file decoding AND encoding lives here + io/ # all file decoding AND encoding lives here, with + # known exceptions to move here in a follow-up: + # `palette`'s JSON is encoded in cli.py (16c), and + # increment 22's GeoJSON writer (on 22's branch) __init__.py - ply.py # arrays -> PLY bytes; takes no path and opens nothing + ply.py # arrays -> PLY bytes; takes no path and opens nothing; + # face_codes= adds the face property + # land_cover_code (16c) vtk_legacy.py # arrays + EdgeVocabulary -> legacy .vtk bytes, for - # ParaView; takes no path and opens nothing + # ParaView; takes no path and opens nothing; + # triangle_codes=, land_cover_codes= add the cell + # array land_cover_code and its field (16c) geotiff.py # TIFF container + GeoKey decoding -> DemTile models.py # Pydantic RasterMeta / DemTile geopackage.py # GeoPackage layer_info / query_features over an diff --git a/src_python/tin_engine/cli.py b/src_python/tin_engine/cli.py index 0dbb0d38..70129b01 100644 --- a/src_python/tin_engine/cli.py +++ b/src_python/tin_engine/cli.py @@ -40,6 +40,7 @@ from __future__ import annotations import importlib.metadata +import json import math import os import shlex @@ -75,6 +76,7 @@ from tin_engine.elevation import Trimmed, trim from tin_engine.feature_input import ( CLASS_MAPS, + ClassMap, FeatureError, FeatureRequest, FeatureSet, @@ -86,7 +88,9 @@ from tin_engine.io.models import DemTile, RasterMeta from tin_engine.io.ply import write_ply from tin_engine.io.vtk_legacy import write_vtk +from tin_engine.landcover import label_triangles from tin_engine.mosaic import Bounds, Seam +from tin_engine.palettes import PALETTES, paraview_preset from tin_engine.raster import to_core from tin_engine.stats import PhaseClock, Refinement, Report, Sizes, _exact, quality, render from tin_engine.viz.fixtures import GALLERY, Fixture @@ -696,9 +700,16 @@ def mesh( ``--stats`` (increment 17) adds a report and changes nothing else: the clock always runs, the quality pass and the report only with the flag. + + A ``--features-map`` with class codes (``corine``, ``corine-water``, + ``clc18_kode``) also gives every triangle its polygon's code, in the cell + array ``land_cover_code`` of a ``.vtk`` or the face property of a ``.ply`` + (increment 16c); ``rasputin palette corine`` writes natural colours for it. """ clock = PhaseClock() dem_run: _DemMesh | None = None + codes: npt.NDArray[np.int32] | None = None + codes_text = "" seams: tuple[Seam, ...] = () if (name is None) == (not dem): raise typer.BadParameter( @@ -839,6 +850,9 @@ def mesh( if notice: fields.append(("features_notice", notice)) comments.append(f"features_notice {notice}") + cmap = CLASS_MAPS[features_map or "property"] + if cmap.codes: + codes, codes_text = _land_cover(dem_run.trimmed, found, cmap, snap_spacing, clock) else: assert name is not None if stride is not None: @@ -898,12 +912,18 @@ def mesh( vocabulary=DEFAULT_VOCABULARY, fields=fields, binary=binary, + triangle_codes=codes, + land_cover_codes=codes_text, ) ] else: encoders = [ lambda: write_ply( - vertices, faces=surface_mesh.triangles, ascii=not binary, comments=comments + vertices, + faces=surface_mesh.triangles, + ascii=not binary, + comments=[*comments, *([f"land_cover_codes {codes_text}"] if codes_text else [])], + face_codes=codes, ), lambda: write_ply( vertices, @@ -924,6 +944,48 @@ def mesh( _write_report(clock, report_target, surface_mesh, dem_run, targets, seams) +def _land_cover( + trimmed: Trimmed, found: FeatureSet, cmap: ClassMap, spacing: float, clock: PhaseClock +) -> tuple[npt.NDArray[np.int32], str]: + """16c, R1-R3: a code per triangle from the coded polygons, timed as the + phase ``land cover``, with its stderr line; the codes and their text.""" + polygons = [(f.polygon, f.code) for f in found.features if f.polygon is not None and f.code] + with clock.phase("land cover"): + mesh = (trimmed.vertices, trimmed.triangles, trimmed.edges) + labels = label_triangles(*mesh, polygons=polygons, margin=2 * spacing) + typer.echo( + f"land cover: {labels.regions} regions, {labels.outside} outside every polygon, " + f"{labels.overlapped} in more than one, {labels.thin} thinner than the snap", + err=True, + ) + text = f"{cmap.codes}, attribute {cmap.attribute}, map {cmap.name}; " + return labels.codes, text + "0 = in no polygon, and every constraint line" + + +@app.command() +def palette( + name: Annotated[str, typer.Argument(help=f"The palette: {', '.join(PALETTES)}.")], + out: Annotated[ + Path | None, typer.Option("--out", help="Write the preset here, not to stdout.") + ] = None, +) -> None: + """Write a ParaView colour preset for ``land_cover_code`` (increment 16c). + + Import it in ParaView's Colour Map Editor (Choose Preset, Import), then + colour by ``land_cover_code``; tick Interpret Values As Categories if the + preset does not. + """ + if name not in PALETTES: + raise typer.BadParameter(f"unknown palette {name}; use {', '.join(PALETTES)}") + table, title = PALETTES[name] + text = json.dumps(paraview_preset(table, title), indent=1) + "\n" + if out is None: + typer.echo(text, nl=False) + return + out.write_text(text) + typer.echo(f"{out}") + + def _report_target( stats: str | None, out_parent: Path | None, label: str, meshes: list[Path] ) -> Path | None: diff --git a/src_python/tin_engine/feature_input.py b/src_python/tin_engine/feature_input.py index c72cc34f..54a21beb 100644 --- a/src_python/tin_engine/feature_input.py +++ b/src_python/tin_engine/feature_input.py @@ -36,7 +36,7 @@ import shapely from pydantic import BaseModel, ConfigDict from pyproj import CRS -from shapely.geometry import LineString, Polygon, shape +from shapely.geometry import LineString, MultiPolygon, Polygon, shape from shapely.geometry.base import BaseGeometry from tin_engine.crs import parse_crs, reprojector, transform_definition @@ -55,6 +55,8 @@ #: WGS 84's smallest and largest radii of curvature, metres (R5, "Long edges"). R_MIN, R_MAX = 6_335_439.0, 6_399_594.0 WATER_CODES = ("511", "512", "521", "522", "523") +#: A map whose values are integer class codes names its code system (16c, R5). +CORINE_CODES = "CORINE Land Cover level-3 code" CORINE_NOTICE = ( "Contains modified CORINE Land Cover 2018 data (version 2020_20u1), (c) European Union, " "Copernicus Land Monitoring Service 2018, European Environment Agency (EEA): clipped and " @@ -70,7 +72,9 @@ class FeatureError(ValueError): class ClassMap(BaseModel): """An attribute's values to vocabulary names (R4). ``otherwise`` is what an - unlisted value gets: names, ``"drop"`` (no constraint) or ``"refuse"``.""" + unlisted value gets: names, ``"drop"`` (no constraint) or ``"refuse"``. + ``codes``, when not empty, names the code system the values are integer + codes of; such a map labels triangles (16c, R5).""" model_config = ConfigDict(frozen=True, extra="forbid") @@ -79,6 +83,7 @@ class ClassMap(BaseModel): classes: Mapping[str, tuple[str, ...]] otherwise: Literal["refuse", "drop"] | tuple[str, ...] = "refuse" notice: str = "" + codes: str = "" def _corine(name: str, attribute: str) -> ClassMap: @@ -89,6 +94,7 @@ def _corine(name: str, attribute: str) -> ClassMap: classes=water, otherwise=("land_cover",), notice=CORINE_NOTICE, + codes=CORINE_CODES, ) @@ -105,6 +111,7 @@ def _corine(name: str, attribute: str) -> ClassMap: classes=dict.fromkeys(WATER_CODES, ("water",)), otherwise="drop", notice=CORINE_NOTICE, + codes=CORINE_CODES, ), "clc18_kode": _corine("clc18_kode", "clc18_kode"), } @@ -128,13 +135,17 @@ class FeatureRequest(BaseModel): class TerrainFeature(BaseModel): """One feature: its fid, mask, and clipped lines in the DEM's CRS; a - closed line (first point repeated) is an unclipped ring.""" + closed line (first point repeated) is an unclipped ring. Under a coded + map (16c, R5) it also keeps its class ``code`` and, if polygonal, its + ``polygon`` moved into the DEM's CRS by the transform its lines took.""" model_config = ConfigDict(frozen=True, arbitrary_types_allowed=True) fid: int | str mask: int lines: tuple[LineString, ...] + code: int | None = None + polygon: Polygon | MultiPolygon | None = None class FeatureSet(BaseModel): @@ -333,18 +344,23 @@ def _take(self, source: FeatureSource, own: str, rows: Iterable[tuple[Any, Any, mask = _mask(cmap, self.vocabulary, value, f"{name}: feature {fid}") if mask is None: continue - moved = None + code = _code(value, f"{name}: feature {fid}") if cmap.codes else None + polygonal = code is not None and geometry.geom_type in ("Polygon", "MultiPolygon") + moved = polygon = None if move is not None and bound is not None and bound[1](geometry): # Pre-clipped where its edges are straight, then moved (R5). chains = pre_clip(geometry, bound[2], bound[0]) kept = [shapely.transform(c, move) for c in chains] + polygon = shapely.transform(geometry, move) if polygonal else None else: # moved first, then pre-clipped exactly in the DEM's CRS moved = shapely.transform(geometry, move) if move is not None else geometry kept = [moved] + polygon = moved if polygonal else None if not all(_finite(g) for g in kept): raise FeatureError(f"{name}: feature {fid} has a vertex with no image in the DEM") if moved is not None: kept = list(pre_clip(moved, region)) + coded: dict[str, Any] = {"code": code, "polygon": polygon} lines = tuple( piece for line in kept @@ -352,9 +368,13 @@ def _take(self, source: FeatureSource, own: str, rows: Iterable[tuple[Any, Any, if isinstance(piece, LineString) and piece.length > 0 ) if lines: - self.features.append(TerrainFeature(fid=fid, mask=mask, lines=lines)) + self.features.append(TerrainFeature(fid=fid, mask=mask, lines=lines, **coded)) # A dropped edge lies outside, so a pre-clipped chain is not covered. self.clipped += not all(self.domain.polygon.covers(g) for g in kept) + elif polygon is not None and polygon.intersects(self.domain.polygon.point_on_surface()): + # No boundary crosses the domain, and a point of it is inside: it + # covers the domain, and labels it (R5). Kept, not clipped. + self.features.append(TerrainFeature(fid=fid, mask=mask, lines=(), **coded)) else: self.outside += 1 self.seconds += time.perf_counter() - t0 @@ -364,6 +384,18 @@ def _finite(geometry: BaseGeometry) -> bool: return bool(np.isfinite(shapely.get_coordinates(geometry)).all()) +def _code(value: Any, what: str) -> int: + """R5: a class code is an integer in 1 .. 2**31 - 1; anything else is refused.""" + try: + code = int(str(value)) if not isinstance(value, list) else 0 + except ValueError: + code = 0 + if not 1 <= code < 2**31: + shown = value[0] if isinstance(value, list) and value else value + raise FeatureError(f"{what}: {shown!r} is not a class code in 1 .. 2**31 - 1") + return code + + def _mask(cmap: ClassMap, vocabulary: EdgeVocabulary, value: Any, what: str) -> int | None: """R4: the mask ``value`` maps to, None for ``drop``; a list is a union.""" names: list[str] = [] diff --git a/src_python/tin_engine/io/ply.py b/src_python/tin_engine/io/ply.py index fa49b630..9fee4f74 100644 --- a/src_python/tin_engine/io/ply.py +++ b/src_python/tin_engine/io/ply.py @@ -61,6 +61,7 @@ def write_ply( ascii: bool = True, comments: Sequence[str] = (), vocabulary: EdgeVocabulary | None = None, + face_codes: npt.ArrayLike | None = None, ) -> bytes: """Encode one mesh as a PLY file. @@ -75,6 +76,9 @@ def write_ply( `feature_bit ` comments sorted by bit and a `feature_vocabulary ` comment (increment 13, ruling 9), after `comments`. + face_codes: `(T,)` land-cover codes, written as the face property + `int land_cover_code` after `vertex_indices` (increment 16c, R3). + The comment naming the code system is the caller's, in `comments`. Raises: ValueError: if `vertices` is not `(N, 3)`; if `faces` and `edges` are @@ -83,7 +87,8 @@ def write_ply( contains a control character, any of which forges a header line in a line-oriented format; if a comment is not ASCII, which the header's encoding cannot carry; or if a mask carries a bit - `vocabulary` does not name. + `vocabulary` does not name; if `face_codes` is given with `edges`, + has not one entry per face, or does not fit int32. """ points = np.ascontiguousarray(vertices, dtype=" bytes: return _lines(" ".join(repr(float(c)) for c in point) for point in points) -def _face_block(faces: npt.NDArray[np.uint32], ascii: bool) -> tuple[list[str], bytes]: - """The 2D mesh's block: one variable-length index list per triangle.""" +def _face_block( + faces: npt.NDArray[np.uint32], codes: npt.ArrayLike | None, ascii: bool +) -> tuple[list[str], bytes]: + """The 2D mesh's block: one variable-length index list per triangle, and + optionally each triangle's land-cover code after it (increment 16c, R3).""" + faces = faces.reshape(-1, 3) declaration = [ f"element face {len(faces)}", f"property list uchar {_INDEX} vertex_indices", ] + columns = [np.full(len(faces), 3), faces] + layout: list[tuple[str, str] | tuple[str, str, int]] = [("n", "u1"), ("v", "= 2**31): + raise ValueError("face_codes must fit int32") + declaration.append("property int land_cover_code") + columns.append(values) + layout.append(("c", " bytes: """Encode one mesh and its constraint edges as legacy VTK 4.2 PolyData. @@ -71,12 +80,18 @@ def write_vtk( vocabulary: what each bit of a mask means. fields: extra dataset strings, such as `("crs", "EPSG:25833")`. binary: write packed big-endian records instead of text. + triangle_codes: `(T,)` land-cover codes (increment 16c, R3), written + as the `int` cell array `land_cover_code`, 0 on every line. + land_cover_codes: what the codes are, written as the dataset string + `land_cover_codes` when non-empty and `triangle_codes` is given. Raises: ValueError: if `vertices` is not `(N, 3)`; if `edge_masks` does not have one entry per edge; if a mask carries a bit the vocabulary does not name; if a field name is outside `^[a-z][a-z0-9_]*$` or - reserved; or if a string is not ASCII or has a control character. + reserved; if a string is not ASCII or has a control character; or + if `triangle_codes` has not one entry per triangle, does not fit + int32, or comes with a vocabulary naming `land_cover_code`. """ points = np.asarray(vertices, dtype=np.float64) if points.ndim != 2 or points.shape[1] != 3: @@ -91,6 +106,17 @@ def write_vtk( raise ValueError(f"field name {name!r} is reserved or not ^[a-z][a-z0-9_]*$") table = sorted((prop.bit, prop.name) for prop in vocabulary.properties) + codes = None + if triangle_codes is not None: + codes = np.asarray(triangle_codes).reshape(-1) + if len(codes) != len(polygons): + raise ValueError(f"triangle_codes has {len(codes)} entries, {len(polygons)} triangles") + if len(codes) and (codes.min() < -(2**31) or codes.max() >= 2**31): + raise ValueError("triangle_codes must fit int32") + if any(name == LAND_COVER_CODE for _, name in table): + raise ValueError(f"the vocabulary names {LAND_COVER_CODE!r}, the codes' array") + if land_cover_codes: + fields = (*fields, ("land_cover_codes", land_cover_codes)) used = {name for mask in np.unique(masks) for name in vocabulary.names(int(mask))} cell_masks = np.concatenate([masks, np.zeros(len(polygons), dtype=np.uint32)]) @@ -105,6 +131,10 @@ def write_vtk( for bit, name in table if name in used ] + if codes is not None: + # Ruling 5: every cell array covers every cell; the lines get 0. + cell_codes = np.concatenate([np.zeros(len(lines), dtype=np.int64), codes]) + features.append(_numeric(LAND_COVER_CODE, cell_codes, "int", binary)) out = [ b"# vtk DataFile Version 4.2\nrasputin mesh\n", diff --git a/src_python/tin_engine/landcover.py b/src_python/tin_engine/landcover.py new file mode 100644 index 00000000..b37d9b17 --- /dev/null +++ b/src_python/tin_engine/landcover.py @@ -0,0 +1,133 @@ +"""A land-cover class per triangle: components first, one point per component. + +Increment 16c (`docs/increments/16c-landcover-labels.md`, R1). Pure: numpy and +shapely in, arrays and counts out; nothing first-party, never `_core`, never a +path. + +The triangles are split into connected components across every edge that is +not a constraint edge (the spread of Shewchuk's Triangle `-A`, blocked by +segments). Each component is then looked up once, at the incentre of its +triangle with the largest inradius: that point is at least `r` from every +constraint edge, so when `r` exceeds `margin` (which bounds how far the noder +moved an input boundary) it lies strictly on one side of every input +boundary. A component whose largest `r` is at most `margin` is still labelled +by its point, and counted `thin`. + +Overlaps (Default D2): the polygon of smallest area wins, ties to the smaller +code, so the result does not depend on the polygons' order. +""" + +from __future__ import annotations + +from collections.abc import Sequence +from dataclasses import dataclass + +import numpy as np +import numpy.typing as npt +import shapely +from shapely.geometry.base import BaseGeometry + + +@dataclass(frozen=True) +class CoverLabels: + """The codes, and the counts of R3's stderr line.""" + + codes: npt.NDArray[np.int32] + regions: int + outside: int + overlapped: int + thin: int + + +def regions(triangles: npt.ArrayLike, edges: npt.ArrayLike) -> npt.NDArray[np.int64]: + """Each triangle's component across unconstrained edges, as the smallest + triangle index in that component.""" + tri = np.asarray(triangles, dtype=np.int64).reshape(-1, 3) + cut = np.asarray(edges, dtype=np.int64).reshape(-1, 2) + n = int(max(tri.max(initial=-1), cut.max(initial=-1))) + 1 + sides = np.concatenate([tri[:, [0, 1]], tri[:, [1, 2]], tri[:, [2, 0]]]) + keys = _keys(sides, n) + owner = np.tile(np.arange(len(tri), dtype=np.int64), 3) + free = ~np.isin(keys, _keys(cut, n)) + order = np.argsort(keys[free], kind="stable") + keys, owner = keys[free][order], owner[free][order] + # A manifold mesh: an interior edge's key appears exactly twice. + pair = np.flatnonzero(keys[1:] == keys[:-1]) + u, v = owner[pair], owner[pair + 1] + parent = np.arange(len(tri), dtype=np.int64) + while True: + ru, rv = parent[u], parent[v] + if np.array_equal(ru, rv): + return parent + # Hook the larger root onto the smaller: labels only fall, no cycle. + np.minimum.at(parent, np.maximum(ru, rv), np.minimum(ru, rv)) + while not np.array_equal(parent, jumped := parent[parent]): + parent = jumped + + +def label_triangles( + vertices: npt.ArrayLike, + triangles: npt.ArrayLike, + edges: npt.ArrayLike, + *, + polygons: Sequence[tuple[BaseGeometry, int]], + margin: float, +) -> CoverLabels: + """A code per triangle from `(polygon, code)` pairs; 0 in no polygon. + + Args: + vertices: `(N, 2)` or `(N, 3)`; only x and y are used. + triangles: `(T, 3)` vertex indices. + edges: `(E, 2)` constraint edges; they block the spread. + polygons: the coded polygons, in the vertices' CRS. + margin: how far a constraint edge may lie from its input boundary. + """ + xy = np.asarray(vertices, dtype=np.float64)[:, :2] + tri = np.asarray(triangles, dtype=np.int64).reshape(-1, 3) + ids = regions(tri, edges) + a, b, c = xy[tri[:, 0]], xy[tri[:, 1]], xy[tri[:, 2]] + la, lb, lc = (np.hypot(*(q - p).T) for p, q in ((b, c), (c, a), (a, b))) + perimeter = la + lb + lc + cross = (b[:, 0] - a[:, 0]) * (c[:, 1] - a[:, 1]) - (c[:, 0] - a[:, 0]) * (b[:, 1] - a[:, 1]) + safe = np.where(perimeter > 0, perimeter, 1.0) + r = np.where(perimeter > 0, np.abs(cross) / safe, 0.0) + centre = (la[:, None] * a + lb[:, None] * b + lc[:, None] * c) / safe[:, None] + # Per component, the triangle of largest r, ties to the lowest index. + order = np.lexsort((np.arange(len(tri)), -r, ids)) + best = order[np.unique(ids[order], return_index=True)[1]] + comp_codes, hits = _lookup(centre[best], polygons) + by_id = np.zeros(len(tri), dtype=np.int32) + by_id[ids[best]] = comp_codes + codes = by_id[ids] + return CoverLabels( + codes=codes, + regions=len(best), + outside=int((hits == 0).sum()), + overlapped=int((hits > 1).sum()), + thin=int((r[best] <= margin).sum()), + ) + + +def _keys(pairs: npt.NDArray[np.int64], n: int) -> npt.NDArray[np.int64]: + """Undirected edge keys, `min * n + max`.""" + return np.minimum(pairs[:, 0], pairs[:, 1]) * n + np.maximum(pairs[:, 0], pairs[:, 1]) + + +def _lookup( + points: npt.NDArray[np.float64], polygons: Sequence[tuple[BaseGeometry, int]] +) -> tuple[npt.NDArray[np.int64], npt.NDArray[np.int64]]: + """Each point's code by D2, and how many polygons it lies in.""" + codes = np.zeros(len(points), dtype=np.int64) + hits = np.zeros(len(points), dtype=np.int64) + if not polygons or not len(points): + return codes, hits + tree = shapely.STRtree([p for p, _ in polygons]) + found, which = tree.query(shapely.points(points), predicate="intersects") + rank = np.array([(p.area, code) for p, code in polygons], dtype=[("a", "f8"), ("c", "i8")]) + # Sort the hits by point, then area, then code: each point's first hit wins. + order = np.lexsort((rank["c"][which], rank["a"][which], found)) + found, which = found[order], which[order] + np.add.at(hits, found, 1) + head = np.unique(found, return_index=True)[1] + codes[found[head]] = rank["c"][which[head]] + return codes, hits diff --git a/src_python/tin_engine/palettes.py b/src_python/tin_engine/palettes.py new file mode 100644 index 00000000..c56250f9 --- /dev/null +++ b/src_python/tin_engine/palettes.py @@ -0,0 +1,86 @@ +"""Natural colours for land-cover codes, and a ParaView preset from them. + +Increment 16c (`docs/increments/16c-landcover-labels.md`, R4). Data plus one +function; nothing first-party, nothing about meshes. The colours follow the +natural-colour convention of Patterson and Kelso (2004): the colour the ground +has from above, not the official CLC legend (which paints peat bogs blue and +the sea almost white). They are taste, and Ola's to change (Default D6). +""" + +from __future__ import annotations + +#: CORINE level-3 code -> (CLC `LABEL3`, `#rrggbb`); 0 is no polygon, and +#: every constraint line in a `.vtk`. +# fmt: off +CORINE_NATURAL: dict[int, tuple[str, str]] = { + 0: ("No polygon; constraint lines", "#3c3c3c"), + 111: ("Continuous urban fabric", "#8c5f5a"), + 112: ("Discontinuous urban fabric", "#a88d86"), + 121: ("Industrial or commercial units", "#8e8a96"), + 122: ("Road and rail networks and associated land", "#6e6e6e"), + 123: ("Port areas", "#7d8796"), + 124: ("Airports", "#a8a8a8"), + 131: ("Mineral extraction sites", "#b49b78"), + 132: ("Dump sites", "#8a7d64"), + 133: ("Construction sites", "#bcae98"), + 141: ("Green urban areas", "#86b86e"), + 142: ("Sport and leisure facilities", "#a3cf7e"), + 211: ("Non-irrigated arable land", "#e6d58c"), + 212: ("Permanently irrigated land", "#d9cc6e"), + 213: ("Rice fields", "#cdd89a"), + 221: ("Vineyards", "#9c7a44"), + 222: ("Fruit trees and berry plantations", "#a9ad5e"), + 223: ("Olive groves", "#8f9a52"), + 231: ("Pastures", "#b4d47a"), + 241: ("Annual crops associated with permanent crops", "#dccf94"), + 242: ("Complex cultivation patterns", "#d2c47c"), + 243: ("Land principally occupied by agriculture, with significant areas of natural vegetation", "#bcc47e"), # noqa: E501 + 244: ("Agro-forestry areas", "#a9b574"), + 311: ("Broad-leaved forest", "#4f8f3f"), + 312: ("Coniferous forest", "#1f5a2e"), + 313: ("Mixed forest", "#357438"), + 321: ("Natural grasslands", "#c2d68a"), + 322: ("Moors and heathland", "#a7c47f"), + 323: ("Sclerophyllous vegetation", "#8a9658"), + 324: ("Transitional woodland-shrub", "#7ea65a"), + 331: ("Beaches, dunes, sands", "#e9ddb2"), + 332: ("Bare rocks", "#8f8f8f"), + 333: ("Sparsely vegetated areas", "#b9b8a0"), + 334: ("Burnt areas", "#4b3f3a"), + 335: ("Glaciers and perpetual snow", "#eef6fb"), + 411: ("Inland marshes", "#6e9470"), + 412: ("Peat bogs", "#8a6642"), + 421: ("Salt marshes", "#7f9f8c"), + 422: ("Salines", "#d8d6cc"), + 423: ("Intertidal flats", "#b3bdb3"), + 511: ("Water courses", "#4c8ec4"), + 512: ("Water bodies", "#3e7bb6"), + 521: ("Coastal lagoons", "#5b93b3"), + 522: ("Estuaries", "#5188b4"), + 523: ("Sea and ocean", "#2b5d8e"), +} +# fmt: on + +#: The palettes `rasputin palette` knows, by name, with their preset names. +PALETTES = {"corine": (CORINE_NATURAL, "rasputin CORINE natural")} + +#: A code the table lacks draws magenta, so it is seen rather than blended in. +_NAN_COLOR = [1.0, 0.0, 1.0] + + +def paraview_preset(table: dict[int, tuple[str, str]], name: str) -> list[dict[str, object]]: + """A categorical ParaView colour preset: `Annotations` (value, label, ...) + and `IndexedColors` (r, g, b in 0 .. 1, one triple per value, in order).""" + annotations: list[str] = [] + colours: list[float] = [] + for code, (label, colour) in sorted(table.items()): + annotations += [str(code), f"{code} {label}"] + colours += [int(colour[i : i + 2], 16) / 255 for i in (1, 3, 5)] + return [ + { + "Name": name, + "Annotations": annotations, + "IndexedColors": colours, + "NanColor": _NAN_COLOR, + } + ] diff --git a/tests/python/landcover_fixtures.py b/tests/python/landcover_fixtures.py new file mode 100644 index 00000000..460fd6d3 --- /dev/null +++ b/tests/python/landcover_fixtures.py @@ -0,0 +1,172 @@ +"""Test support for increment 16c: the land-cover oracle and the spread check. + +`docs/increments/16c-landcover-labels.md`, R1 ("The oracle (tests only)") and +the invariants I1 and I2. Everything here is built from the triangles and the +input polygons, never from the producer's records: no component ids, no chosen +points. It borrows only the producer's *predicate* (a point against the input +polygons, `intersects`) and its overlap rule (smallest area, ties to the +smaller code; Default D2). + +- `landcover_oracle`: each triangle's centroid tested against the polygons, + for every triangle whose inradius exceeds `1.5 * margin` (the centroid is + then further than `margin` from every input boundary, so the answer is + exact). The legacy's per-centre test (`legacy/rasputin/gml_repository.py`, + `land_cover`), vectorised, with 0 for a centre in no polygon. +- `spread_violations`: I1. Every interior edge that is not a constraint edge + has the same code on both sides. Arrays only, no geometry. +- `vtk_labels`: I3, I1 and I2 on one `.vtk` file, for the CLI suites. +""" + +from __future__ import annotations + +from collections.abc import Sequence +from dataclasses import dataclass + +import numpy as np +import numpy.typing as npt +import shapely +from shapely.geometry.base import BaseGeometry + +from vtkread import VtkFile, lines_as_array, polygons_as_array + +#: The noder's snap spacing at the CLI's default (`cli.DEFAULT_SNAP_SPACING`). +SNAP = 1e-3 +#: R1: `margin = 2h`. +MARGIN = 2 * SNAP + +Coded = Sequence[tuple[BaseGeometry, int]] + + +def xy_of(vertices: npt.ArrayLike) -> np.ndarray: + return np.asarray(vertices, dtype=np.float64)[:, :2] + + +def inradius(vertices: npt.ArrayLike, triangles: npt.ArrayLike) -> np.ndarray: + """`2 * area / perimeter` per triangle, in x and y.""" + xy = xy_of(vertices) + tri = np.asarray(triangles, dtype=np.int64).reshape(-1, 3) + a, b, c = xy[tri[:, 0]], xy[tri[:, 1]], xy[tri[:, 2]] + area = 0.5 * np.abs( + (b[:, 0] - a[:, 0]) * (c[:, 1] - a[:, 1]) - (c[:, 0] - a[:, 0]) * (b[:, 1] - a[:, 1]) + ) + perimeter = np.hypot(*(b - c).T) + np.hypot(*(c - a).T) + np.hypot(*(a - b).T) + return np.divide(2 * area, perimeter, out=np.zeros_like(area), where=perimeter > 0) + + +def centroids(vertices: npt.ArrayLike, triangles: npt.ArrayLike) -> np.ndarray: + xy = xy_of(vertices) + tri = np.asarray(triangles, dtype=np.int64).reshape(-1, 3) + return xy[tri].mean(axis=1) + + +def pick(hits: Sequence[int], polygons: Coded) -> int: + """D2: the smallest area wins, ties to the smaller code; none is 0.""" + if not hits: + return 0 + return min((polygons[i][0].area, polygons[i][1]) for i in hits)[1] + + +def codes_at(points: np.ndarray, polygons: Coded) -> tuple[np.ndarray, np.ndarray]: + """Each point's code by D2, and how many polygons it is in.""" + codes = np.zeros(len(points), dtype=np.int64) + counts = np.zeros(len(points), dtype=np.int64) + if not len(polygons) or not len(points): + return codes, counts + tree = shapely.STRtree([p for p, _ in polygons]) + found, hit = tree.query(shapely.points(points), predicate="intersects") + per_point: dict[int, list[int]] = {} + for i, j in zip(found.tolist(), hit.tolist(), strict=True): + per_point.setdefault(i, []).append(j) + for i, hits in per_point.items(): + codes[i] = pick(hits, polygons) + counts[i] = len(hits) + return codes, counts + + +@dataclass(frozen=True) +class Oracle: + """The expected code of each triangle, and which triangles it is exact for.""" + + codes: np.ndarray + checked: np.ndarray + + def mismatches(self, produced: npt.ArrayLike) -> list[tuple[int, int, int]]: + got = np.asarray(produced).reshape(-1) + bad = np.flatnonzero(self.checked & (got != self.codes)) + return [(int(t), int(got[t]), int(self.codes[t])) for t in bad[:20]] + + +def landcover_oracle( + vertices: npt.ArrayLike, triangles: npt.ArrayLike, polygons: Coded, margin: float = MARGIN +) -> Oracle: + """R1's oracle: centroids against the polygons, where `r > 1.5 * margin`.""" + codes, _ = codes_at(centroids(vertices, triangles), polygons) + checked = inradius(vertices, triangles) > 1.5 * margin + return Oracle(codes=codes, checked=checked) + + +def edge_keys(pairs: npt.ArrayLike, n: int) -> np.ndarray: + e = np.asarray(pairs, dtype=np.int64).reshape(-1, 2) + return np.minimum(e[:, 0], e[:, 1]) * n + np.maximum(e[:, 0], e[:, 1]) + + +def spread_violations( + triangles: npt.ArrayLike, constraints: npt.ArrayLike, codes: npt.ArrayLike +) -> list[tuple[int, int, int, int]]: + """I1: `(t, u, code_t, code_u)` for each unconstrained interior edge whose + two triangles differ in code.""" + tri = np.asarray(triangles, dtype=np.int64).reshape(-1, 3) + got = np.asarray(codes).reshape(-1) + n = int(tri.max()) + 1 if len(tri) else 0 + sides = np.concatenate([tri[:, [0, 1]], tri[:, [1, 2]], tri[:, [2, 0]]]) + keys = edge_keys(sides, n) + owner = np.tile(np.arange(len(tri)), 3) + free = ~np.isin(keys, edge_keys(constraints, n)) + keys, owner = keys[free], owner[free] + order = np.argsort(keys, kind="stable") + keys, owner = keys[order], owner[order] + pair = np.flatnonzero(keys[1:] == keys[:-1]) + t, u = owner[pair], owner[pair + 1] + bad = np.flatnonzero(got[t] != got[u]) + return [(int(t[k]), int(u[k]), int(got[t[k]]), int(got[u[k]])) for k in bad[:20]] + + +def interior_edge_count(triangles: npt.ArrayLike) -> int: + """How many edges have two triangles: the probe I1 needs to mean anything.""" + tri = np.asarray(triangles, dtype=np.int64).reshape(-1, 3) + n = int(tri.max()) + 1 + keys = edge_keys(np.concatenate([tri[:, [0, 1]], tri[:, [1, 2]], tri[:, [2, 0]]]), n) + _, counts = np.unique(keys, return_counts=True) + return int((counts == 2).sum()) + + +def triangle_areas(vertices: npt.ArrayLike, triangles: npt.ArrayLike) -> np.ndarray: + xy = xy_of(vertices) + tri = np.asarray(triangles, dtype=np.int64).reshape(-1, 3) + a, b, c = xy[tri[:, 0]], xy[tri[:, 1]], xy[tri[:, 2]] + return 0.5 * np.abs( + (b[:, 0] - a[:, 0]) * (c[:, 1] - a[:, 1]) - (c[:, 0] - a[:, 0]) * (b[:, 1] - a[:, 1]) + ) + + +def vtk_labels( + vtk: VtkFile, polygons: Coded, margin: float = MARGIN +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """I3, I1 and I2 on one `.vtk`; returns its lines, triangles and the + triangles' codes.""" + lines, triangles = lines_as_array(vtk), polygons_as_array(vtk) + array = vtk.cell_array("land_cover_code") + values = np.asarray(array.values) + # I3: an int scalar, one value per cell, 0 on every line. + assert (array.type_name, array.components) == ("int", 1) + assert len(values) == vtk.cells == len(lines) + len(triangles) + assert not values[: len(lines)].any() + codes = values[len(lines) :] + # I1: the spread, on the arrays alone. + assert interior_edge_count(triangles) > 0 + assert spread_violations(triangles, lines, codes) == [] + # I2: the oracle, from the input polygons; it must reach most triangles. + oracle = landcover_oracle(vtk.points, triangles, polygons, margin) + assert oracle.checked.sum() > 0.9 * len(triangles) + assert oracle.mismatches(codes) == [] + return lines, triangles, codes diff --git a/tests/python/test_cli_mesh_features.py b/tests/python/test_cli_mesh_features.py index 57fda6ea..a67f8016 100644 --- a/tests/python/test_cli_mesh_features.py +++ b/tests/python/test_cli_mesh_features.py @@ -73,6 +73,7 @@ from gpkg_fixtures import ( DTM10, EXTRACT, + EXTRACT_TABLE, LEGACY_GML, OLA_EUROPE, OLA_NORWAY, @@ -82,6 +83,7 @@ needs_rtree, write_gpkg, ) +from landcover_fixtures import triangle_areas, vtk_labels from plyread import read_ply from test_cli_mesh_dem import USAGE, invoke, write_tiff from test_cli_mesh_domain import COLS, ROWS, SQUARE, geojson, quarter_circle @@ -589,7 +591,12 @@ def _max_error(tile: Any, points: np.ndarray, triangles: np.ndarray, domain: Pol class TestCommittedExtract: """Acceptance 16b-1/2 in CI: the quarter circle on the committed tile with the committed extract, `--features-map corine`, at 10 m (M5's reference: - about 54 840 triangles). I7: the tolerance holds, every node covered.""" + about 54 840 triangles). I7: the tolerance holds, every node covered. + + Increment 16c adds I1, I2, I3 and I5 on the same run (the `test_16c_` + tests, committed red at `196147e`: the file had no `land_cover_code`; + green since `0487ed0`); its I4 is + `test_i3_rows_in_reverse_order_give_the_same_file`, unchanged.""" @pytest.fixture(scope="class") def run(self, tmp_path_factory: pytest.TempPathFactory) -> tuple[VtkFile, str, Path]: @@ -659,6 +666,35 @@ def test_i3_rows_in_reverse_order_give_the_same_file(self, run: Any, tmp_path: P assert code == 0, output assert target.read_bytes() == first.read_bytes() + @needs_codecs + def test_16c_i1_i2_i3_land_cover_codes(self, run: Any) -> None: + """Increment 16c on the real mesh: every cell carries a code, 0 on the + lines (I3); the spread (I1); the centroid oracle against the extract's + polygons moved to EPSG:25833 (I2). The extract covers the quarter + circle, so no triangle is 0.""" + vtk, _, _ = run + _, _, codes = vtk_labels(vtk, _extract_polygons()) + assert 0 not in set(codes.tolist()) + + @needs_codecs + def test_16c_i5_class_areas_match_the_clipped_extract(self, run: Any) -> None: + """I5: per code, the triangles' xy area equals the area of each + extract polygon (moved to EPSG:25833) intersected with the quarter + circle, within 0.01 % of the domain's area.""" + vtk, _, _ = run + codes = np.asarray(vtk.cell_array("land_cover_code").values)[len(lines_as_array(vtk)) :] + triangles = polygons_as_array(vtk) + areas = triangle_areas(vtk.points, triangles) + domain = Polygon(quarter_circle()) + expected: dict[int, float] = {} + for polygon, code in _extract_polygons(): + expected[code] = expected.get(code, 0.0) + polygon.intersection(domain).area + got = {int(c): float(areas[codes == c].sum()) for c in np.unique(codes)} + assert len(expected) >= 5 + for code in set(expected) | set(got): + gap = abs(got.get(code, 0.0) - expected.get(code, 0.0)) + assert gap <= 1e-4 * domain.area, (code, got.get(code), expected.get(code)) + @needs_codecs def test_the_legacy_gml(self, tmp_path: Path) -> None: """Q6 (b): the legacy GML (EPSG:4326) with `clc18_kode`, same domain.""" @@ -680,6 +716,14 @@ def test_the_legacy_gml(self, tmp_path: Path) -> None: assert V.mask("land_cover", "water") in set(edges(vtk)[1].tolist()) +def _extract_polygons() -> list[tuple[Any, int]]: + """The committed extract's polygons in the DEM's CRS, with their codes (16c).""" + return [ + (ff.moved(geometry, "EPSG:3035"), int(code)) + for _, geometry, code in ff.extract_features(EXTRACT, EXTRACT_TABLE) + ] + + # ------------------------------------------------------------ Ola's local data diff --git a/tests/python/test_cli_mesh_landcover.py b/tests/python/test_cli_mesh_landcover.py new file mode 100644 index 00000000..65dee2c8 --- /dev/null +++ b/tests/python/test_cli_mesh_landcover.py @@ -0,0 +1,542 @@ +"""`rasputin mesh ... --features-map corine`: a land-cover code per cell (16c). + +`docs/increments/16c-landcover-labels.md`: R3 (what the files carry), R5 (what +`feature_input` keeps), R6 (the CLI) and "Tests for @tester", through the +command on the synthetic `bumpy` DEM of `test_cli_mesh_features.py`, with +GeoJSON features carrying `Code_18`. For every fixture: I1 (the spread, from +the file's `LINES` and triangles), I2 (the centroid oracle of +`landcover_fixtures`, from the input polygons), I3 (one value per cell, 0 on +every line), and the codes the design names for it. + +Pinned beyond the design's text (the design leaves them open): + +- The stderr line is matched as `land cover: regions, outside every + polygon, in more than one, thinner than the snap` (R3's words). +- The `land_cover_codes` string contains `CORINE Land Cover level-3 code`, + `Code_18`, `corine` and `0 =`; the PLY comment is `land_cover_codes + `. +- A refused code is a usage error (exit 2) naming the feature and the value, + and `--features`, as every other `FeatureError` is (16b). +- A `LineString` under a coded map keeps its code and has no polygon. +- `rasputin mesh --help` names `land_cover_code` and `rasputin palette`. + +Committed red at `196147e`: no code wrote `land_cover_code`, +`land_cover_codes` or the stderr line; `rasputin palette` was no command; +`ClassMap` had no `codes` and `TerrainFeature` no `code` or `polygon`; a +covering polygon was dropped as outside; and `corine` accepted any value. +Every test failed on one of those. They landed in `0487ed0` and the suite has +been green since. +""" + +from __future__ import annotations + +import json +import re +from dataclasses import dataclass +from pathlib import Path +from typing import Any + +import numpy as np +import pytest +import shapely +from numpy.testing import assert_array_equal +from shapely.geometry import LineString, Polygon, box +from shapely.geometry.base import BaseGeometry +from typer.testing import CliRunner + +import feature_fixtures as ff +from feature_fixtures import Feat, domain_of, write_geojson +from geotiff_fixtures import TIE_X, TIE_Y, micro_tiff +from landcover_fixtures import MARGIN, landcover_oracle, spread_violations, vtk_labels +from plyread import read_ply +from test_cli_mesh import plain +from test_cli_mesh_dem import USAGE, invoke, write_tiff +from test_cli_mesh_domain import COLS, ROWS, SQUARE, geojson +from tin_engine.cli import app +from tin_engine.palettes import CORINE_NATURAL, paraview_preset +from vtkread import VtkFile, read_vtk + +runner = CliRunner(env={"NO_COLOR": "1", "TERM": "dumb"}) + +LAND_COVER_LINE = re.compile( + r"land cover: (?P\d+) regions, (?P\d+) outside every polygon, " + r"(?P\d+) in more than one, (?P\d+) thinner than the snap" +) +FEATURES_LINE = re.compile( + r"\b(?P\d+) features kept, (?P\d+) dropped outside, (?P\d+) clipped, " + r"(?P\d+) empty skipped\b" +) +CODES_SYSTEM = "CORINE Land Cover level-3 code" +PRESET_NAME = "rasputin CORINE natural" + + +def rel(x: float, y: float) -> tuple[float, float]: + return (TIE_X + x, TIE_Y + y) + + +def rect(x0: float, y0: float, x1: float, y1: float) -> Polygon: + (ax, ay), (bx, by) = rel(x0, y0), rel(x1, y1) + return box(ax, ay, bx, by) + + +def coded(fid: Any, geometry: BaseGeometry, code: Any) -> Feat: + return Feat(fid, geometry, {"Code_18": code}) + + +# Inside `SQUARE` (x 12.3 .. 187.7, y -73.3 .. -6.7 from the tie point), off +# every DEM node. +WEST = rect(40.3, -60.2, 100.7, -20.4) +EAST = rect(100.7, -60.2, 160.9, -20.4) +EAST_OFFSET = rect(100.7001, -60.2, 160.9, -20.4) # 0.1 mm from WEST's east side +LAKE = rect(80.1, -50.2, 120.4, -30.3) +FOREST = rect(30.3, -65.1, 170.2, -15.3) +HOLED_FOREST = Polygon(FOREST.exterior, [LAKE.exterior]) +ROAD = LineString([rel(30.1, -40.3), rel(170.2, -40.3)]) +COVER = rect(-50.3, -150.1, 300.2, 50.4) +#: The domain with a notch cut down from its north side, x 80 .. 120. +NOTCHED: list[tuple[float, float]] = [ + rel(12.3, -73.3), + rel(187.7, -72.9), + rel(186.1, -6.7), + rel(120.3, -6.9), + rel(119.9, -40.1), + rel(80.1, -40.3), + rel(79.7, -7.0), + rel(13.9, -7.1), +] +BAND = rect(50.2, -30.1, 150.3, -15.2) + + +@pytest.fixture +def bumpy(tmp_path: Path) -> Path: + array = np.random.default_rng(16).uniform(0.0, 50.0, (ROWS, COLS)).astype(np.float32) + return write_tiff(tmp_path / "bumpy.tif", micro_tiff(array)) + + +@pytest.fixture +def plain_square(tmp_path: Path) -> Path: + return geojson(tmp_path / "square.geojson", SQUARE) + + +def run( + tmp_path: Path, tif: Path, domain: Path, *extra: str, out: str = "x.vtk" +) -> tuple[int, str, Path]: + target = tmp_path / out + code, output = invoke( + "--dem", str(tif), "--domain", str(domain), "--tolerance", "1", "--out", str(target), *extra + ) + return code, output, target + + +@dataclass(frozen=True) +class Labelled: + """A `.vtk` written with land-cover codes, and what I1-I3 need of it.""" + + vtk: VtkFile + output: str + points: np.ndarray + lines: np.ndarray + triangles: np.ndarray + codes: np.ndarray + + def inside(self, geometry: BaseGeometry) -> np.ndarray: + centroids = self.points[self.triangles][:, :, :2].mean(axis=1) + return np.asarray(shapely.contains_xy(geometry, *centroids.T)) + + def stderr_counts(self) -> dict[str, int]: + found = [m.groupdict() for m in LAND_COVER_LINE.finditer(self.output)] + assert len(found) == 1, self.output + return {k: int(v) for k, v in found[0].items()} + + +def labelled_vtk(vtk: VtkFile, output: str, polygons: list[tuple[BaseGeometry, int]]) -> Labelled: + """I1, I2 and I3 on one file (`landcover_fixtures.vtk_labels`).""" + lines, triangles, codes = vtk_labels(vtk, polygons, MARGIN) + return Labelled(vtk, output, vtk.points, lines, triangles, codes) + + +def polygons_of(features: list[Feat]) -> list[tuple[BaseGeometry, int]]: + return [ + (f.geometry, int(f.properties["Code_18"])) + for f in features + if isinstance(f.geometry, Polygon) + ] + + +def mesh_corine( + tmp_path: Path, tif: Path, features: list[Feat], domain: Path | None = None +) -> Labelled: + path = write_geojson(tmp_path / "corine.geojson", features) + domain = domain or geojson(tmp_path / "square.geojson", SQUARE) + code, output, target = run( + tmp_path, tif, domain, "--features", str(path), "--features-map", "corine" + ) + assert code == 0, output + return labelled_vtk(read_vtk(target.read_bytes()), output, polygons_of(features)) + + +def text_field(vtk: VtkFile, name: str) -> str: + (value,) = vtk.field_data[name].values + return str(value) + + +# ------------------------------------------------------------ the fixtures + + +class TestFixtures: + """The design's seven labelling fixtures, end to end.""" + + def test_two_squares_side_by_side(self, tmp_path: Path, bumpy: Path) -> None: + got = mesh_corine(tmp_path, bumpy, [coded(1, WEST, "311"), coded(2, EAST, "512")]) + assert set(got.codes[got.inside(WEST)].tolist()) == {311} + assert set(got.codes[got.inside(EAST)].tolist()) == {512} + assert set(got.codes[~got.inside(WEST.union(EAST))].tolist()) == {0} + # The triangles along the shared side carry both codes. + seam = LineString([rel(100.7, -55.1), rel(100.7, -25.3)]) + cells = shapely.polygons(got.points[got.triangles][:, :, :2]) + touching = np.asarray(shapely.distance(seam, cells)) <= 1e-9 + assert set(got.codes[touching].tolist()) == {311, 512} + counts = got.stderr_counts() + assert counts["o"] >= 1 and counts["v"] == 0 + + def test_a_forest_with_an_empty_hole(self, tmp_path: Path, bumpy: Path) -> None: + got = mesh_corine(tmp_path, bumpy, [coded("f", HOLED_FOREST, "312")]) + assert got.inside(LAKE).any() + assert set(got.codes[got.inside(LAKE)].tolist()) == {0} + assert set(got.codes[got.inside(HOLED_FOREST)].tolist()) == {312} + + def test_a_lake_in_the_hole_of_a_holed_forest(self, tmp_path: Path, bumpy: Path) -> None: + got = mesh_corine( + tmp_path, bumpy, [coded("f", HOLED_FOREST, "312"), coded("l", LAKE, "512")] + ) + assert set(got.codes[got.inside(LAKE)].tolist()) == {512} + assert set(got.codes[got.inside(HOLED_FOREST)].tolist()) == {312} + assert got.stderr_counts()["v"] == 0 + + def test_a_lake_in_an_unholed_forest(self, tmp_path: Path, bumpy: Path) -> None: + """D2: the smaller polygon wins, and stderr counts the overlap.""" + got = mesh_corine(tmp_path, bumpy, [coded("f", FOREST, "312"), coded("l", LAKE, "512")]) + assert set(got.codes[got.inside(LAKE)].tolist()) == {512} + assert set(got.codes[got.inside(FOREST.difference(LAKE))].tolist()) == {312} + assert got.stderr_counts()["v"] > 0 + + def test_a_polygon_clipped_by_the_domain_into_two_pieces( + self, tmp_path: Path, bumpy: Path + ) -> None: + domain = geojson(tmp_path / "notched.geojson", NOTCHED) + got = mesh_corine(tmp_path, bumpy, [coded("b", BAND, "324")], domain=domain) + west = got.inside(BAND.intersection(rect(12, -80, 80.1, 0))) + east = got.inside(BAND.intersection(rect(119.9, -80, 190, 0))) + assert west.any() and east.any() + assert set(got.codes[west | east].tolist()) == {324} + + def test_a_road_across_a_polygon(self, tmp_path: Path, bumpy: Path) -> None: + got = mesh_corine(tmp_path, bumpy, [coded("f", WEST, "311"), coded("r", ROAD, "122")]) + north = got.inside(WEST.intersection(rect(0, -40.3, 200, 0))) + south = got.inside(WEST.intersection(rect(0, -80, 200, -40.3))) + assert north.any() and south.any() + assert set(got.codes[north | south].tolist()) == {311} + + def test_a_polygon_covering_the_domain(self, tmp_path: Path, bumpy: Path) -> None: + """R5: kept with no lines, counted kept and not clipped.""" + got = mesh_corine(tmp_path, bumpy, [coded("c", COVER, "333")]) + assert set(got.codes.tolist()) == {333} + found = [m.groupdict() for m in FEATURES_LINE.finditer(got.output)] + assert found == [{"n": "1", "o": "0", "c": "0", "e": "0"}], got.output + assert got.stderr_counts()["o"] == 0 + + def test_boundaries_a_tenth_of_a_millimetre_apart(self, tmp_path: Path, bumpy: Path) -> None: + """Less than the snap apart: the run succeeds and, away from the seam, + each square carries its own code (I2).""" + got = mesh_corine(tmp_path, bumpy, [coded(1, WEST, "311"), coded(2, EAST_OFFSET, "512")]) + assert set(got.codes[got.inside(rect(40.3, -60.2, 100.6, -20.4))].tolist()) == {311} + assert set(got.codes[got.inside(rect(100.8, -60.2, 160.9, -20.4))].tolist()) == {512} + + def test_a_code_given_as_a_json_number(self, tmp_path: Path, bumpy: Path) -> None: + got = mesh_corine(tmp_path, bumpy, [coded(1, WEST, 311)]) + assert set(got.codes[got.inside(WEST)].tolist()) == {311} + + +# ------------------------------------------------------------ the record + + +class TestRecord: + @pytest.fixture + def two_squares(self, tmp_path: Path) -> Path: + return write_geojson( + tmp_path / "two.geojson", [coded(1, WEST, "311"), coded(2, EAST, "512")] + ) + + def test_the_codes_field( + self, tmp_path: Path, bumpy: Path, plain_square: Path, two_squares: Path + ) -> None: + code, output, target = run( + tmp_path, + bumpy, + plain_square, + "--features", + str(two_squares), + "--features-map", + "corine", + ) + assert code == 0, output + text = text_field(read_vtk(target.read_bytes()), "land_cover_codes") + for words in (CODES_SYSTEM, "Code_18", "corine", "0 ="): + assert words in text, text + + @pytest.mark.parametrize("binary", [False, True], ids=["ascii", "binary"]) + def test_every_cell_carries_a_code_in_both_encodings( + self, tmp_path: Path, bumpy: Path, plain_square: Path, two_squares: Path, binary: bool + ) -> None: + extra = ["--binary"] if binary else [] + code, output, target = run( + tmp_path, + bumpy, + plain_square, + "--features", + str(two_squares), + "--features-map", + "corine", + *extra, + ) + assert code == 0, output + got = labelled_vtk(read_vtk(target.read_bytes()), output, [(WEST, 311), (EAST, 512)]) + assert {311, 512} <= set(got.codes.tolist()) + assert list(got.vtk.scalars) == ["feature_mask"] + + def test_no_coded_map_no_codes( + self, tmp_path: Path, bumpy: Path, plain_square: Path, two_squares: Path + ) -> None: + """D3: `property` writes no labels, nor does a run without features; + `corine` over the same polygons does (the control).""" + by_property = write_geojson( + tmp_path / "p.geojson", + [Feat(1, WEST, {"property": "land_cover"}), Feat(2, EAST, {"property": "water"})], + ) + runs = { + "none": (), + "property": ("--features", str(by_property)), + "corine": ("--features", str(two_squares), "--features-map", "corine"), + } + files = {} + for name, extra in runs.items(): + (tmp_path / name).mkdir() + code, output, target = run(tmp_path / name, bumpy, plain_square, *extra) + assert code == 0, output + files[name] = (read_vtk(target.read_bytes()), output) + for name in ("none", "property"): + vtk, output = files[name] + with pytest.raises(KeyError): + vtk.cell_array("land_cover_code") + assert "land_cover_codes" not in vtk.field_data + assert not LAND_COVER_LINE.search(output), output + vtk, output = files["corine"] + assert "land_cover_codes" in vtk.field_data + assert LAND_COVER_LINE.search(output), output + + def test_the_stats_phase( + self, tmp_path: Path, bumpy: Path, plain_square: Path, two_squares: Path + ) -> None: + code, output, _ = run( + tmp_path, + bumpy, + plain_square, + "--features", + str(two_squares), + "--features-map", + "corine", + "--stats", + "-", + ) + assert code == 0, output + assert re.search(r"\|\s*land cover\s*\|", output), output + + def test_the_help_names_the_array_and_the_palette(self) -> None: + result = runner.invoke(app, ["mesh", "--help"]) + text = plain(result.output) + assert "land_cover_code" in text and "rasputin palette" in text, text + + +# ------------------------------------------------------------ refusals + + +class TestRefusals: + @pytest.mark.parametrize( + "value", + ["forest", "0", "-311", "2147483648", "3.5", ["311", "312"]], + ids=["word", "zero", "negative", "over-int32", "decimal", "list"], + ) + def test_a_code_that_is_not_a_class_code_is_refused( + self, tmp_path: Path, bumpy: Path, plain_square: Path, value: Any + ) -> None: + """R5: `int(str(value))` in 1 .. 2**31 - 1, else refused naming the + feature and the value; a list is refused too.""" + path = write_geojson(tmp_path / "bad.geojson", [coded("f-77", WEST, value)]) + code, output, target = run( + tmp_path, bumpy, plain_square, "--features", str(path), "--features-map", "corine" + ) + assert code == USAGE, output + assert "f-77" in output and "--features" in output, output + shown = value[0] if isinstance(value, list) else value + assert str(shown) in output, output + assert not target.exists() + + +# ------------------------------------------------------------ .ply + + +class TestPly: + def test_the_face_file_carries_the_codes_and_the_edge_file_not( + self, tmp_path: Path, bumpy: Path, plain_square: Path + ) -> None: + path = write_geojson( + tmp_path / "two.geojson", [coded(1, WEST, "311"), coded(2, EAST, "512")] + ) + edges_path = tmp_path / "e.ply" + code, output, target = run( + tmp_path, + bumpy, + plain_square, + "--features", + str(path), + "--features-map", + "corine", + "--out-edges", + str(edges_path), + out="x.ply", + ) + assert code == 0, output + header, data = read_ply(target.read_bytes()) + face = header.element("face") + assert [p.name for p in face.properties] == ["vertex_indices", "land_cover_code"] + assert face.properties[1].type_name == "int" + texts = [c for c in header.comments if c.startswith("land_cover_codes ")] + assert len(texts) == 1 and CODES_SYSTEM in texts[0], header.comments + + edge_header, edge_data = read_ply(edges_path.read_bytes()) + assert all(p.name != "land_cover_code" for e in edge_header.elements for p in e.properties) + assert not any(c.startswith("land_cover_codes") for c in edge_header.comments) + + # I1 and I2 on the face file, with the edge file's constraints. + vertex = data["vertex"] + points = np.column_stack([vertex["x"], vertex["y"], vertex["z"]]) + triangles = np.stack(list(data["face"]["vertex_indices"])).astype(np.int64) + codes = np.asarray(data["face"]["land_cover_code"]) + lines = np.column_stack([edge_data["edge"]["vertex1"], edge_data["edge"]["vertex2"]]) + assert spread_violations(triangles, lines, codes) == [] + oracle = landcover_oracle(points, triangles, [(WEST, 311), (EAST, 512)], MARGIN) + assert oracle.mismatches(codes) == [] + assert {311, 512} <= set(codes.tolist()) + + +# ------------------------------------------------------------ VTK readback + + +def test_vtk_reads_the_array_back(tmp_path: Path, bumpy: Path, plain_square: Path) -> None: + """`vtkPolyDataReader` sees `land_cover_code` with one value per cell, + 0 on the lines (which VTK numbers first). Skipped without `vtk`, as + `test_io_vtk_readback.py` is; CI runs it in the `viewer` extra's step.""" + vtk = pytest.importorskip("vtk") + path = write_geojson(tmp_path / "two.geojson", [coded(1, WEST, "311"), coded(2, EAST, "512")]) + code, output, target = run( + tmp_path, bumpy, plain_square, "--features", str(path), "--features-map", "corine" + ) + assert code == 0, output + reader = vtk.vtkPolyDataReader() + errors: list[str] = [] + reader.AddObserver("ErrorEvent", lambda _obj, _event: errors.append("error")) + reader.SetFileName(str(target)) + reader.ReadAllFieldsOn() + reader.Update() + assert not errors + polydata = reader.GetOutput() + array = polydata.GetCellData().GetArray("land_cover_code") + assert array is not None + values = np.array([array.GetValue(i) for i in range(array.GetNumberOfTuples())]) + assert len(values) == polydata.GetNumberOfCells() + lines = polydata.GetNumberOfLines() + assert lines > 0 and not values[:lines].any() + assert {311, 512} <= set(values[lines:].tolist()) + ours = read_vtk(target.read_bytes()).cell_array("land_cover_code").values + assert_array_equal(values, np.asarray(ours)) + + +# ------------------------------------------------------------ the palette + + +def test_palette_out_writes_the_preset(tmp_path: Path) -> None: + """R4: `rasputin palette corine --out FILE` writes the ParaView preset.""" + target = tmp_path / "corine_natural.json" + result = runner.invoke(app, ["palette", "corine", "--out", str(target)]) + assert result.exit_code == 0, result.output + assert json.loads(target.read_text()) == paraview_preset(CORINE_NATURAL, PRESET_NAME) + + +# ------------------------------------------------- what feature_input keeps (R5) + + +class TestWhatFeatureInputKeeps: + DOMAIN = domain_of(Polygon(SQUARE)) + + def open(self, tmp_path: Path, features: list[Feat], map_name: str, **kw: Any) -> Any: + path = write_geojson(tmp_path / f"{map_name}.geojson", features, **kw) + return ff.open_one(path, self.DOMAIN, map_name) + + @pytest.mark.parametrize( + ("name", "codes"), + [ + ("corine", CODES_SYSTEM), + ("corine-water", CODES_SYSTEM), + ("clc18_kode", CODES_SYSTEM), + ("property", ""), + ], + ) + def test_the_maps_name_their_code_system(self, name: str, codes: str) -> None: + assert ff.feature_input().CLASS_MAPS[name].codes == codes + + def test_a_coded_polygon_keeps_its_code_and_polygon(self, tmp_path: Path) -> None: + fs = self.open(tmp_path, [coded(1, WEST, "311")], "corine") + (feature,) = fs.features + assert feature.code == 311 + assert feature.polygon is not None + assert feature.polygon.equals(WEST) + + def test_a_coded_line_keeps_its_code_and_no_polygon(self, tmp_path: Path) -> None: + fs = self.open(tmp_path, [coded("r", ROAD, "122")], "corine") + (feature,) = fs.features + assert (feature.code, feature.polygon) == (122, None) + + def test_the_property_map_keeps_neither(self, tmp_path: Path) -> None: + fs = self.open(tmp_path, [Feat(1, WEST, {"property": "land_cover"})], "property") + (feature,) = fs.features + assert (feature.code, feature.polygon) == (None, None) + + def test_corine_water_keeps_the_water_only(self, tmp_path: Path) -> None: + fs = self.open( + tmp_path, [coded("f", FOREST, "312"), coded("l", LAKE, "512")], "corine-water" + ) + (feature,) = fs.features + assert (feature.fid, feature.code) == ("l", 512) + + def test_a_reprojected_polygon_is_moved_into_the_dems_crs(self, tmp_path: Path) -> None: + lonlat = ff.moved(WEST, "EPSG:25833", "EPSG:4326") + fs = self.open(tmp_path, [coded(1, lonlat, "311")], "corine", crs=None) + (feature,) = fs.features + expected = ff.moved(lonlat, "EPSG:4326") + assert feature.polygon.hausdorff_distance(expected) <= 1e-6 + assert feature.polygon.hausdorff_distance(WEST) <= 1e-3 + + def test_a_covering_polygon_is_kept_and_a_disjoint_one_dropped(self, tmp_path: Path) -> None: + far = rect(900.3, -900.1, 950.7, -850.2) + fs = self.open(tmp_path, [coded("c", COVER, "333"), coded("x", far, "311")], "corine") + (feature,) = fs.features + assert (feature.fid, feature.code, feature.lines) == ("c", 333, ()) + assert feature.polygon.equals(COVER) + assert (fs.outside, fs.clipped) == (1, 0) + + def test_an_uncoded_map_still_drops_a_covering_polygon(self, tmp_path: Path) -> None: + """16b's behaviour for `property`; the coded control keeps it.""" + uncoded = self.open(tmp_path, [Feat("c", COVER, {"property": "land_cover"})], "property") + assert (len(uncoded.features), uncoded.outside) == (0, 1) + kept = self.open(tmp_path, [coded("c", COVER, "333")], "corine") + assert len(kept.features) == 1 diff --git a/tests/python/test_io_ply.py b/tests/python/test_io_ply.py index 22745137..c09dcef6 100644 --- a/tests/python/test_io_ply.py +++ b/tests/python/test_io_ply.py @@ -440,3 +440,55 @@ def test_the_only_first_party_import_is_the_vocabulary(self) -> None: def test_it_is_deterministic(self) -> None: assert write_ply(VERTICES, faces=FACES) == write_ply(VERTICES, faces=FACES) + + +class TestFaceCodes: + """Increment 16c, R3: `write_ply(..., face_codes=...)` puts a face property + `int land_cover_code` after `vertex_indices`, on the face file only. + The header comment `land_cover_codes ...` is the caller's (it arrives in + `comments`), and `test_cli_mesh_landcover.py` pins it in the CLI's file. + + Committed red at `196147e`: `face_codes` was no parameter of `write_ply`, + so every call with it failed on `TypeError`. Green since `0487ed0`. + """ + + CODES = np.array([311, 2**31 - 1], dtype=np.int64) + + @pytest.mark.parametrize("ascii", [True, False], ids=["ascii", "binary"]) + def test_the_face_property_follows_the_index_list(self, ascii: bool) -> None: + header, _ = read_ply(write_ply(VERTICES, faces=FACES, face_codes=self.CODES, ascii=ascii)) + face = header.element("face") + assert [p.name for p in face.properties] == ["vertex_indices", "land_cover_code"] + assert face.properties[1].type_name == "int" + assert not face.properties[1].is_list + + @pytest.mark.parametrize("ascii", [True, False], ids=["ascii", "binary"]) + def test_the_codes_and_the_faces_round_trip(self, ascii: bool) -> None: + _, data = read_ply(write_ply(VERTICES, faces=FACES, face_codes=self.CODES, ascii=ascii)) + assert_array_equal(data["face"]["land_cover_code"], self.CODES) + assert_array_equal(np.stack(list(data["face"]["vertex_indices"])), FACES) + assert_array_equal(vertex_array(data), VERTICES) + + def test_ascii_and_binary_agree(self) -> None: + _, text = read_ply(write_ply(VERTICES, faces=FACES, face_codes=self.CODES)) + _, packed = read_ply(write_ply(VERTICES, faces=FACES, face_codes=self.CODES, ascii=False)) + assert_array_equal(text["face"]["land_cover_code"], packed["face"]["land_cover_code"]) + + def test_only_when_given(self) -> None: + coded, _ = read_ply(write_ply(VERTICES, faces=FACES, face_codes=self.CODES)) + plain, _ = read_ply(write_ply(VERTICES, faces=FACES)) + assert len(coded.element("face").properties) == 2 + assert [p.name for p in plain.element("face").properties] == ["vertex_indices"] + + def test_refused_with_edges(self) -> None: + with pytest.raises(ValueError, match="face_codes"): + write_ply(VERTICES, edges=EDGES, face_codes=self.CODES) + + @pytest.mark.parametrize("count", [len(FACES) - 1, len(FACES) + 1]) + def test_one_code_per_face(self, count: int) -> None: + with pytest.raises(ValueError, match="face"): + write_ply(VERTICES, faces=FACES, face_codes=np.zeros(count, dtype=np.int64)) + + def test_a_code_outside_int32_is_refused(self) -> None: + with pytest.raises(ValueError, match="int32"): + write_ply(VERTICES, faces=FACES, face_codes=np.array([311, 2**31], dtype=np.int64)) diff --git a/tests/python/test_io_vtk_legacy.py b/tests/python/test_io_vtk_legacy.py index 483261d4..d4d08398 100644 --- a/tests/python/test_io_vtk_legacy.py +++ b/tests/python/test_io_vtk_legacy.py @@ -518,3 +518,116 @@ def test_point_data_sits_between_the_cells_and_the_cell_data(self) -> None: def test_no_dataset_field_is_named_elevation(self) -> None: assert "elevation" not in read_vtk(write()).field_data + + +class TestLandCoverCode: + """Increment 16c, R3: `write_vtk(..., triangle_codes=...)` writes a cell + array `land_cover_code`, `int`, 0 on every `LINES` cell and each + triangle's code after them, in the cell `FIELD` block beside the + per-feature arrays (whose count includes it), or in a block of its own + when there are none. `feature_mask` stays the only `SCALARS` (13, ruling + 5). `land_cover_codes` joins `RESERVED`. + + Not pinned here: how the text of the `land_cover_codes` dataset string + reaches the writer. R3 reserves the name, so `fields` cannot carry it, + and the design names no parameter for it; `test_cli_mesh_landcover.py` + pins the string in the file. + + Committed red at `196147e`: `triangle_codes` was no parameter of + `write_vtk`, so every call with it failed on `TypeError`, and + `land_cover_codes` was not reserved. Green since `0487ed0`. + """ + + CODES = np.array([311, 0, 512, 2**31 - 1], dtype=np.int64) + + def coded(self, binary: bool) -> VtkFile: + return read_vtk(write(triangle_codes=self.CODES, binary=binary)) + + def test_lines_carry_0_then_each_triangle_its_code(self, binary: bool) -> None: + coded = self.coded(binary) + expected = np.concatenate([np.zeros(len(EDGES), dtype=np.int64), self.CODES]) + assert_array_equal(coded.cell_array("land_cover_code").values, expected) + + def test_it_is_an_int_scalar_covering_every_cell(self, binary: bool) -> None: + coded = self.coded(binary) + array = coded.cell_array("land_cover_code") + assert (array.type_name, array.components) == ("int", 1) + assert len(array.values) == len(EDGES) + len(TRIANGLES) + + def test_it_sits_in_the_features_block_beside_the_per_feature_arrays( + self, binary: bool + ) -> None: + coded = self.coded(binary) + assert list(coded.cell_fields) == ["features"] + assert set(coded.cell_fields["features"]) == {"river", "road", "railway", "land_cover_code"} + + def test_feature_mask_stays_the_only_scalar(self, binary: bool) -> None: + coded = self.coded(binary) + assert list(coded.scalars) == ["feature_mask"] + assert_array_equal(coded.cell_array("feature_mask").values, EXPECTED_MASK) + + def test_with_no_feature_arrays_it_is_written_alone(self, binary: bool) -> None: + parsed = read_vtk( + write( + edge_masks=np.zeros(len(EDGES), dtype=np.uint32), + triangle_codes=self.CODES, + binary=binary, + ) + ) + assert len(parsed.cell_fields) == 1 + (arrays,) = parsed.cell_fields.values() + assert list(arrays) == ["land_cover_code"] + + def test_with_no_edges_it_is_the_codes(self, binary: bool) -> None: + parsed = read_vtk( + write( + edges=np.zeros((0, 2), dtype=np.uint32), + edge_masks=np.zeros(0, dtype=np.uint32), + triangle_codes=self.CODES, + binary=binary, + ) + ) + assert_array_equal(parsed.cell_array("land_cover_code").values, self.CODES) + + def test_only_when_given(self) -> None: + with_codes = read_vtk(write(triangle_codes=self.CODES)) + without = read_vtk(write()) + assert "land_cover_code" in with_codes.cell_fields["features"] + with pytest.raises(KeyError): + without.cell_array("land_cover_code") + + def test_ascii_and_binary_agree(self) -> None: + text = read_vtk(write(triangle_codes=self.CODES)) + packed = read_vtk(write(triangle_codes=self.CODES, binary=True)) + assert_array_equal( + text.cell_array("land_cover_code").values, packed.cell_array("land_cover_code").values + ) + + def test_binary_is_big_endian_int32(self) -> None: + blob = write( + edges=np.zeros((0, 2), dtype=np.uint32), + edge_masks=np.zeros(0, dtype=np.uint32), + triangle_codes=self.CODES, + binary=True, + ) + assert struct.pack(">4i", *self.CODES.tolist()) in blob + + @pytest.mark.parametrize("count", [len(TRIANGLES) - 1, len(TRIANGLES) + 1]) + def test_one_code_per_triangle(self, count: int) -> None: + with pytest.raises(ValueError, match="triangle"): + write(triangle_codes=np.zeros(count, dtype=np.int64)) + + def test_a_code_outside_int32_is_refused(self) -> None: + with pytest.raises(ValueError, match="int32"): + write(triangle_codes=np.array([311, 0, 512, 2**31], dtype=np.int64)) + + def test_a_vocabulary_naming_land_cover_code_is_refused(self) -> None: + clash = EdgeVocabulary(properties=(EdgeProperty(name="land_cover_code", bit=0),)) + masks = np.array([1, 0, 0], dtype=np.uint32) + write(vocabulary=clash, edge_masks=masks) # without codes: an ordinary name + with pytest.raises(ValueError, match="land_cover_code"): + write(vocabulary=clash, edge_masks=masks, triangle_codes=self.CODES) + + def test_land_cover_codes_is_reserved(self) -> None: + with pytest.raises(ValueError, match="land_cover_codes"): + write(fields=(("land_cover_codes", "x"),)) diff --git a/tests/python/test_landcover.py b/tests/python/test_landcover.py new file mode 100644 index 00000000..5dc8811e --- /dev/null +++ b/tests/python/test_landcover.py @@ -0,0 +1,453 @@ +"""`tin_engine.landcover`: a land-cover class per triangle (increment 16c, R1). + +`docs/increments/16c-landcover-labels.md`, "Tests for @tester": `regions` and +`label_triangles` on hand-made meshes, no `_core` call, no file. Every mesh +here is a rectilinear grid whose lines carry every input boundary, so the +constraint edges are exactly the grid edges lying on the input linework (the +domain's boundary, each polygon's rings, each road), found by distance and not +by the producer. Coordinates sit at UTM 33N magnitudes, off round numbers. + +Each labelling fixture is checked three ways: against the centroid oracle of +`landcover_fixtures` (I2), for I1 (the spread), and against the codes the +design names for it. The fixtures are the design's CLI list as pure meshes: +two squares side by side, a forest with an empty hole, a lake in a holed and +in an unholed forest, a polygon clipped by the domain into two pieces, a road +across a polygon, a polygon covering the domain, boundaries 0.1 mm apart, +determinism, and overlaps (smallest area, ties to the smaller code; D2). + +Pinned beyond the design's text (the design leaves them open): + +- `label_triangles(vertices, triangles, edges, polygons=..., margin=...)`, + `vertices` `(N, 3)` (the trimmed mesh's; only x and y are used), and + `polygons` a sequence of `(geometry, code)` pairs. +- `CoverLabels.regions`, `.outside`, `.overlapped` and `.thin` are counts + (ints), the numbers of R3's stderr line; `.codes` is the `(T,)` array. +- `regions(triangles, edges)` takes no vertex count. + +Committed red at `196147e`: `tin_engine.landcover` did not exist yet, so +every test failed on `ModuleNotFoundError`. The module landed in `0487ed0` +and the suite has been green since. +""" + +from __future__ import annotations + +from collections.abc import Sequence +from dataclasses import dataclass +from typing import Any + +import numpy as np +import pytest +import shapely +from numpy.testing import assert_array_equal +from shapely.geometry import LineString, MultiLineString, Polygon, box +from shapely.geometry.base import BaseGeometry + +from importscan import first_party_imports +from landcover_fixtures import ( + MARGIN, + edge_keys, + interior_edge_count, + landcover_oracle, + spread_violations, +) +from tin_engine import landcover + +X0, Y0 = 500_000.3, 6_600_000.7 + + +def rect(x0: float, y0: float, x1: float, y1: float) -> Polygon: + return box(X0 + x0, Y0 + y0, X0 + x1, Y0 + y1) + + +def line(*points: tuple[float, float]) -> LineString: + return LineString([(X0 + x, Y0 + y) for x, y in points]) + + +# ---------------------------------------------------------------- the meshes + + +@dataclass(frozen=True) +class Mesh: + """A trimmed mesh as the writers see it: `(N, 3)` vertices, triangles, and + the constraint edges.""" + + vertices: np.ndarray + triangles: np.ndarray + edges: np.ndarray + + @property + def centroids(self) -> np.ndarray: + return self.vertices[self.triangles][:, :, :2].mean(axis=1) + + def inside(self, geometry: BaseGeometry) -> np.ndarray: + """Triangles whose centroid is strictly inside `geometry`.""" + return np.asarray(shapely.contains_xy(geometry, *self.centroids.T)) + + +def grid_mesh( + xs: Sequence[float], + ys: Sequence[float], + lines: Sequence[BaseGeometry], + domain: Polygon | None = None, +) -> Mesh: + """A rectilinear grid, each cell split on its rising diagonal, keeping the + cells whose centre is in `domain` (default: the whole grid). The + constraint edges are the grid edges lying on `lines` or on the domain's + boundary, within a micrometre at both ends and the midpoint.""" + gx, gy = np.meshgrid(np.asarray(xs) + X0, np.asarray(ys) + Y0, indexing="xy") + vertices = np.column_stack([gx.ravel(), gy.ravel(), np.zeros(gx.size)]) + nx = len(xs) + if domain is None: + domain = box(X0 + xs[0], Y0 + ys[0], X0 + xs[-1], Y0 + ys[-1]) + triangles = [] + for j in range(len(ys) - 1): + for i in range(nx - 1): + v00, v10 = j * nx + i, j * nx + i + 1 + v01, v11 = v00 + nx, v10 + nx + centre = vertices[[v00, v11], :2].mean(axis=0) + if domain.contains(shapely.Point(*centre)): + triangles += [(v00, v10, v11), (v00, v11, v01)] + tri = np.array(triangles, dtype=np.int64) + sides = np.concatenate([tri[:, [0, 1]], tri[:, [1, 2]], tri[:, [2, 0]]]) + _, first = np.unique(edge_keys(sides, len(vertices)), return_index=True) + candidates = sides[np.sort(first)] + linework = shapely.union_all([domain.boundary, *(_linework(g) for g in lines)]) + a, b = vertices[candidates[:, 0], :2], vertices[candidates[:, 1], :2] + near = [ + np.asarray(shapely.distance(linework, shapely.points(p))) <= 1e-6 + for p in (a, b, (a + b) / 2) + ] + on = near[0] & near[1] & near[2] + return Mesh(vertices=vertices, triangles=tri, edges=candidates[on]) + + +def _linework(geometry: BaseGeometry) -> BaseGeometry: + if isinstance(geometry, Polygon): + return MultiLineString([r.coords for r in (geometry.exterior, *geometry.interiors)]) + return geometry + + +def label(mesh: Mesh, polygons: Sequence[tuple[BaseGeometry, int]], margin: float = MARGIN) -> Any: + return landcover.label_triangles( + mesh.vertices, mesh.triangles, mesh.edges, polygons=list(polygons), margin=margin + ) + + +def checked_label( + mesh: Mesh, polygons: Sequence[tuple[BaseGeometry, int]], margin: float = MARGIN +) -> Any: + """`label_triangles`, with I1 and I2 asserted on the way out.""" + labels = label(mesh, polygons, margin) + codes = np.asarray(labels.codes) + assert codes.shape == (len(mesh.triangles),) + assert interior_edge_count(mesh.triangles) > 0 + assert spread_violations(mesh.triangles, mesh.edges, codes) == [] + oracle = landcover_oracle(mesh.vertices, mesh.triangles, polygons, margin) + assert oracle.checked.any() + assert oracle.mismatches(codes) == [] + return labels + + +STEPS = [float(v) for v in range(21)] +DOMAIN = rect(0, 0, 20, 20) + + +# ------------------------------------------------------------------- regions + + +def two_triangles() -> tuple[np.ndarray, np.ndarray]: + return np.array([[0, 1, 2], [0, 2, 3]]), np.array([[0, 1], [1, 2], [2, 3], [3, 0]]) + + +def strip(cells: int) -> np.ndarray: + """`2 * cells` triangles in a row: bottom row 0..cells, top row after it.""" + top = cells + 1 + out = [] + for i in range(cells): + out += [(i, i + 1, top + i + 1), (i, top + i + 1, top + i)] + return np.array(out, dtype=np.int64) + + +class TestRegions: + def test_an_unconstrained_shared_edge_joins_two_triangles(self) -> None: + triangles, hull = two_triangles() + assert_array_equal(landcover.regions(triangles, hull), [0, 0]) + + def test_a_constrained_shared_edge_separates_them(self) -> None: + triangles, hull = two_triangles() + cut = np.vstack([hull, [[0, 2]]]) + assert_array_equal(landcover.regions(triangles, cut), [0, 1]) + + def test_a_constraint_given_backwards_still_blocks(self) -> None: + triangles, hull = two_triangles() + cut = np.vstack([hull, [[2, 0]]]) + assert_array_equal(landcover.regions(triangles, cut), [0, 1]) + + def test_no_constraints_at_all(self) -> None: + triangles, _ = two_triangles() + assert_array_equal(landcover.regions(triangles, np.zeros((0, 2), dtype=np.int64)), [0, 0]) + + @pytest.mark.parametrize("order", ["forward", "reversed", "shuffled"]) + def test_a_strip_of_1000_triangles_is_one_component(self, order: str) -> None: + """The union-find loops to a fixed point: a chain 1 000 long is one + component whatever order its triangles are listed in.""" + triangles = strip(500) + assert len(triangles) == 1000 + if order == "reversed": + triangles = triangles[::-1] + elif order == "shuffled": + triangles = triangles[np.random.default_rng(16).permutation(len(triangles))] + ids = np.asarray(landcover.regions(triangles, np.zeros((0, 2), dtype=np.int64))) + assert ids.shape == (1000,) + assert set(ids.tolist()) == {0} + + def test_ids_are_the_smallest_triangle_index_of_each_component(self) -> None: + cut = line((10, 0), (10, 20)) + mesh = grid_mesh(STEPS, STEPS, [cut]) + order = np.random.default_rng(3).permutation(len(mesh.triangles)) + triangles = mesh.triangles[order] + ids = np.asarray(landcover.regions(triangles, mesh.edges)) + west = np.asarray(shapely.contains_xy(rect(0, 0, 10, 20), *mesh.centroids[order].T)) + assert west.any() and (~west).any() + assert set(ids[west].tolist()) == {int(np.flatnonzero(west).min())} + assert set(ids[~west].tolist()) == {int(np.flatnonzero(~west).min())} + + def test_permuting_and_flipping_the_constraint_rows_changes_nothing(self) -> None: + mesh = grid_mesh(STEPS, STEPS, [rect(4, 4, 16, 16), line((0, 10), (20, 10))]) + rng = np.random.default_rng(7) + shuffled = mesh.edges[rng.permutation(len(mesh.edges))] + flipped = shuffled[:, ::-1] + first = np.asarray(landcover.regions(mesh.triangles, mesh.edges)) + assert len(set(first.tolist())) == 4 + assert_array_equal(landcover.regions(mesh.triangles, shuffled), first) + assert_array_equal(landcover.regions(mesh.triangles, flipped), first) + + +# ---------------------------------------------------------- label_triangles + + +class TestLabelBasics: + def test_codes_are_int32_one_per_triangle(self) -> None: + mesh = grid_mesh(STEPS, STEPS, [rect(4, 4, 16, 16)]) + labels = label(mesh, [(rect(4, 4, 16, 16), 311)]) + codes = np.asarray(labels.codes) + assert codes.dtype == np.int32 + assert codes.shape == (len(mesh.triangles),) + + def test_a_square_cut_by_a_constrained_diagonal(self) -> None: + """Each side of the diagonal gets its own polygon's code.""" + diagonal = line((0, 0), (20, 20)) + lower = Polygon([(X0, Y0), (X0 + 20, Y0), (X0 + 20, Y0 + 20)]) + upper = Polygon([(X0, Y0), (X0 + 20, Y0 + 20), (X0, Y0 + 20)]) + mesh = grid_mesh(STEPS, STEPS, [diagonal]) + labels = checked_label(mesh, [(lower, 311), (upper, 512)]) + codes = np.asarray(labels.codes) + assert set(codes[mesh.inside(lower)].tolist()) == {311} + assert set(codes[mesh.inside(upper)].tolist()) == {512} + assert (labels.regions, labels.outside, labels.overlapped) == (2, 0, 0) + + def test_no_polygons_is_all_zero_and_every_region_outside(self) -> None: + mesh = grid_mesh(STEPS, STEPS, [rect(4, 4, 16, 16), line((0, 10), (20, 10))]) + labels = checked_label(mesh, []) + assert set(np.asarray(labels.codes).tolist()) == {0} + assert labels.regions == 4 + assert labels.outside == labels.regions + assert labels.overlapped == 0 + + def test_a_sliver_component_is_counted_thin_and_still_labelled(self) -> None: + """A 1 mm sliver beside a big triangle, the shared edge constrained: + two components, the sliver's largest inradius about 0.5 mm <= margin.""" + vertices = np.array( + [ + [X0, Y0, 0.0], + [X0 + 100, Y0, 0.0], + [X0 + 100, Y0 + 100, 0.0], + [X0 + 100.001, Y0 + 50, 0.0], + ] + ) + triangles = np.array([[0, 1, 2], [1, 3, 2]]) + edges = np.array([[0, 1], [1, 3], [3, 2], [2, 0], [1, 2]]) + cover = rect(-10, -10, 110, 110) + labels = landcover.label_triangles( + vertices, triangles, edges, polygons=[(cover, 311)], margin=MARGIN + ) + assert_array_equal(labels.codes, [311, 311]) + assert (labels.regions, labels.thin) == (2, 1) + finer = landcover.label_triangles( + vertices, triangles, edges, polygons=[(cover, 311)], margin=1e-5 + ) + assert finer.thin == 0 + + +# ------------------------------------------------------------ the fixtures + + +class TestFixtures: + def test_two_squares_side_by_side(self) -> None: + west, east = rect(2, 4, 10, 16), rect(10, 4, 18, 16) + mesh = grid_mesh(STEPS, STEPS, [west, east]) + labels = checked_label(mesh, [(west, 311), (east, 512)]) + codes = np.asarray(labels.codes) + assert set(codes[mesh.inside(west)].tolist()) == {311} + assert set(codes[mesh.inside(east)].tolist()) == {512} + assert set(codes[~mesh.inside(west.union(east))].tolist()) == {0} + # The triangles touching the shared side carry both codes. + seam = line((10, 5), (10, 15)) + touching = ( + np.asarray( + shapely.distance(seam, shapely.polygons(mesh.vertices[mesh.triangles][:, :, :2])) + ) + == 0 + ) + assert set(codes[touching].tolist()) == {311, 512} + assert (labels.regions, labels.outside, labels.overlapped, labels.thin) == (3, 1, 0, 0) + + def test_a_forest_with_an_empty_hole(self) -> None: + hole = rect(8, 8, 12, 12) + forest = Polygon(rect(2, 2, 18, 18).exterior, [hole.exterior]) + mesh = grid_mesh(STEPS, STEPS, [forest]) + labels = checked_label(mesh, [(forest, 312)]) + codes = np.asarray(labels.codes) + assert set(codes[mesh.inside(hole)].tolist()) == {0} + assert set(codes[mesh.inside(forest)].tolist()) == {312} + assert (labels.regions, labels.outside) == (3, 2) + + def test_a_lake_filling_the_hole_of_a_holed_forest(self) -> None: + lake = rect(8, 8, 12, 12) + forest = Polygon(rect(2, 2, 18, 18).exterior, [lake.exterior]) + mesh = grid_mesh(STEPS, STEPS, [forest, lake]) + labels = checked_label(mesh, [(forest, 312), (lake, 512)]) + codes = np.asarray(labels.codes) + assert set(codes[mesh.inside(lake)].tolist()) == {512} + assert set(codes[mesh.inside(forest)].tolist()) == {312} + assert labels.overlapped == 0 + + def test_a_lake_inside_an_unholed_forest_goes_to_the_smaller(self) -> None: + """D2: the lake's component is in both polygons; the lake is smaller.""" + lake = rect(8, 8, 12, 12) + forest = rect(2, 2, 18, 18) + mesh = grid_mesh(STEPS, STEPS, [forest, lake]) + labels = checked_label(mesh, [(forest, 312), (lake, 512)]) + codes = np.asarray(labels.codes) + assert set(codes[mesh.inside(lake)].tolist()) == {512} + assert set(codes[mesh.inside(forest.difference(lake))].tolist()) == {312} + assert labels.overlapped == 1 + + def test_a_polygon_clipped_by_the_domain_into_two_pieces(self) -> None: + """A band crossing a notch cut into the domain from the north: its two + pieces are separate components, and both carry its code.""" + domain = rect(0, 0, 20, 20).difference(rect(8, 10, 12, 20)) + band = rect(4, 12, 16, 16) + mesh = grid_mesh(STEPS, STEPS, [band], domain=domain) + labels = checked_label(mesh, [(band, 324)]) + codes = np.asarray(labels.codes) + west, east = mesh.inside(rect(4, 12, 8, 16)), mesh.inside(rect(12, 12, 16, 16)) + assert west.any() and east.any() + assert set(codes[west | east].tolist()) == {324} + ids = np.asarray(landcover.regions(mesh.triangles, mesh.edges)) + assert set(ids[west].tolist()).isdisjoint(ids[east].tolist()) + assert (labels.regions, labels.outside) == (3, 1) + + def test_a_road_across_a_polygon(self) -> None: + """The road is a constraint and no polygon: both halves are the forest's.""" + forest = rect(4, 4, 16, 16) + road = line((2, 10), (18, 10)) + mesh = grid_mesh(STEPS, STEPS, [forest, road]) + labels = checked_label(mesh, [(forest, 311)]) + codes = np.asarray(labels.codes) + north, south = mesh.inside(rect(4, 10, 16, 16)), mesh.inside(rect(4, 4, 16, 10)) + assert set(codes[north | south].tolist()) == {311} + ids = np.asarray(landcover.regions(mesh.triangles, mesh.edges)) + assert set(ids[north].tolist()).isdisjoint(ids[south].tolist()) + assert (labels.regions, labels.outside) == (3, 1) + + def test_a_polygon_covering_the_domain(self) -> None: + cover = rect(-5, -5, 25, 25) + mesh = grid_mesh(STEPS, STEPS, []) + labels = checked_label(mesh, [(cover, 333)]) + assert set(np.asarray(labels.codes).tolist()) == {333} + assert (labels.regions, labels.outside, labels.overlapped) == (1, 0, 0) + + def test_boundaries_a_tenth_of_a_millimetre_apart(self) -> None: + """Two squares whose shared side is offset by 0.1 mm, less than the + snap, with nothing merged (the pure form of the CLI fixture, where + the noder may merge them): the sliver column between the two sides + opens onto the outside at both ends, so it is part of the outside + component; away from the seam each square has its own code (I2).""" + xs = [*STEPS[:11], 10.0001, *STEPS[11:]] + west, east = rect(2, 4, 10, 16), rect(10.0001, 4, 18, 16) + mesh = grid_mesh(xs, STEPS, [west, east]) + labels = checked_label(mesh, [(west, 311), (east, 512)]) + codes = np.asarray(labels.codes) + assert set(codes[mesh.inside(rect(2, 4, 9.99, 16))].tolist()) == {311} + assert set(codes[mesh.inside(rect(10.01, 4, 18, 16))].tolist()) == {512} + assert (labels.regions, labels.outside, labels.thin) == (3, 1, 0) + gap = mesh.inside(rect(10, 4.5, 10.0001, 15.5)) + assert gap.any() and set(codes[gap].tolist()) == {0} + + +class TestOverlaps: + """D2: the smallest area wins, ties to the smaller code, in any order.""" + + @pytest.mark.parametrize("reverse", [False, True]) + def test_equal_areas_go_to_the_smaller_code(self, reverse: bool) -> None: + square = rect(4, 4, 16, 16) + mesh = grid_mesh(STEPS, STEPS, [square]) + polygons = [(square, 412), (rect(4, 4, 16, 16), 322)] + labels = checked_label(mesh, polygons[::-1] if reverse else polygons) + assert set(np.asarray(labels.codes)[mesh.inside(square)].tolist()) == {322} + assert labels.overlapped == 1 + + @pytest.mark.parametrize("reverse", [False, True]) + def test_nested_polygons_give_the_inner_code(self, reverse: bool) -> None: + inner, outer = rect(6, 6, 14, 14), rect(-1, -1, 21, 21) + mesh = grid_mesh(STEPS, STEPS, [inner]) + polygons = [(outer, 312), (inner, 512)] + labels = checked_label(mesh, polygons[::-1] if reverse else polygons) + codes = np.asarray(labels.codes) + assert set(codes[mesh.inside(inner)].tolist()) == {512} + assert set(codes[~mesh.inside(inner)].tolist()) == {312} + assert (labels.regions, labels.overlapped, labels.outside) == (2, 1, 0) + + +class TestDeterminism: + """R1: the codes are a function of the triangles, the constraint edges and + the polygons as a set.""" + + @pytest.fixture + def case(self) -> tuple[Mesh, list[tuple[BaseGeometry, int]]]: + lake = rect(8, 8, 12, 12) + forest = rect(2, 2, 18, 18) + bog = rect(12, 2, 18, 6) + road = line((0, 15), (20, 15)) + mesh = grid_mesh(STEPS, STEPS, [forest, lake, bog, road]) + return mesh, [(forest, 312), (lake, 512), (bog, 412)] + + def test_the_same_call_twice_is_the_same_array(self, case: Any) -> None: + mesh, polygons = case + first, second = label(mesh, polygons), label(mesh, polygons) + assert np.asarray(first.codes).tobytes() == np.asarray(second.codes).tobytes() + + def test_permuted_triangles_edges_and_polygons_give_the_same_labels(self, case: Any) -> None: + mesh, polygons = case + base = np.asarray(label(mesh, polygons).codes) + rng = np.random.default_rng(29) + order = rng.permutation(len(mesh.triangles)) + rows = rng.permutation(len(mesh.edges)) + shuffled = Mesh(mesh.vertices, mesh.triangles[order], mesh.edges[rows][:, ::-1]) + again = label(shuffled, polygons[::-1]) + assert_array_equal(np.asarray(again.codes), base[order]) + first = label(mesh, polygons) + assert (again.regions, again.outside, again.overlapped, again.thin) == ( + first.regions, + first.outside, + first.overlapped, + first.thin, + ) + + +class TestPurity: + def test_it_imports_nothing_first_party(self) -> None: + """R1's boundary: numpy and shapely, nothing first-party, never `_core`.""" + module = landcover + assert first_party_imports(module) == set() diff --git a/tests/python/test_palettes.py b/tests/python/test_palettes.py new file mode 100644 index 00000000..8e64aca6 --- /dev/null +++ b/tests/python/test_palettes.py @@ -0,0 +1,221 @@ +"""`tin_engine.palettes` and `rasputin palette`: natural colours for CORINE (16c, R4). + +`docs/increments/16c-landcover-labels.md`, R4 and "Tests for @tester". The +colours are Ola's to change (D6), so they are pinned by family, as hue, +saturation and lightness ranges, never as exact triples: forest green, heath a +lighter green than conifers, bare rock grey, peat bog brown, every water class +blue, glaciers near white and not warm. + +Pinned beyond the design's text: + +- An `Annotations` value is compared as `int(value)`, so the preset may carry + codes as strings (ParaView's own presets do) or as numbers. +- Each annotation's label contains the class's `LABEL3` name. +- `NanColor` is present and magenta (R4: "A code the table lacks draws in + `NanColor`, magenta"). + +Committed red at `196147e`: `tin_engine.palettes` did not exist yet and +`rasputin palette` was no command, so every test failed on +`ModuleNotFoundError` or on the CLI's usage error for an unknown command. Both +landed in `0487ed0` and the suite has been green since. +""" + +from __future__ import annotations + +import colorsys +import json +import re +import sqlite3 +from contextlib import closing +from pathlib import Path +from typing import Any + +import pytest +from typer.testing import CliRunner + +from gpkg_fixtures import OLA_NORWAY +from tin_engine import palettes +from tin_engine.cli import app + +runner = CliRunner(env={"NO_COLOR": "1", "TERM": "dumb"}) + +#: R4's "in the Norway extract" column: 34 codes. +NORWAY = frozenset( + { + *(111, 112, 121, 122, 123, 124, 131, 132, 133, 141, 142), + *(211, 222, 231, 242, 243), + *(311, 312, 313, 321, 322, 324, 331, 332, 333, 334, 335), + *(411, 412, 423), + *(511, 512, 522, 523), + } +) + +#: The 44 CLC level-3 classes (Kosztra et al. 2019). +CLC = frozenset( + { + *(111, 112, 121, 122, 123, 124, 131, 132, 133, 141, 142), + *(211, 212, 213, 221, 222, 223, 231, 241, 242, 243, 244), + *(311, 312, 313, 321, 322, 323, 324, 331, 332, 333, 334, 335), + *(411, 412, 421, 422, 423), + *(511, 512, 521, 522, 523), + } +) + +PRESET_NAME = "rasputin CORINE natural" +HEX = re.compile(r"#[0-9a-fA-F]{6}") + + +def table() -> dict[int, tuple[str, str]]: + return dict(palettes.CORINE_NATURAL) + + +def rgb(code: int) -> tuple[float, float, float]: + """The table's colour for `code`, as r, g, b in 0 .. 1.""" + text = table()[code][1] + assert HEX.fullmatch(text), text + r, g, b = (int(text[i : i + 2], 16) / 255 for i in (1, 3, 5)) + return r, g, b + + +def hls(code: int) -> tuple[float, float, float]: + """Hue in degrees, lightness and saturation (HLS), each 0 .. 1 but hue.""" + h, lightness, s = colorsys.rgb_to_hls(*rgb(code)) + return h * 360, lightness, s + + +# --------------------------------------------------------------- the table + + +class TestTheTable: + def test_every_clc_class_and_zero_has_a_colour(self) -> None: + assert set(table()) == CLC | {0} + + def test_every_norway_code_has_a_colour(self) -> None: + assert len(NORWAY) == 34 and NORWAY <= CLC and len(CLC) == 44 + assert NORWAY - set(table()) == set() + + def test_entries_are_a_name_and_a_hex_colour(self) -> None: + for code, (name, colour) in table().items(): + assert isinstance(name, str) and name, code + assert HEX.fullmatch(colour), (code, colour) + + def test_the_names_are_label3(self) -> None: + # Spot checks from `Legend/CLC_legend.csv` as R4 quotes it. + assert table()[312][0] == "Coniferous forest" + assert table()[412][0] == "Peat bogs" + assert table()[335][0] == "Glaciers and perpetual snow" + assert table()[512][0] == "Water bodies" + + def test_the_colours_are_distinct(self) -> None: + colours = [c.lower() for _, c in table().values()] + assert len(set(colours)) == len(colours) + + +class TestFamilies: + """D6: the families, by hue ranges, not by exact triple.""" + + @pytest.mark.parametrize("code", [311, 312, 313]) + def test_forests_are_green(self, code: int) -> None: + hue, _, saturation = hls(code) + assert 75 <= hue <= 165 and saturation >= 0.2, hls(code) + + def test_heath_is_green_and_lighter_than_conifers(self) -> None: + hue, lightness, saturation = hls(322) + assert 60 <= hue <= 150 and saturation >= 0.15, hls(322) + assert lightness > hls(312)[1] + + def test_bare_rock_is_grey(self) -> None: + _, lightness, saturation = hls(332) + assert saturation <= 0.08 and 0.25 <= lightness <= 0.75, hls(332) + + def test_peat_bogs_are_brown(self) -> None: + hue, lightness, saturation = hls(412) + assert 15 <= hue <= 50 and saturation >= 0.2 and lightness <= 0.55, hls(412) + + @pytest.mark.parametrize("code", sorted(c for c in CLC if c // 100 == 5)) + def test_water_is_blue(self, code: int) -> None: + hue, _, saturation = hls(code) + assert 190 <= hue <= 235 and saturation >= 0.25, hls(code) + + def test_glaciers_are_near_white_and_not_warm(self) -> None: + r, _, b = rgb(335) + _, lightness, _ = hls(335) + assert lightness >= 0.9 and b >= r, rgb(335) + + def test_the_official_legend_is_not_used_for_bogs_or_the_sea(self) -> None: + """Direction 2: the official legend paints peat bogs blue (077-077-255) + and the sea almost white (230-242-255).""" + assert table()[412][1].lower() != "#4d4dff" + assert table()[523][1].lower() != "#e6f2ff" + + +# ---------------------------------------------------------- the preset + + +class TestPreset: + def preset(self) -> dict[str, Any]: + out = palettes.paraview_preset(palettes.CORINE_NATURAL, "a name") + assert isinstance(out, list) and len(out) == 1 + assert isinstance(out[0], dict) + return out[0] + + def test_it_carries_the_name(self) -> None: + preset = self.preset() + assert preset["Name"] == "a name" + + def test_annotations_and_colours_have_one_entry_per_code(self) -> None: + preset = self.preset() + entries = len(table()) + assert len(preset["Annotations"]) == 2 * entries + assert len(preset["IndexedColors"]) == 3 * entries + assert all(0.0 <= float(c) <= 1.0 for c in preset["IndexedColors"]) + + def test_they_are_in_the_same_order(self) -> None: + preset = self.preset() + values = [int(v) for v in preset["Annotations"][0::2]] + labels = [str(v) for v in preset["Annotations"][1::2]] + colours = preset["IndexedColors"] + assert sorted(values) == sorted(table()) + for k, code in enumerate(values): + assert table()[code][0] in labels[k], (code, labels[k]) + got = [float(c) for c in colours[3 * k : 3 * k + 3]] + assert got == pytest.approx(rgb(code), abs=1e-6), code + + def test_the_nan_colour_is_magenta(self) -> None: + preset = self.preset() + r, g, b = (float(c) for c in preset["NanColor"]) + assert r >= 0.8 and b >= 0.8 and g <= 0.3, preset["NanColor"] + + def test_it_is_json(self) -> None: + preset = self.preset() + assert json.loads(json.dumps(preset)) == preset + + +# ------------------------------------------------------------- the command + + +class TestCommand: + def test_stdout_is_the_preset(self) -> None: + result = runner.invoke(app, ["palette", "corine"]) + assert result.exit_code == 0, result.output + expected = palettes.paraview_preset(palettes.CORINE_NATURAL, PRESET_NAME) + assert json.loads(result.stdout) == expected + + def test_an_unknown_name_is_a_usage_error_listing_the_known(self) -> None: + result = runner.invoke(app, ["palette", "nosuch"]) + assert result.exit_code == 2, result.output + assert "nosuch" in result.output and "corine" in result.output + assert "No such command" not in result.output + + +# --------------------------------------------------------- Ola's local data + + +@pytest.mark.skipif(not OLA_NORWAY.exists(), reason=f"{OLA_NORWAY} is not here") +def test_every_code_in_olas_norway_extract_has_a_colour() -> None: + """R4's query, rerun: every `code_18` in Ola's Norway GeoPackage.""" + uri = Path(OLA_NORWAY).resolve().as_uri() + "?mode=ro" + with closing(sqlite3.connect(uri, uri=True)) as con: + codes = {int(c) for (c,) in con.execute("SELECT DISTINCT code_18 FROM corine2018")} + assert codes == NORWAY + assert codes <= set(table())