22 (1/2): auto-catchment, the fine catchment

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; 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 |
+| — | 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` |
 | — | `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

@@ -12,6 +12,7 @@
 #include <terrain/core/point.hpp>
 #include <terrain/core/pslg.hpp>
 #include <terrain/core/pslg_builder.hpp>
+#include <terrain/hydrology/upstream.hpp>
 #include <terrain/noding/node.hpp>
 #include <terrain/noding/noded_pslg_builder.hpp>
 #include <terrain/predicates/default_kernel.hpp>
@@ -291,6 +292,13 @@ template <typename T>
     return BoundRasterView{array, terrain::raster::RasterView<T>{geometry, a.data(), sentinel}};
 }
 
+// The flood's outcome with the shape its mask is viewed in.
+struct BoundUpstream {
+    terrain::hydrology::UpstreamOutcome outcome;
+    std::size_t rows;
+    std::size_t cols;
+};
+
 }  // namespace
 
 PYBIND11_MODULE(_core, m) {
@@ -945,5 +953,56 @@ constraint_feet inserts, for a worst node close to a constraint segment, its
 foot on the segment instead; off by default.
 A refused input comes back as a status; a mis-shaped array is a ValueError.
 Releases the GIL.
+)doc");
+
+    py::class_<BoundUpstream>(m, "UpstreamOutcome", R"doc(
+What upstream() returned: the catchment's node mask, its size and bounds, and
+whether it may continue past the window's edge or past NoData.
+)doc")
+        .def_property_readonly(
+            "mask",
+            [](const py::object& self) {
+                const auto& b = self.cast<const BoundUpstream&>();
+                const auto cols = static_cast<py::ssize_t>(b.cols);
+                return readonly_view<std::uint8_t>(self, b.outcome.mask.data(),
+                                                   {static_cast<py::ssize_t>(b.rows), cols},
+                                                   {cols, 1});
+            },
+            "Read-only (rows, cols) uint8: 1 for a node in the catchment, 0 otherwise.")
+        .def_property_readonly("nodes_in", [](const BoundUpstream& b) { return b.outcome.nodes_in; })
+        .def_property_readonly("row_min", [](const BoundUpstream& b) { return b.outcome.row_min; })
+        .def_property_readonly("row_max", [](const BoundUpstream& b) { return b.outcome.row_max; })
+        .def_property_readonly("col_min", [](const BoundUpstream& b) { return b.outcome.col_min; })
+        .def_property_readonly("col_max", [](const BoundUpstream& b) { return b.outcome.col_max; })
+        .def_property_readonly("touches_edge",
+                               [](const BoundUpstream& b) { return b.outcome.touches_edge; })
+        .def_property_readonly("touches_nodata",
+                               [](const BoundUpstream& b) { return b.outcome.touches_nodata; });
+
+    m.def(
+        "upstream",
+        [](const BoundRasterView& raster, const py::object& seed) {
+            using U8 = py::array_t<std::uint8_t, py::array::c_style | py::array::forcecast>;
+            const auto s = U8::ensure(seed);
+            const auto g = std::visit([](const auto& v) -> const auto& { return v.geometry(); },
+                                      raster.view);
+            if (!s || s.ndim() != 2 || static_cast<std::size_t>(s.shape(0)) != g.rows()
+                || static_cast<std::size_t>(s.shape(1)) != g.cols())
+                throw py::value_error(std::format(
+                    "upstream: the seed mask's shape must be the raster's ({}, {})", g.rows(),
+                    g.cols()));
+            const std::span<const std::uint8_t> seeds{s.data(), g.size()};
+            // Every buffer read below is held by a local or by `raster`.
+            const py::gil_scoped_release unlocked;
+            return BoundUpstream{
+                std::visit([&](const auto& v) { return terrain::hydrology::upstream(v, seeds); },
+                           raster.view),
+                g.rows(), g.cols()};
+        },
+        py::arg("view"), py::arg("seed"), R"doc(
+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");
 }