22 (2/2): auto-catchment, the outline reduction (stacked on #110)

ROADMAP.md

@@ -49,7 +49,7 @@ increment that most needs a picture to check against
 | 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; **a requirement, not an option (Ola, 2026-09-29): "we will get an extreme amount of points in the catchment polygon. This must be taken into account, or the number of triangles will be overwhelming."** A cell outline has a vertex at every cell step (tens of thousands for a mid-sized catchment at 10 m, far more for Glomma), and every boundary edge is a constraint the mesh must honour; 16b measured what dense constraints cost (CORINE borders make the 10 m mesh 3.6x larger). The design must state how the outline is reduced, to what horizontal tolerance, and what it guarantees (a simple polygon, the pour point inside); the fine outline is defined on the DEM's node lattice (Ola, 2026-09-29: "Isn't the DEM points, really?"), but the reduced polygon's vertices need not be DEM nodes: off-node domain vertices are supported since 16 (Ola: "We have aleady established that the polygons and DEM don't need to match"), so vertices may be placed to keep the area; and the reduced polygon keeps approximately the fine one's area (Ola, 2026-09-29: "the resulting polygon has approximately the same area as the fine one"), since catchment area drives runoff volume. Literature to check: area-preserving polyline simplification (e.g. Bose et al. 2006; Kronenfeld et al. 2020, segment collapse); 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) | **increment 22, PR 1 built** (2026-09-29, overnight): the fine catchment, net 619 lines. Bygdin: 304.91 km² against NVE's 305.54 km² (−0.21 %), in 6.5 s. PR 2, the outline reduction, follows as a stacked PR. Designed by `@architect`, 2026-09-29, Bygdin first, moved ahead of 20c for Ola's night run: a lake polygon (CORINE) as the seed, a Priority-Flood labelling the lake's catchment, a marching-squares outline, area-preserving segment collapse to `--outline-tolerance` (default twice the cell). Two PRs on one branch; defaults for Ola to confirm are marked in the file | `docs/increments/22-auto-catchment.md` |
+| — | 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; **a requirement, not an option (Ola, 2026-09-29): "we will get an extreme amount of points in the catchment polygon. This must be taken into account, or the number of triangles will be overwhelming."** A cell outline has a vertex at every cell step (tens of thousands for a mid-sized catchment at 10 m, far more for Glomma), and every boundary edge is a constraint the mesh must honour; 16b measured what dense constraints cost (CORINE borders make the 10 m mesh 3.6x larger). The design must state how the outline is reduced, to what horizontal tolerance, and what it guarantees (a simple polygon, the pour point inside); the fine outline is defined on the DEM's node lattice (Ola, 2026-09-29: "Isn't the DEM points, really?"), but the reduced polygon's vertices need not be DEM nodes: off-node domain vertices are supported since 16 (Ola: "We have aleady established that the polygons and DEM don't need to match"), so vertices may be placed to keep the area; and the reduced polygon keeps approximately the fine one's area (Ola, 2026-09-29: "the resulting polygon has approximately the same area as the fine one"), since catchment area drives runoff volume. Literature to check: area-preserving polyline simplification (e.g. Bose et al. 2006; Kronenfeld et al. 2020, segment collapse); 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) | **built as increment 22** (2026-09-29, overnight), in two PRs: PR 1, the fine catchment, net 619 lines; PR 2, the outline reduction, net 406 lines. Bygdin: 304.91 km² against NVE's 305.54 km² (−0.21 %), 99.1/99.3 % node overlap with NVE's polygon, 17 812 → 740 outline vertices at 20 m with the area kept, catchment in 6.5 s (`docs/benchmarks/2026-09-29/bygdin.md`). The reduction bounds the distance to the fine outline at its vertices only (measured on Bygdin: 19.92 m at 20 m). Parked for later (Ola, 2026-09-29): "It could be useful to have the option of some tolerance, doing more adaptive coarsening, but write it down and leave it for now": an enforced whole-outline tolerance as an option, which also belongs with 20c's input coarsening. Designed by `@architect`, 2026-09-29, Bygdin first, moved ahead of 20c for Ola's night run: a lake polygon (CORINE) as the seed, a Priority-Flood labelling the lake's catchment, a marching-squares outline, area-preserving segment collapse to `--outline-tolerance` (default twice the cell). Two PRs on one branch; defaults for Ola to confirm are marked in the file | `docs/increments/22-auto-catchment.md` |
 | — | `raster/`: grid-to-world geometry and bilinear sampling | shipped (`7785fea`), **no record** | none — predates the protocol |
 
 ## What stands between here and an operational MVP

bindings/core.cpp

@@ -20,6 +20,7 @@
 #include <terrain/raster/sample.hpp>
 #include <terrain/raster/view.hpp>
 #include <terrain/refinement/refine.hpp>
+#include <terrain/vector_simplify/area_collapse.hpp>
 
 #include <cstddef>
 #include <cstdint>
@@ -299,6 +300,17 @@ struct BoundUpstream {
     std::size_t cols;
 };
 
+// An (N, 2) float64 array as Point2 pairs, copied; any other shape is a ValueError.
+[[nodiscard]] std::vector<Point2> xy_points(const py::object& a, const char* what) {
+    const auto xy = py::array_t<double, py::array::c_style | py::array::forcecast>::ensure(a);
+    if (!xy || xy.ndim() != 2 || xy.shape(1) != 2)
+        throw py::value_error(std::format("reduce_ring: {} must have shape (N, 2)", what));
+    std::vector<Point2> out(static_cast<std::size_t>(xy.shape(0)));
+    for (std::size_t i = 0; i < out.size(); ++i)
+        out[i] = Point2{xy.at(static_cast<py::ssize_t>(i), 0), xy.at(static_cast<py::ssize_t>(i), 1)};
+    return out;
+}
+
 }  // namespace
 
 PYBIND11_MODULE(_core, m) {
@@ -1004,5 +1016,52 @@ The catchment of a seed mask: every node draining into a seed (Priority-Flood).
 
 seed is (rows, cols), non-zero (or True) for a seed; any other shape is a
 ValueError. NoData is never in. Releases the GIL.
+)doc");
+
+    using terrain::vector_simplify::ReduceOutcome;
+    using terrain::vector_simplify::ReduceStatus;
+    py::enum_<ReduceStatus>(m, "ReduceStatus", R"doc(
+Why reduce_ring reduced a ring or did not: Ok, or a refusal of the input.
+)doc")
+        .value("Ok", ReduceStatus::Ok)
+        .value("InvalidTolerance", ReduceStatus::InvalidTolerance)
+        .value("NotCounterClockwise", ReduceStatus::NotCounterClockwise)
+        .value("TooFewVertices", ReduceStatus::TooFewVertices);
+
+    py::class_<ReduceOutcome>(m, "ReduceOutcome", R"doc(
+What reduce_ring returned: the reduced ring, open, a status, and the counts.
+)doc")
+        .def_property_readonly(
+            "ring",
+            [](const py::object& self) {
+                return point_view(self, self.cast<const ReduceOutcome&>().ring);
+            },
+            "Read-only (M, 2) float64, open: the first vertex is not repeated.")
+        .def_readonly("status", &ReduceOutcome::status)
+        .def_property_readonly("collinear", [](const ReduceOutcome& o) { return o.counts.collinear; })
+        .def_property_readonly("collapses", [](const ReduceOutcome& o) { return o.counts.collapses; })
+        .def_property_readonly("rejected_crossing",
+                               [](const ReduceOutcome& o) { return o.counts.rejected_crossing; })
+        .def_property_readonly("rejected_seed",
+                               [](const ReduceOutcome& o) { return o.counts.rejected_seed; })
+        .def_property_readonly("rejected_tolerance",
+                               [](const ReduceOutcome& o) { return o.counts.rejected_tolerance; });
+
+    m.def(
+        "reduce_ring",
+        [](const py::object& ring, double tolerance, const py::object& keep) {
+            const std::vector<Point2> points = xy_points(ring, "ring");
+            const std::vector<Point2> kept = xy_points(keep, "keep");
+            const py::gil_scoped_release unlocked;
+            return terrain::vector_simplify::reduce_ring<terrain::pred::DefaultKernel>(
+                points, tolerance, kept);
+        },
+        py::arg("ring"), py::arg("tolerance"), py::arg("keep"), R"doc(
+Reduce an open counter-clockwise ring to a tolerance, keeping its area and the
+keep-points inside (area-preserving segment collapse, exact crossing tests).
+
+ring is (N, 2) and keep (K, 2), float64-convertible; any other shape is a
+ValueError. Tolerance 0 only drops exactly collinear vertices. A refused input
+comes back as a status. Releases the GIL.
 )doc");
 }

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

@@ -0,0 +1,68 @@
+# Increment 22 acceptance: Bygdin (@perf, 2026-09-29): summary
+
+**Verdict: ACCEPTED** against the design's acceptance, which asks for the
+area to be within 2 % of NVE. The reduced catchment is **304.9096 km²
+against NVE's 305.54 km², -0.206 %**.
+
+There is **no performance baseline**: this is the first `catchment` run and
+the first mesh of this domain. Its speed figures are the baseline for the
+next run.
+
+- **Run conditions**: AC power at the start and end of every block. The
+  commit is `6757925`, with a Release `_core`. Method, tables and raw logs:
+  `bygdin/README.md`.
+
+**`rasputin catchment`** (seed 8.5425 61.3512, CORINE lake, DTM10):
+
+- **Windows**: three (4.35 M, 8.15 M and 16.94 M nodes). The design
+  expected two.
+- **Result**: 3,049,095 nodes. The fine outline has 17,812 vertices and the
+  reduced one 740, at 20 m. Both are 304.909550 km²; they differ by
+  -3e-07 m².
+- **Hausdorff distance, fine to reduced**: 19.92 m at 20 m.
+- **Times, median of 3**: flood 4.68 s, trace 0.13 s, reduce 0.04 s, wall
+  6.53 s.
+- **Memory**: max RSS 989 MB (peak footprint 920 MB).
+
+`--outline-tolerance` ladder (the area is equal at every step; all four
+outlines are valid):
+
+| tolerance | vertices | Hausdorff |
+|---|---|---|
+| 0 | 5,513 | 0 m |
+| 10 | 1,382 | 9.99 m |
+| 20 | 740 | 19.92 m |
+| 50 | 295 | 49.85 m |
+
+**Against NVE's polygon** (delfelt 1187, fetched without GDAL and not
+committed; the README gives the command and the sha256), comparing DEM nodes:
+
+- 99.12 % of NVE's nodes are in ours;
+- 99.33 % of ours are in NVE's;
+- 2.05 km² are ours only, and 2.67 km² are NVE's only.
+
+This is the design's prototype figure, reproduced by the built code.
+
+**`rasputin mesh --domain`** (10 threads, median of 3):
+
+| domain | tol | triangles | refine | total | max RSS | worst angle | max degree |
+|---|---|---|---|---|---|---|---|
+| reduced | 10 m | 52,198 | 0.034 s | 0.36 s | 646 MB | 0.356° | 12 |
+| fine | 10 m | 72,292 | 0.043 s | 0.39 s | 647 MB | 0.655° | 13 |
+| reduced | 1 m | 1,129,026 | 0.497 s | 0.93 s | 925 MB | 0.046° | 18 |
+| fine | 1 m | 1,137,703 | 0.483 s | 0.92 s | 942 MB | 0.083° | 16 |
+| reduced + CORINE | 10 m | 83,171 | 0.046 s | 1.97 s | 776 MB | 0.0041° | 16 |
+
+- **Checks**: every mesh is within its tolerance, has 0 uncovered nodes and
+  0 constrained-Delaunay violations.
+- **The reduction** saves 27.8 % of the triangles at 10 m and 0.8 % at 1 m.
+- **Land cover inside the catchment** (CORINE 2018):
+  - sparse vegetation (333) 42 %;
+  - bare rock (332) 21 %;
+  - heath (322) 17 %;
+  - water (512) 16 %;
+  - glacier (335) 2.4 %.
+
+**For the main session**:
+
+- No defects were found.