From 1c6ead27186a2e0ad70520e8c759509e523a992a Mon Sep 17 00:00:00 2001 From: Ola Skavhaug Date: Tue, 29 Sep 2026 00:38:19 +0200 Subject: [PATCH 01/14] ROADMAP: auto-catchment must reduce its outline (Ola) Co-Authored-By: Claude Opus 5.5 --- ROADMAP.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ROADMAP.md b/ROADMAP.md index 12348301..289a4ded 100644 --- a/ROADMAP.md +++ b/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 | planned: ruled by Ola 2026-09-28 (16b's Q2) as a separate increment, next after 16b | `docs/increments/16b-terrain-polygons.md` (Q2, R6) | | 16d | A record of the input polylines in the mesh. Today a constraint edge keeps only its type bits; the source feature, its attributes (a river's name or order, a road's class) and which edges belong to one polyline are lost, and where edges from two features merge only the union of their bits survives. Proposed: a feature table (one row per input feature: source ID and attributes) and, per constraint edge, the list of features it came from (a list, since a merged edge belongs to several). The noder already tracks each piece's input chain; the open part is the output: a per-edge list in `.vtk`/`.ply`, or a side file | proposed by Ola 2026-09-28; to be designed by `@architect` once Ola has said what it is for; after 16b | - | | 20c | Soft quality criterion: a penalty that each Steiner node or constraint split must pay for in angle gained, instead of 20's hard 25°; applied at the start and during DEM refinement; may split constraint segments when that improves the mesh. Ola's rulings on 20's C1-C3 | to design after 16b (`@architect` measures cost against 20 first) | `docs/increments/20-start-quality.md` (Ola's rulings) | -| — | Auto-catchment: the watershed upstream of a coordinate, computed from the DEM and handed to `--domain`, so a catchment no longer has to be supplied as a file (Ola, 2026-09-27: "not far into the future"). The textbook route is depression handling (Priority-Flood, Barnes, Lehman and Mulla 2014), D8 flow directions (O'Callaghan and Mark 1984) and accumulation, the pour point snapped to the strongest flow nearby, the upstream cells traced and their outline turned into a polygon; the literature check is `@architect`'s. Open for its design: whether it runs in the C++ core (a 10 m tile is 25 M cells); how a stair-stepped cell outline becomes a domain polygon, which meets input coarsening; and that a real catchment crosses tile edges, so it needs gap 6 (a DEM in several tiles) first. Legacy has nothing on it (`grep -rliE "watershed|flow.?acc|flow.?dir|pour.?point|catchment" legacy` returns no files) | to design; after gap 6, which it needs; placed after 20c, can move ahead of it on Ola's word | none yet | +| — | 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); and that a real catchment crosses tile edges, so it needs gap 6 (a DEM in several tiles) first. Legacy has nothing on it (`grep -rliE "watershed|flow.?acc|flow.?dir|pour.?point|catchment" legacy` returns no files) | to design; after gap 6, which it needs; placed after 20c, can move ahead of it on Ola's word | none yet | | — | `raster/`: grid-to-world geometry and bilinear sampling | shipped (`7785fea`), **no record** | none — predates the protocol | ## What stands between here and an operational MVP From f8236735c89e756532c0bfbc1b24efc7eaf584a1 Mon Sep 17 00:00:00 2001 From: Ola Skavhaug Date: Tue, 29 Sep 2026 00:42:39 +0200 Subject: [PATCH 02/14] ROADMAP: auto-catchment outline through DEM nodes, area kept (Ola) Co-Authored-By: Claude Opus 5.5 --- ROADMAP.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ROADMAP.md b/ROADMAP.md index 289a4ded..9118016c 100644 --- a/ROADMAP.md +++ b/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 | planned: ruled by Ola 2026-09-28 (16b's Q2) as a separate increment, next after 16b | `docs/increments/16b-terrain-polygons.md` (Q2, R6) | | 16d | A record of the input polylines in the mesh. Today a constraint edge keeps only its type bits; the source feature, its attributes (a river's name or order, a road's class) and which edges belong to one polyline are lost, and where edges from two features merge only the union of their bits survives. Proposed: a feature table (one row per input feature: source ID and attributes) and, per constraint edge, the list of features it came from (a list, since a merged edge belongs to several). The noder already tracks each piece's input chain; the open part is the output: a per-edge list in `.vtk`/`.ply`, or a side file | proposed by Ola 2026-09-28; to be designed by `@architect` once Ola has said what it is for; after 16b | - | | 20c | Soft quality criterion: a penalty that each Steiner node or constraint split must pay for in angle gained, instead of 20's hard 25°; applied at the start and during DEM refinement; may split constraint segments when that improves the mesh. Ola's rulings on 20's C1-C3 | to design after 16b (`@architect` measures cost against 20 first) | `docs/increments/20-start-quality.md` (Ola's rulings) | -| — | Auto-catchment: the watershed upstream of a coordinate, computed from the DEM and handed to `--domain`, so a catchment no longer has to be supplied as a file (Ola, 2026-09-27: "not far into the future"). The textbook route is depression handling (Priority-Flood, Barnes, Lehman and Mulla 2014), D8 flow directions (O'Callaghan and Mark 1984) and accumulation, the pour point snapped to the strongest flow nearby, the upstream cells traced and their outline turned into a polygon; the literature check is `@architect`'s. Open for its design: whether it runs in the C++ core (a 10 m tile is 25 M cells); how a stair-stepped cell outline becomes a domain polygon, which meets input coarsening; **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); 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 outline runs through the catchment's outermost DEM nodes (8-connected), and the reduction keeps a subset of them, so every vertex stays a DEM node (Ola, 2026-09-29: "Isn't the DEM points, really?"); 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) | to design; after gap 6, which it needs; placed after 20c, can move ahead of it on Ola's word | none yet | | — | `raster/`: grid-to-world geometry and bilinear sampling | shipped (`7785fea`), **no record** | none — predates the protocol | ## What stands between here and an operational MVP From 5504d60dd979499d9aee90954cca87fcb0dfab02 Mon Sep 17 00:00:00 2001 From: Ola Skavhaug Date: Tue, 29 Sep 2026 00:44:29 +0200 Subject: [PATCH 03/14] ROADMAP: auto-catchment's reduced outline may use off-node vertices (Ola's correction) Co-Authored-By: Claude Opus 5.5 --- ROADMAP.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ROADMAP.md b/ROADMAP.md index 9118016c..cd49e64c 100644 --- a/ROADMAP.md +++ b/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 | planned: ruled by Ola 2026-09-28 (16b's Q2) as a separate increment, next after 16b | `docs/increments/16b-terrain-polygons.md` (Q2, R6) | | 16d | A record of the input polylines in the mesh. Today a constraint edge keeps only its type bits; the source feature, its attributes (a river's name or order, a road's class) and which edges belong to one polyline are lost, and where edges from two features merge only the union of their bits survives. Proposed: a feature table (one row per input feature: source ID and attributes) and, per constraint edge, the list of features it came from (a list, since a merged edge belongs to several). The noder already tracks each piece's input chain; the open part is the output: a per-edge list in `.vtk`/`.ply`, or a side file | proposed by Ola 2026-09-28; to be designed by `@architect` once Ola has said what it is for; after 16b | - | | 20c | Soft quality criterion: a penalty that each Steiner node or constraint split must pay for in angle gained, instead of 20's hard 25°; applied at the start and during DEM refinement; may split constraint segments when that improves the mesh. Ola's rulings on 20's C1-C3 | to design after 16b (`@architect` measures cost against 20 first) | `docs/increments/20-start-quality.md` (Ola's rulings) | -| — | Auto-catchment: the watershed upstream of a coordinate, computed from the DEM and handed to `--domain`, so a catchment no longer has to be supplied as a file (Ola, 2026-09-27: "not far into the future"). The textbook route is depression handling (Priority-Flood, Barnes, Lehman and Mulla 2014), D8 flow directions (O'Callaghan and Mark 1984) and accumulation, the pour point snapped to the strongest flow nearby, the upstream cells traced and their outline turned into a polygon; the literature check is `@architect`'s. Open for its design: whether it runs in the C++ core (a 10 m tile is 25 M cells); how a stair-stepped cell outline becomes a domain polygon, which meets input coarsening; **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 outline runs through the catchment's outermost DEM nodes (8-connected), and the reduction keeps a subset of them, so every vertex stays a DEM node (Ola, 2026-09-29: "Isn't the DEM points, really?"); 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) | 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) | to design; after gap 6, which it needs; placed after 20c, can move ahead of it on Ola's word | none yet | | — | `raster/`: grid-to-world geometry and bilinear sampling | shipped (`7785fea`), **no record** | none — predates the protocol | ## What stands between here and an operational MVP From 684f1f52634ea595b4a19dd0572453cb5f737f53 Mon Sep 17 00:00:00 2001 From: Ola Skavhaug Date: Tue, 29 Sep 2026 01:09:45 +0200 Subject: [PATCH 04/14] 22 design: auto-catchment (Bygdin first) Co-Authored-By: Claude Opus 5.5 --- ROADMAP.md | 2 +- docs/increments/22-auto-catchment.md | 651 +++++++++++++++++++++++++++ 2 files changed, 652 insertions(+), 1 deletion(-) create mode 100644 docs/increments/22-auto-catchment.md diff --git a/ROADMAP.md b/ROADMAP.md index cd49e64c..c70a0225 100644 --- a/ROADMAP.md +++ b/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 | planned: ruled by Ola 2026-09-28 (16b's Q2) as a separate increment, next after 16b | `docs/increments/16b-terrain-polygons.md` (Q2, R6) | | 16d | A record of the input polylines in the mesh. Today a constraint edge keeps only its type bits; the source feature, its attributes (a river's name or order, a road's class) and which edges belong to one polyline are lost, and where edges from two features merge only the union of their bits survives. Proposed: a feature table (one row per input feature: source ID and attributes) and, per constraint edge, the list of features it came from (a list, since a merged edge belongs to several). The noder already tracks each piece's input chain; the open part is the output: a per-edge list in `.vtk`/`.ply`, or a side file | proposed by Ola 2026-09-28; to be designed by `@architect` once Ola has said what it is for; after 16b | - | | 20c | Soft quality criterion: a penalty that each Steiner node or constraint split must pay for in angle gained, instead of 20's hard 25°; applied at the start and during DEM refinement; may split constraint segments when that improves the mesh. Ola's rulings on 20's C1-C3 | to design after 16b (`@architect` measures cost against 20 first) | `docs/increments/20-start-quality.md` (Ola's rulings) | -| — | Auto-catchment: the watershed upstream of a coordinate, computed from the DEM and handed to `--domain`, so a catchment no longer has to be supplied as a file (Ola, 2026-09-27: "not far into the future"). The textbook route is depression handling (Priority-Flood, Barnes, Lehman and Mulla 2014), D8 flow directions (O'Callaghan and Mark 1984) and accumulation, the pour point snapped to the strongest flow nearby, the upstream cells traced and their outline turned into a polygon; the literature check is `@architect`'s. Open for its design: whether it runs in the C++ core (a 10 m tile is 25 M cells); how a stair-stepped cell outline becomes a domain polygon, which meets input coarsening; **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) | 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) | **designed as increment 22** (`@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 diff --git a/docs/increments/22-auto-catchment.md b/docs/increments/22-auto-catchment.md new file mode 100644 index 00000000..5228c918 --- /dev/null +++ b/docs/increments/22-auto-catchment.md @@ -0,0 +1,651 @@ +# Increment 22 — auto-catchment: the catchment of a lake, from the DEM (Bygdin first) + +Status: **designed** (`@architect`, 2026-09-29), on branch +`increment22-autocatchment` off master `b4847d7`. Nothing is built. Ola was +asleep while this was written; every choice he would normally make is marked +"Default (main session / @architect, 2026-09-29), for Ola to confirm", with +the alternative, so the loop can run tonight. + +**Closes.** The auto-catchment row of `ROADMAP.md`: a catchment computed from +the DEM, reduced to a polygon a mesh can afford, and handed to `--domain`. +After this increment: + +```sh +rasputin catchment --dem ../rasputin_data/DTM10_UTM33_20260925 \ + --seed 8.5425 61.3512 \ + --lakes ../rasputin_data/corine2018_dtm10_utm33.gpkg --lakes-layer corine2018 \ + --out bygdin.geojson +rasputin mesh --dem ../rasputin_data/DTM10_UTM33_20260925 \ + --domain bygdin.geojson --tolerance 1 --out bygdin.vtk +``` + +**Not closed.** Stream extraction and Strahler order (`auto_catchments.md`). +Flow accumulation, and with it snapping a pour point to the strongest flow +nearby. A catchment whose area runs off the data (a border with Sweden, the +sea): it is refused, not written. Holes in the catchment are filled, not kept. +Parallel flooding (Barnes 2016). Moving the outline tracer into C++ (below, +"Where it runs"). The depression case under "The window" is a known limit. + +## Ola's requirements, quoted + +From the ROADMAP row, 2026-09-29: + +1. "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." The traced outline is reduced before meshing, to a stated + horizontal tolerance. +2. "the resulting polygon has approximately the same area as the fine one". + The reduction keeps the area. Off-node vertices ("We have aleady + established that the polygons and DEM don't need to match") make exact + area possible. +3. The reduced polygon is simple (no self-crossings) and contains the seed. + +And on what the DEM is: "Isn't the DEM points, really?" The catchment below is +a set of DEM nodes, and its outline is drawn between nodes, never around cell +squares. + +Tonight's goal, verbatim: "I would really like you to work while I sleep +tonight, to get a first working example of the autocatchment. If possible I +would like the Norwegian lake Bygdin to be the seed for the calculation." + +## What the data says about Bygdin (measured 2026-09-29) + +These were measured before the design, because two of the facts in the brief +did not hold, and the design depends on the third. + +- **The coordinates in the brief are off.** UTM33 (177934, 6818339), given as + the dam, is on a slope at 1580 m north of Vinstre, 10 km east of Bygdin and + outside its catchment. NVE's regulation dam, "BYGDIN REGULERINGSDAM" (dam + 1192), is at (168087, 6815337) in EPSG:25833, about 8.792 E, 61.330 N. The + lake runs from x 142.3 km to 168.5 km, so its west end is in tile + `6801_3`, not `6801_2`: the catchment crosses a tile edge, as 15a intended. +- **The lake is not one flat in DTM10.** Of the 399,371 DEM nodes inside + CORINE's Bygdin polygon, 101,058 are exactly 1057.4 m; the rest scatter + between 1055.4 and 1059.7 m (1st and 99th percentile). The largest + connected exactly-flat region is 7.1 km² of a 40 km² lake. A seed defined as + "the flat under the point" would miss most of the lake. +- **CORINE has the lake.** Feature `fid` 54101 of `corine2018` in + `corine2018_dtm10_utm33.gpkg`, class 512, 39.94 km², one part, no holes. + NVE gives 40.03 km² at the highest regulated level. +- **NVE publishes the catchment.** NVE's reservoir catchment ("delfelt") 1187, + "BYGDIN", has `delfeltAreal_km2` = **305.54**, and its upstream list is + itself alone, so it is the whole catchment. Source: NVE's map service, + `https://gis3.nve.no/map/rest/services/Mapservices/VassdragsreguleringVannkraft/MapServer/8/query?where=delfeltNr%3D1187&outFields=*&f=json`, + fetched 2026-09-29; the same query with `returnGeometry=true&outSR=25833&f=geojson` + returns the polygon (1376 vertices, 305.539 km² by shapely). The reservoir + record is layer 9 (`magasinNavn='BYGDIN'`, 40.03 km²). Norwegian Wikipedia's + Bygdin article gives 305.59 km² citing NVE's REGINE units (2015); the NVE + service is the number used here. +- **The method below was prototyped** (throwaway Python in the scratchpad, not + committed) over a 4201 x 3101-node window of 10 m DTM10: 3,049,095 nodes, + 304.91 km², against NVE's 305.54 km² (-0.21 %); 99.1 % of NVE's nodes are + in ours and 99.3 % of ours in NVE's. The traced outline had 17,804 + vertices. shapely's topology-preserving Douglas-Peucker, used only as a + rough gauge, took it to 924 vertices at 20 m. The acceptance run measures + these numbers again with the real code. + +## Prior art: legacy and literature + +**Legacy.** Nothing to carry over. The ROADMAP's grep still returns no files: + +``` +$ grep -rliE "watershed|flow.?acc|flow.?dir|pour.?point|catchment" legacy +$ echo $? +1 +``` + +For the reduction: + +``` +$ grep -rliE "simplif|visvalingam|douglas|peucker" legacy +legacy/rasputin/triangulate_dem.h +legacy/rasputin/geo_tiff_reader.py +legacy/rasputin/mesh.py +legacy/rasputin/application.py +legacy/tests/test_gml_repository.py +legacy/bindings.cpp +``` + +Every hit is CGAL's surface-mesh edge collapse (Lindstrom-Turk cost), which +simplifies a triangle mesh, not a polygon. Nothing is carried over. + +**Literature.** Each reference below was checked against Crossref on +2026-09-29 (title, authors, venue, DOI); abstracts were read where Crossref or +Semantic Scholar had them. No paper was read in full tonight, and where the +design leans on a detail beyond the abstract, it says so. + +- **Barnes, Lehman and Mulla 2014a**, "Priority-Flood: An optimal + depression-filling and watershed-labeling algorithm for digital elevation + models", *Computers & Geosciences* 62:117-127, doi:10.1016/j.cageo.2013.04.024 + (open copy arXiv:1511.04463). Floods the DEM inward from its edges with a + priority queue; the abstract says it "can also be adapted to label + watersheds and determine flow directions". **This increment builds on it.** + What differs: the flood labels one catchment (the lake's) rather than every + edge outlet's, and ties are broken first in, first out by a push counter so + the result is deterministic. Barnes's faster variant (a plain queue inside + depressions) is an optimisation that can come later. +- **Metz, Mitasova and Harmon 2011**, "Efficient extraction of drainage + networks from massive, radar-based elevation models with least cost path + search", *HESS* 15:667-678, doi:10.5194/hess-15-667-2011. GRASS's + `r.watershed`: routing by least-cost search instead of filling first. The + flood here routes the same way in effect: each node drains to the node that + flooded it. +- **O'Callaghan and Mark 1984**, "The extraction of drainage networks from + digital elevation data", *Computer Vision, Graphics, and Image Processing*, + doi:10.1016/S0734-189X(84)80047-X. Crossref lists it as volume 27(2), page + 247; it is commonly cited as 28(3):323-344. Cite it by DOI until someone + resolves which. D8: each node drains to its steepest neighbour, slope + measured with the diagonal's length. **Departure, and why:** the flood + drains each node to its *lowest* filled neighbour, not its steepest, so the + diagonal distance plays no part. The two differ only where a diagonal + neighbour is lower but less steep, which moves a divide by a node here and + there. In exchange there is one pass instead of three (fill, flats, + directions), and no flat resolution at all (next point). D8 comes when + stream extraction needs accumulation. +- **Garbrecht and Martz 1997**, "The assignment of drainage direction over + flat surfaces in raster digital elevation models", *J. Hydrology* + 193:204-213, doi:10.1016/S0022-1694(96)03138-1; and **Barnes, Lehman and + Mulla 2014b**, "An efficient assignment of drainage direction over flat + surfaces in raster digital elevation models", *Computers & Geosciences* + 62:128-135, doi:10.1016/j.cageo.2013.01.009. Both give flats realistic flow + paths (towards lower terrain, away from higher). **Not needed here, and + why:** which catchment a flat belongs to depends only on where it spills, + and the flood gets that exactly; the path water takes across the flat does + not change the answer. The lake itself is the seed, so its surface, flat or + noisy, is never routed at all. The one case the flow path decides is a flat + that spills over two outlets at the same level: the flood splits it by + distance in steps from each outlet. That is recorded as a limit, not + solved. +- **Wang and Liu 2006** (*IJGIS* 20:193-213, doi:10.1080/13658810500433453) + and **Lindsay 2016** (*Hydrological Processes* 30:846-857, + doi:10.1002/hyp.10648): filling and breaching. Priority-Flood fills; breach + is out of scope, as `auto_catchments.md` has it. +- **Kong and Rosenfeld 1989**, "Digital topology: Introduction and survey", + *CVGIP* 48:357-393, doi:10.1016/0734-189X(89)90147-3. The (8, 4) pairing: + the catchment is 8-connected (water steps diagonally), so the outside must + be taken 4-connected for the outline to be a simple curve. The tracer's + saddle rule below is that pairing. **Lorensen and Cline 1987**, "Marching + cubes", *SIGGRAPH Computer Graphics* 21(4):163-169, doi:10.1145/37402.37422, + is the family the tracer belongs to (its 2-D case, marching squares). +- **Kronenfeld, Stanislawski, Buttenfield and Brockmeyer 2020**, + "Simplification of polylines by segment collapse: minimizing areal + displacement while preserving area", *International Journal of Cartography* + 6(1):22-46, doi:10.1080/23729333.2019.1631535 (online 2019). APSC. From its + abstract: segments are collapsed to Steiner points in priority order, with + placement and displacement functions chosen so that area is preserved + exactly, and "self-intersections can be avoided by testing for + intersections with two new line segments associated with each segment + collapse operation". **The reduction builds on it.** What differs: (a) the + stopping rule is a horizontal tolerance against the fine outline, because + that is what Ola asked for, where APSC targets a vertex count or scale; + (b) the priority is that same deviation, not areal displacement; (c) the + placement below is ours. The paper is paywalled and its placement function + was not read. The area guarantee does not depend on it: any point on the + area-preserving line keeps the area. (d) A containment check for the seed + is added. +- **Buchin, Meulemans, van Renssen and Speckmann 2016**, "Area-preserving + simplification and schematization of polygonal subdivisions", *ACM TSAS* + 2(1), doi:10.1145/2818373. The edge-move, which also keeps area and + topology. The alternative if APSC's placement proves poor on lattice + outlines. +- **Bose, Cabello, Cheong, Gudmundsson, van Kreveld and Speckmann 2006**, + "Area-preserving approximations of polygonal paths", *J. Discrete + Algorithms* 4(4):554-566, doi:10.1016/j.jda.2005.06.008. The optimisation + version (fewest vertices). This increment is greedy and claims no minimum. +- **Visvalingam and Whyatt 1993** (*Cartographic J.* 30:46-51, + doi:10.1179/000870493786962263) and **Saalfeld 1999** (*CaGIS* 26:7-18, + doi:10.1559/152304099782424901): the vertex-removal and topology-consistent + simplifiers. Neither keeps area; they are why off-node vertices matter. + +**Novelty.** None is claimed. Seeding a watershed with a polygon of pour +points is standard in GIS tools; lake-as-seed with Priority-Flood labelling, +a marching-squares outline and APSC with a tolerance band is a combination of +published parts. If a later write-up wants to claim something (for example +the tolerance-band guarantee against the fine outline), the check is still to +be done. + +## The design + +### Data flow + +``` +cli.py catchment + | parses flags into a CatchmentRequest (frozen Pydantic), builds the + | repository (paths stop here), calls delineate(), writes GeoJSON + v +catchment.py delineate(request, repository) -> Catchment [no paths] + 1. seed point -> DEM CRS (crs.reprojector); lake polygon containing it + -> DEM CRS (seed.py, shapely/pyproj) + 2. window = seed bounds grown by the margin + 3. loop: plan_mosaic + assemble (15a) -> DemTile, read-only + seed mask: DEM nodes inside the lake (shapely.contains_xy) + _core.upstream(raster_view, seed_mask) -> mask, bbox, flags [C++] + contained? done : grow the window and repeat + 4. outline.trace(mask) -> rings, lattice units (numpy, no _core) + pick the ring around the seed; drop the rest, count them + 5. _core.reduce_ring(ring, tolerance, keep=[seed]) [C++] + 6. -> Catchment(fine, reduced, counts, areas, windows) frozen +``` + +The C++ core sees an elevation array, a byte mask, and a ring of doubles in +metres. No path, file or CRS crosses (CLAUDE.md §2, I/O boundary). + +### Where it runs + +- **The flood is C++.** It visits every node of the window: 12-13 M nodes for + Bygdin, far more for Glomma. The Python prototype took 54 s on its + 13 M-node Bygdin window at 10 m; a priority queue in C++ should take a + few seconds. + `include/terrain/hydrology/upstream.hpp`, header-only, a template on the + `RasterSource` concept, bound over the existing `RasterView` variant, GIL + released. +- **The reduction is C++**, `include/terrain/vector_simplify/area_collapse.hpp`, + because its crossing tests need the exact kernel (`noding/intersect.hpp`'s + `classify`, `DefaultKernel` in the binding). Both module names are the + ones `project_structure.md` already plans. +- **The tracer is Python and numpy**, `src_python/tin_engine/outline.py`, + never importing `_core`. It is raster topology on a byte mask, O(boundary) + after one vectorised pass, and the prototype traced Bygdin in 0.09 s. Doing + it in C++ would add a binding and a C++ suite tonight for no measured gain. + Default (main session / @architect, 2026-09-29), for Ola to confirm; + alternative: C++ next to the flood, when a Glomma-size run says so. +- **Planning and I/O are Python**, as 15a's `mosaic.py` and `dem_input.py` + are. `catchment.py` takes a `DemRepository` (the Protocol in + `io/repository.py`) rather than paths, so tests pass an in-memory one; the + five lines in `open_dem` that build a repository from sources move to a + helper in `dem_input.py` that both commands call. + +### The seed + +**Default (main session / @architect, 2026-09-29), for Ola to confirm: the +seed is a lake polygon, and the catchment is every DEM node that drains into +any node inside it.** The user gives a point (`--seed X Y`, in `--seed-crs`, +default EPSG:4326 so lon lat), and a polygon source (`--lakes`); the polygon +containing the point is the lake. Its DEM nodes are all seeds. For Bygdin: +`--seed 8.5425 61.3512` (mid-lake, UTM33 about (155000, 6819000), on the +1057.4 m surface, inside CORINE feature 54101). + +Why this and not the alternatives: + +- *A pour point at the dam, snapped to the strongest flow nearby*: needs flow + accumulation (a second pass and a stored order), a snap radius that can + jump to the wrong channel, and flat routing across the lake to reach the + outlet. More code and more ways to be wrong. The brief's dam coordinate + shows how easily the point itself is off. +- *The DEM flat under the point*: measured above, DTM10's Bygdin is not one + flat; the flat would be 7 km² of 40. +- *The lake polygon*: no snapping, no flats, no dependence on where exactly + the outlet is. Its cost: the polygon comes from outside the DEM. If it + reaches past the outlet, whatever drains into that stretch of river joins + the catchment. CORINE's Bygdin polygon ends at x 168460, NVE's catchment at + x 168680 and the dam at x 168087, so a sliver of river below the dam may be + included. The acceptance run measures it (nodes in ours and not in NVE's). + +The source is any polygon layer: a GeoPackage table (`--lakes-layer`, read +with 16b's `io/geopackage.py`, box query on the seed point, in the layer's +own CRS) or a GeoJSON `FeatureCollection`. A multipolygon contributes the +part containing the point. No class filter: the polygon under the point is +the one meant. Refusals: the point in no polygon, or in two (overlapping +input); both name the point in the source's CRS. + +**Without `--lakes`**, the seed is the one DEM node nearest the point, a pour +point with no snapping. It is there because it costs five lines and gives the +synthetic tests a seed without a polygon file. stderr says that a pour point +must lie on the flow line. Snapping is a later increment. + +### Flow and membership: one flood + +`upstream(z, seed)`: + +- Every valid node on the window's edge, and every valid node with a NoData + 8-neighbour, is an outlet: pushed with its own z at the start. NoData nodes + (the sentinel or NaN, as `RasterView::is_nodata` says) are never pushed and + never in the catchment. +- Pop the lowest key (level, push counter). For each 8-neighbour not yet + reached: its level is max(its z, the popped level); its label is *in* if it + is a seed, else the popped node's label; push it. +- Outlets are labelled *in* if they are seeds, *out* otherwise. + +The catchment is the nodes labelled *in*. A node is in it exactly when the +chain of nodes that flooded it, which is its drainage path over the filled +surface, passes through the lake. Levels are `double` (a float32 z is exact +in it). The counter makes equal levels first in, first out, so the result +depends on nothing but the input. + +Memory: the elevation array is 15a's, borrowed. The flood's own is one byte +per node (unreached, out, in), which is also the returned mask, plus the queue +(24 bytes per entry; the frontier, not the window, in practice). The seed mask +is one more byte per node, owned by numpy. + +Returned: the mask (a numpy `uint8` array the outcome owns), the number of +nodes in, the bounding rows and columns of the in-nodes, and two flags: +`touches_edge` (an in-node on the window's edge) and `touches_nodata` (an +in-node with a NoData neighbour). A seed mask whose shape is not the raster's +is a `ValueError` in the binding. + +### The window + +A read region grown until the catchment clearly lies inside it, the same +re-plan pattern as 15b's `_domain_plan`: + +1. Start from the seed's bounds (the lake polygon's, or the point) grown by + the margin, `WINDOW_MARGIN_M = 2000` metres. Default (main session / + @architect, 2026-09-29), for Ola to confirm; no flag tonight. +2. Plan and assemble that box (15a), flood it. +3. If the in-nodes' bounds grown by the margin fit inside the window, stop. + Otherwise the new box is the old box joined with the in-nodes' bounds grown + by twice the margin, the margin doubles, and the loop repeats from 2. The + box only grows and the margin doubles, so the total work is within about + twice the last flood's. +4. Refuse, with the side named, when the catchment reaches a side the data + cannot extend (the planned window is smaller than the box asked for on + that side), or touches NoData: the catchment is truncated, and a truncated + catchment is a wrong one. Default (main session / @architect, + 2026-09-29), for Ola to confirm; alternative: write it with a warning + under an `--allow-truncated` flag. +5. Memory: before each flood, refuse if nodes x (itemsize + 2) exceeds half + the physical memory (`mosaic.physical_memory`, dtype-aware as 15a's cap + is), naming the window's size. 15a's own cap on the array still applies + inside `plan_mosaic`. + +**Known limit.** Not touching the window's edge does not prove the catchment +complete. A closed depression that straddles the window's edge drains out +through the edge in the window, but in the full DEM it may fill and spill into +the catchment. The margin makes this unlikely and each growth step makes it +less likely; nothing here rules it out. The acceptance run's comparison +against NVE is the check on Bygdin. + +For Bygdin the first box is the lake's bounds plus 2 km: x 140.3-170.5 km, +y 6811.2-6825.6 km. NVE's catchment reaches y 6832.3 km, so one growth step is +expected, to a window of roughly 4100 x 3000 nodes (12 M, 49 MB at float32). + +### The fine outline + +`outline.trace(mask) -> list of rings`, each a closed ring of `(row, col)` in +lattice units, as halves (every vertex is the midpoint of a lattice edge +between an in-node and an out-node). This is marching squares on the 0/1 +node values: in each square of four nodes, a boundary segment joins the +midpoints of its in/out edges, oriented with the in-nodes on the left. The +one ambiguous square, two in-nodes on one diagonal, is resolved by keeping the +two in-nodes connected (the (8, 4) pairing): the catchment is 8-connected +because water steps diagonally. The mask is padded by one row and column of +out-nodes so every ring closes. + +**The guarantee.** Every ring is a simple closed polygon, and no two rings +share a point. Each boundary lattice edge has exactly one midpoint, used by +exactly one incoming and one outgoing segment; two segments in one square +never meet (the saddle case cuts two opposite corners); segments of +different squares meet only at shared midpoints. Every in-node lies strictly +inside an odd number of rings and every out-node inside an even number, at a +distance of at least a quarter of the cell's diagonal. Outer rings are +counter-clockwise in the world frame (x east, y north), holes clockwise. + +**Holes and pinches.** A pinch (two parts meeting diagonally) is not a +special case: the saddle rule joins them, so the traced ring is simple. A +spur one node wide becomes a thin, simple ring around it. **Holes are filled**: +the fine outline is the one outer ring that contains the seed point (the +user's point with a lake, the pour node without), and every other ring is +dropped. stderr reports how many outer rings were dropped with their node +counts, and how many holes were filled with their area. Default (main session +/ @architect, 2026-09-29), for Ola to confirm; alternative: keep holes as +domain holes (16 supports them), which then need reducing too. If no ring +contains the seed point, `CatchmentError` (a data oddity: the point in a gap +the flood did not reach). + +**Area.** The fine outline's area is close to the in-node count times the +cell area: a straight side sits half a cell outside the in-nodes, and each +convex corner cuts off an eighth of a cell. On Bygdin (prototype) 304.909 km² +against 304.9095 km². The fine area is the reference the reduction keeps. + +Why not a polygon through the boundary nodes (the brief's "8-connected, +diagonal steps"): it lies half a cell inside the node area along the whole +perimeter (up to about 0.7 km² on Bygdin: half a cell times the fine outline's 140 km), and it +touches itself at every pinch, which needs a repair pass with its own proofs. +The midpoint ring needs none. Default (main session / @architect, +2026-09-29), for Ola to confirm. + +The tracer converts to world coordinates only at the end: `x = x_min + col * +delta_x`, `y = y_max - row * delta_y`; halves of a 10 m lattice on 5 m +multiples are exact in doubles. + +### The reduction + +`reduce_ring(ring, tolerance, keep) -> ReduceOutcome` in +`vector_simplify/area_collapse.hpp`. Input: a simple counter-clockwise ring of +`Point2` in metres (Python subtracts the window's lower-left corner first, so +coordinates stay below about 10^5 and areas lose nothing), a tolerance in +metres, and points that must stay inside (the seed point). Output: the reduced +ring, a status (`Ok`, `InvalidTolerance` for negative or non-finite, +`NotCounterClockwise`, `TooFewVertices` below 4), and counts (collinear +vertices dropped, collapses made, candidates rejected for crossing, for the +seed, for tolerance). + +1. **Collinear pass.** Drop every vertex exactly collinear with its two + neighbours and between them (exact `orient2d` is zero). The lattice outline + is full of them. This changes neither area nor shape. +2. **Candidates.** For each edge B-C with neighbours A before and D after, + the collapse replaces B and C by one new point E. E lies on the line + parallel to A-D at the signed distance that keeps the area of A-B-C-D + equal to that of A-E-D (APSC's area rule). On that line, E is where it + meets line A-B or line C-D, whichever gives the smaller deviation (ties to + A-B); if both are parallel to it, the foot of B-C's midpoint. +3. **Deviation.** Each current edge carries the range of fine-ring vertices it + stands for (from the original fine ring, before the collinear pass). For a + candidate, the range is from A-B's start to C-D's end. Its deviation is the + larger of: the greatest distance from a fine vertex in that range to the + chain A-E-D, and E's distance to the fine segments in the range. A + candidate is admissible when its deviation is at most the tolerance. The + reference is always the fine ring, so errors never accumulate. +4. **Order.** A heap on (deviation, id of B), where ids are the fine + indices and then a counter for new points. Least deviation first, so the + outline moves as little as possible for each vertex it loses. +5. **Checks before a collapse**, all with the exact kernel on the doubles as + given: + - *Simple*: A-E and E-D do not meet any current edge except at A (with the + edge ending at A) and at D (with the edge starting at D), and do not + overlap those two. A uniform grid of current edges (bucket side the + tolerance, at least one cell) finds the candidates; it is updated on + every collapse. + - *Seed inside*: each keep-point has winding number zero around the closed + loop A-B-C-D-E-A and lies on neither new edge. (The difference between + its winding in the old ring and in the new one is exactly that loop's.) + - A rejected candidate waits until one of its four vertices changes. +6. **Apply** the best admissible candidate, re-evaluate the candidates whose + four vertices include A, E or D, and repeat until none is admissible or 4 + vertices remain. Stale heap entries are skipped by a per-vertex version. + +With tolerance 0 only the collinear pass runs. + +**Guarantees**, each with the check the tests make: + +- *Area*: equal to the fine ring's up to rounding. Tested as + `|A_reduced - A_fine| <= 1e-9 * A_fine`; the acceptance run reports the + difference in m². +- *Simple*: no two non-adjacent edges meet, adjacent edges share only their + vertex. Tested by shapely `is_valid` and, on small cases, a brute-force + pairwise test with exact orientation. +- *Seed inside*: the seed point is strictly inside (shapely `contains`). +- *Tolerance*: every fine vertex is within the tolerance of the reduced ring, + and every reduced vertex within the tolerance of the fine ring. So the + symmetric Hausdorff distance between the two is at most the tolerance plus + one cell: tested with shapely `hausdorff_distance(..., densify=0.05)`. +- *Deterministic*: serial, ordered by (deviation, id); the same input gives + the same bits. +- Not guaranteed: the fewest vertices (it is greedy), or that every + catchment node is inside (a node within the tolerance of the boundary may + fall either side). + +**The tolerance**: `--outline-tolerance METRES`, default twice the DEM's cell +(20 m on DTM10). Default (main session / @architect, 2026-09-29), for Ola to +confirm. Divides are not known better than a cell or two (ours and NVE's disagree by +about 1.5 % of the area in nodes). The rough gauge above gives about 900 +vertices for Bygdin at 20 m, from 17,800. + +**Locality and parallelism** (the geometry skill asks): every collapse is +local (four vertices, a grid query); the heap is global and the loop serial. +At about 18,000 vertices that is milliseconds. It parallelises later by +cutting the ring into pieces fixed at their ends, if Glomma needs it. + +### The command + +``` +rasputin catchment --dem PATH [--dem PATH ...] --seed X Y [--seed-crs CRS] + [--lakes PATH [--lakes-layer NAME]] + [--outline-tolerance METRES] --out FILE.geojson + [--out-parent DIR] +``` + +- `--dem` as `mesh` has it: one directory or several files (15a). +- `--seed X Y`, two numbers, like `--bbox`'s four. `--seed-crs` is anything + pyproj reads, default `EPSG:4326` (so `LON LAT`), moved into the DEM's CRS + by `crs.reprojector` (15b's one `from_crs` site, `always_xy`). All tiles in + one CRS, or refuse, as 15b does. +- `--lakes`: `.gpkg` (with `--lakes-layer` when it has more than one feature + table, as the CORINE extract does), `.geojson` or `.json`. Its CRS is the + file's own (GeoPackage `srs_id`, GeoJSON `crs` member or WGS 84). +- `--out`: `.geojson` or `.json`, resolved and checked by the same + `_destination` as `mesh`. Written: a `FeatureCollection` of one `Feature`, + the reduced polygon in the DEM's CRS, with a `crs` member naming it + (`EPSG:25833`), which `domain.read_domain` already reads. Properties: seed + point and CRS, tolerance, node count, fine and reduced vertex counts, fine + and reduced areas in m², the windows' sizes. Coordinates written with + `repr` precision, so the area survives the round trip. +- `--outline-tolerance 0` writes the fine outline (collinear vertices only + removed), for comparison. +- stderr, one line each: every window (box, nodes, flood seconds, and + "grown: touches north" or "contained"); the seed (lake polygon with its + area and seed-node count, or the pour node); the catchment (nodes, node + area); the fine outline (vertices, area, rings dropped, holes filled); the + reduced outline (vertices, area, the difference in m² and relative, the + tolerance, seconds). Areas in km² with enough digits to see the + difference. +- Every refusal is a non-zero exit and no file, as `mesh`'s are. + +`catchment.delineate` is blocking (the C++ calls release the GIL); an async +caller runs it in `asyncio.to_thread`, as 16b's `open_features` is used. The +request and the result are frozen; the CLI is the only place with paths. + +### New and changed files + +| File | What | Production lines (estimate) | +|---|---|---| +| `include/terrain/hydrology/upstream.hpp` | the flood, `UpstreamOutcome` | 110 | +| `include/terrain/vector_simplify/area_collapse.hpp` | the reduction, `ReduceOutcome`, edge grid | 260 | +| `bindings/core.cpp` | `upstream`, `reduce_ring`, two outcome classes | 80 | +| `src_python/tin_engine/_core.pyi` | their stubs | 30 | +| `src_python/tin_engine/outline.py` | the tracer | 70 | +| `src_python/tin_engine/catchment.py` | request, seed, window loop, result | 170 | +| `src_python/tin_engine/dem_input.py` | repository helper split out | 10 | +| `src_python/tin_engine/cli.py` | `catchment` command, report, writer | 100 | +| `project_structure.md` | the two C++ modules and two Python modules | docs | + +About 830 lines, over the 700 ceiling (CLAUDE.md §2), so two PRs on this +branch. + +### The PR split + +- **PR 1, the fine catchment** (about 480 lines): `upstream.hpp` and its + binding, `outline.py`, `catchment.py`, the repository helper, and the + `catchment` command writing the fine outline. It already answers "what is + Bygdin's catchment" and can be compared against NVE. Red, green, review. +- **PR 2, the reduction** (about 350 lines): `area_collapse.hpp` and its + binding, and the command reducing by default with + `--outline-tolerance`. Red, green, review, then the Bygdin acceptance run. + +Both go on `increment22-autocatchment`, PR 2 stacked on PR 1. Neither touches +refine or mesh code, so the 1 m benchmark and scaling sweep (README, rule 2) +do not apply; the Bygdin run below is this increment's acceptance. + +## The red suites + +Lean: no throwaway implementations, no mutation round. Default (main session +/ @architect, 2026-09-29), for Ola to confirm; the brute-force oracle and the +invariants on random inputs are the defence. None is named invariant-critical +for mutation testing tonight. + +**PR 1, C++ (Catch2), `tests/cpp/unit/hydrology_upstream.cpp`:** + +- *Brute-force oracle*. For random DEMs with no interior pit and no ties, + built as z = 10 x (steps to the nearest edge) + a distinct fraction per + node, every node drains to its lowest neighbour, and the flood's catchment + of a random seed set must equal the set of nodes whose descent reaches it, + node for node. (On such a DEM the flood pops in global z order, so each + node is flooded by its lowest neighbour; the oracle is exact, not + approximate.) A few hundred small grids, seeds of one node and of blobs. +- *Depressions*: a closed bowl whose spill leads into the seed's valley is + wholly in; one that spills elsewhere is wholly out. +- *A flat lake on a plateau*: the lake nodes are the seed; a flat shelf + beside it that spills into it is in; a shelf that spills away is out. +- *V-valley*: two planes meeting at a valley whose outlet is the seed; the + catchment is the valley's side slopes up to the ridge lines, exactly. +- *Edge and NoData*: in-nodes on the window's edge set `touches_edge`; + an in-node beside NoData sets `touches_nodata`; NoData nodes are never in; + a seed on the edge is allowed; the bounds are the in-nodes' bounds. +- *Invariants on random DEMs with pits and flats*: seeds are in; every + in-node has an 8-path to a seed through in-nodes; running twice gives the + same mask. + +**PR 1, Python:** + +- `test_outline.py` (no `_core`): a single in-node gives one diamond of 4 + vertices and area half a cell; a 2 x 2 block; an L; a diagonal pair is one + ring (the saddle rule); a ring of nodes gives an outer ring and a clockwise + hole; random masks: every ring valid and simple, rings pairwise disjoint, + every in-node inside an odd number of rings and every out-node an even + number; mask on the window's edge closes (padding). +- `test_catchment.py` with an in-memory repository (as + `mosaic_fixtures.py` builds): a lake polygon from GeoJSON seeds the nodes + inside it; the point in no polygon, or in two, is refused; `--seed-crs` + 4326 lands on the same node as the DEM-CRS point; a catchment larger than + the first window grows and ends equal, node for node, to one flood over the + whole raster; reaching the data's edge, or NoData, is refused naming the + side; the memory cap refuses (monkeypatched `physical_memory`); holes are + filled and extra rings dropped, and the counts say so. +- `test_cli_catchment.py`: a synthetic tiled DEM directory and a GeoJSON lake; + the output reads back through `read_domain` in the DEM's CRS; `rasputin + mesh --dem ... --domain out.geojson --tolerance ...` succeeds on it; stderr + carries the window, fine vertex count and fine area lines; a wrong `--out` + suffix and `--lakes-layer` without `--lakes` are refused. + +**PR 2, C++ and Python:** + +- `tests/cpp/unit/area_collapse.cpp`: status for a negative, NaN or infinite + tolerance, a clockwise ring, fewer than 4 vertices; tolerance 0 removes + only collinear vertices; a traced rectangle of nodes reduces to at most 8 + vertices at a tolerance of one cell, with its area unchanged; a keep-point + half a cell inside a notch that a collapse would cut off stays inside; a + thin corridor (two long sides 1.5 cells apart) at a tolerance of 5 cells + must not cross itself; the same input twice gives equal bits. +- `test_core_reduce.py`, on rings traced from random blobs and from the + PR 1 synthetic catchments: the five guarantees above (area, simple, seed + inside, tolerance at vertices, Hausdorff within tolerance plus a cell), and + determinism. +- `test_cli_catchment.py` gains: the default reduces; the reduced line + reports vertices, area and the difference; `--outline-tolerance 0` gives + the fine ring. + +## Acceptance: Bygdin, end to end + +Run by `@perf` after PR 2 is green, recorded under +`docs/benchmarks/2026-09-29/bygdin/README.md` with the commands, the commit, +`pmset -g batt`, and the numbers: + +1. `rasputin catchment` as at the top. Record every window (box, nodes, + seconds), node count, fine and reduced vertex counts, fine and reduced + areas and their difference, the times of flood, trace and reduction. +2. **Sanity on the area**: against NVE's 305.54 km² (delfelt 1187, the URL + above). Expected within 2 %; the prototype gave -0.21 %. Outside 2 % is a + finding to explain before merge, not a pass. Also the node overlap with + NVE's polygon (fetched by the URL, not committed), both ways. +3. `rasputin mesh --dem ... --domain bygdin.geojson --tolerance 1` and + `--tolerance 10`, each with `--stats`: triangles, vertices, time. And once + at `--tolerance 10` with the fine outline (`--outline-tolerance 0`), to + show what the reduction saves, which is Ola's first requirement. +4. A picture is optional; the `.vtk` opens in ParaView. + +## What each persona reads + +`@tester` and `@developer`: this file, then `docs/increments/15-dem-mosaic.md` +(the plan and the re-plan loop) and `docs/increments/16-domain-polygon.md` +(what `--domain` accepts). `@developer` also reads `noding/intersect.hpp` for +`classify` and `raster/view.hpp`. From 6a700761510db706c7cf39902464513c3eb7e2c9 Mon Sep 17 00:00:00 2001 From: Ola Skavhaug Date: Tue, 29 Sep 2026 01:10:16 +0200 Subject: [PATCH 05/14] Retrospective agenda: a design-step prototype Co-Authored-By: Claude Opus 5.5 --- docs/retrospectives/next.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/docs/retrospectives/next.md b/docs/retrospectives/next.md index f3c2f1ad..e59d0c00 100644 --- a/docs/retrospectives/next.md +++ b/docs/retrospectives/next.md @@ -58,3 +58,7 @@ bug without a red test. `@perf` declined, citing its own rule (a failing test in `tests/python/test_bench.py` first). The brief crossed the line, the persona held it. The fix (create the mesh's parent directory for a `--label` with `/`) waits for a `@tester` red step. + +2026-09-29, night: `@architect` built an uncommitted Python prototype of +the 22 design to check its area against NVE (304.91 against 305.54 km²). +Useful evidence, but it is implementation work in the design step. From 1e1b3bb3735d0a53df6aaf0fe01a38fb3aa732e8 Mon Sep 17 00:00:00 2001 From: Ola Skavhaug Date: Tue, 29 Sep 2026 01:30:07 +0200 Subject: [PATCH 06/14] 22 PR1 red: the fine catchment The red suite for increment 22 PR 1 (docs/increments/22-auto-catchment.md, "The red suites", PR 1): the flood, the outline, the window loop and the `rasputin catchment` command. No production code. - tests/cpp/unit/test_hydrology_upstream.cpp: upstream equals lowest-neighbour descent node for node on 300 random no-pit, no-tie DEMs (an oracle in the file); depressions, a flat lake with shelves, a V-valley (exact); NoData never in; seed on the edge; wrong seed size; 1-node rasters; invariants on 400 random DEMs with pits, flats and NoData (seeds in, 8-path to a seed, bounds, determinism, catchments of a union are the union). Registered behind if(EXISTS hydrology/upstream.hpp), increment 12's scaffold: REMOVE AT GREEN. - test_core_upstream.py: the binding (mask uint8 (rows, cols), bounds, flags, ValueError on a seed of another shape). - test_outline.py: diamond, octagon, L, saddle, pinch, holes clockwise, and 40 random masks against the guarantee and a square-by-square area oracle. - test_catchment.py (in-memory repository): lake seeds, refusals for no lake and two lakes, WGS 84 and EPSG:3035 inputs, the grown window equal to one flood over the whole raster, the edge (north, east), NoData and memory refusals, holes filled and pieces dropped with counts, to_thread. - test_cli_catchment.py: round trip through read_domain and mesh --domain, stderr lines, --seed-crs default, GeoJSON and GeoPackage lakes, refusals, and Bygdin on Ola's data (skipped without it) within 2 % of 305.54 km^2. Red because include/terrain/hydrology/upstream.hpp, _core.upstream, tin_engine.outline, tin_engine.catchment and the catchment command do not exist. Every fixture's expected value was checked in the scratchpad with the design's flood written out, not committed. The main build exits 0, ctest 803/803, and the rest of pytest, ruff, format and mypy are green. Names the design leaves open are chosen in each file's docstring. Design gap: edge and NoData-adjacent nodes are outlets and labelled in only when seeds, so touches_edge/touches_nodata as worded fire only for seeds; the C++ suite pins only what any repair keeps, and the Python suite pins the refusals by behaviour. Co-Authored-By: Claude Opus 5.5 --- tests/cpp/CMakeLists.txt | 22 + tests/cpp/unit/test_hydrology_upstream.cpp | 589 +++++++++++++++++++++ tests/python/catchment_fixtures.py | 139 +++++ tests/python/test_catchment.py | 352 ++++++++++++ tests/python/test_cli_catchment.py | 283 ++++++++++ tests/python/test_core_upstream.py | 122 +++++ tests/python/test_outline.py | 253 +++++++++ 7 files changed, 1760 insertions(+) create mode 100644 tests/cpp/unit/test_hydrology_upstream.cpp create mode 100644 tests/python/catchment_fixtures.py create mode 100644 tests/python/test_catchment.py create mode 100644 tests/python/test_cli_catchment.py create mode 100644 tests/python/test_core_upstream.py create mode 100644 tests/python/test_outline.py diff --git a/tests/cpp/CMakeLists.txt b/tests/cpp/CMakeLists.txt index fb00cc92..f8c91a8d 100644 --- a/tests/cpp/CMakeLists.txt +++ b/tests/cpp/CMakeLists.txt @@ -311,3 +311,25 @@ add_terrain_backend_test(test_mesh_lattice_incircle unit/test_mesh_lattice_incir # threads, so it is not in the TSan job. add_terrain_backend_test(prop_noding_verifier_sweep property/prop_noding_verifier_sweep.cpp) + +# The fine catchment's flood, increment 22 PR 1 +# (docs/increments/22-auto-catchment.md, "Flow and membership: one flood" and +# "The red suites"). Not named invariant-critical; the descent oracle and the +# invariants on random DEMs are the defence. It starts no threads. +# +# RED UNTIL hydrology/upstream.hpp EXISTS (increment 12's scaffold, 39dba00): +# while the header is absent the target is EXCLUDE_FROM_ALL and not registered +# with ctest, so `cmake --build build` and `ctest` stay green; see the red with +# `cmake --build build --target test_hydrology_upstream`. Once the header lands, +# re-running `cmake -S . -B build` registers the suite with no edit here. +# REMOVE AT GREEN: replace the if/else with the plain add_terrain_test line. +if(EXISTS "${PROJECT_SOURCE_DIR}/include/terrain/hydrology/upstream.hpp") + add_terrain_test(test_hydrology_upstream unit/test_hydrology_upstream.cpp) +else() + message(STATUS "test_hydrology_upstream: include/terrain/hydrology/upstream.hpp absent; " + "suite built only on request (increment 22 red step)") + add_executable(test_hydrology_upstream EXCLUDE_FROM_ALL unit/test_hydrology_upstream.cpp) + target_link_libraries(test_hydrology_upstream PRIVATE + terrain_headers terrain_test_support Catch2::Catch2WithMain) + target_compile_options(test_hydrology_upstream PRIVATE -Wall -Wextra -Wpedantic -Werror) +endif() diff --git a/tests/cpp/unit/test_hydrology_upstream.cpp b/tests/cpp/unit/test_hydrology_upstream.cpp new file mode 100644 index 00000000..97e20b46 --- /dev/null +++ b/tests/cpp/unit/test_hydrology_upstream.cpp @@ -0,0 +1,589 @@ +// Increment 22, PR 1 (docs/increments/22-auto-catchment.md, "Flow and +// membership: one flood" and "The red suites", PR 1 C++): the flood that +// labels every DEM node draining into a seed set. +// +// The design names this file tests/cpp/unit/hydrology_upstream.cpp; it carries +// the test_ prefix every other suite in tests/cpp/unit/ has. +// +// Interface assumed (the design fixes the header, the template on the +// RasterSource concept, and what is returned; the names are chosen here and +// stated in the handback): +// +// #include +// namespace terrain::hydrology { +// struct UpstreamOutcome { +// std::vector mask; // row-major, geometry().size() bytes: +// // 1 for a node in the catchment, 0 for +// // every other (out, unreached, NoData) +// std::size_t nodes_in; // the number of 1s +// std::size_t row_min, row_max, col_min, col_max; // the in-nodes' bounds, +// // inclusive; meaningful when nodes_in > 0 +// bool touches_edge; // see below +// bool touches_nodata; +// }; +// template +// UpstreamOutcome upstream(const R& z, std::span seed); +// } +// +// `seed` is row-major, one byte per node, non-zero for a seed. A span whose +// size is not geometry().size() is std::invalid_argument (the binding's +// ValueError). A NoData seed is ignored: NoData is never in. +// +// THE FLAGS ARE PINNED ONLY WHERE THE DESIGN'S WORDS AND ANY REPAIR OF THEM +// AGREE. As written, every edge node and every NoData-adjacent node is an +// outlet and "outlets are labelled in if they are seeds", so "an in-node on +// the window's edge" can only be a seed; a catchment truncated by the window +// shows as in-nodes one node inside the edge, not on it. That is reported as a +// design gap. This suite pins the cases true under either reading: a seed on +// the edge sets touches_edge; a catchment with no in-node within one node of +// the edge (or within two of NoData) sets neither flag; an in-node on the edge +// (beside NoData) implies the flag. The Python window suite pins the +// behaviour that matters, truncation refused, whatever the flag means. +// +// THE ORACLE for the no-pit, no-tie family is lowest-neighbour descent, +// computed here and sharing nothing with the header: on +// z = 10 * (steps to the nearest edge) + a distinct fraction in [0, 1) +// every interior node has a strictly lower neighbour, levels equal z, the +// flood pops in global z order, and so each node is flooded by its lowest +// neighbour. The flood's catchment must equal descent's node for node. + +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using terrain::hydrology::upstream; +using terrain::hydrology::UpstreamOutcome; +using terrain::raster::CellIndex; +using terrain::raster::Raster; +using terrain::raster::RasterGeometry; +using terrain::raster::RasterView; + +namespace { + +constexpr float kNaN = std::numeric_limits::quiet_NaN(); +constexpr float kSentinel = -9999.0F; + +using Mask = std::vector; + +RasterGeometry grid(std::size_t rows, std::size_t cols) { + // UTM-scale origin, unequal spacing: nothing here depends on either, and a + // flood that did would show it. + return RasterGeometry{500000.0, 6600000.0, 10.0, 5.0, cols, rows}; +} + +struct Dem { + std::size_t rows; + std::size_t cols; + std::vector z; + + [[nodiscard]] std::size_t at(std::size_t r, std::size_t c) const { return r * cols + c; } + [[nodiscard]] Raster raster(std::optional nodata = std::nullopt) const { + return Raster{grid(rows, cols), z, nodata}; + } +}; + +Mask none(const Dem& d) { return Mask(d.rows * d.cols, 0); } + +UpstreamOutcome flood(const Dem& d, const Mask& seed, std::optional nodata = std::nullopt) { + return upstream(d.raster(nodata), std::span{seed}); +} + +bool invalid(const Dem& d, std::size_t i, std::optional nodata) { + const float v = d.z[i]; + return std::isnan(v) || (nodata && v == *nodata); +} + +template +void each_neighbour(const Dem& d, std::size_t r, std::size_t c, F&& f) { + for (int dr = -1; dr <= 1; ++dr) + for (int dc = -1; dc <= 1; ++dc) { + if (dr == 0 && dc == 0) + continue; + const auto rr = static_cast(r) + dr; + const auto cc = static_cast(c) + dc; + if (rr < 0 || cc < 0 || rr >= static_cast(d.rows) + || cc >= static_cast(d.cols)) + continue; + f(static_cast(rr), static_cast(cc)); + } +} + +bool on_edge(const Dem& d, std::size_t r, std::size_t c) { + return r == 0 || c == 0 || r + 1 == d.rows || c + 1 == d.cols; +} + +// ---- the no-pit, no-tie family and its descent oracle ---------------------- + +Dem no_pit_dem(std::size_t rows, std::size_t cols, std::mt19937& rng) { + Dem d{rows, cols, std::vector(rows * cols)}; + std::vector order(rows * cols); + std::iota(order.begin(), order.end(), std::size_t{0}); + std::shuffle(order.begin(), order.end(), rng); + const auto n = static_cast(rows * cols); + for (std::size_t r = 0; r < rows; ++r) + for (std::size_t c = 0; c < cols; ++c) { + const std::size_t steps = std::min({r, c, rows - 1 - r, cols - 1 - c}); + // A distinct fraction per node, k / n with k a permutation: exact in + // float32 for every grid here (n <= 900, steps <= 14). + d.z[d.at(r, c)] = 10.0F * static_cast(steps) + + static_cast(order[d.at(r, c)]) / n; + } + return d; +} + +// In iff the descent chain from the node (the node itself included) meets a +// seed. Edge nodes end their chain: they are outlets. +Mask descent_oracle(const Dem& d, const Mask& seed) { + const std::size_t n = d.rows * d.cols; + std::vector next(n, n); + for (std::size_t r = 0; r < d.rows; ++r) + for (std::size_t c = 0; c < d.cols; ++c) { + if (on_edge(d, r, c)) + continue; + float best = std::numeric_limits::infinity(); + each_neighbour(d, r, c, [&](std::size_t rr, std::size_t cc) { + if (d.z[d.at(rr, cc)] < best) { + best = d.z[d.at(rr, cc)]; + next[d.at(r, c)] = d.at(rr, cc); + } + }); + REQUIRE(best < d.z[d.at(r, c)]); // the family's premise: no pit + } + Mask in(n, 0); + for (std::size_t i = 0; i < n; ++i) + for (std::size_t j = i; j != n; j = next[j]) + if (seed[j]) { + in[i] = 1; + break; + } + return in; +} + +// ---- invariants any outcome must satisfy ----------------------------------- + +void check_invariants(const Dem& d, const Mask& seed, const UpstreamOutcome& out, + std::optional nodata = std::nullopt) { + const std::size_t n = d.rows * d.cols; + REQUIRE(out.mask.size() == n); + std::size_t count = 0; + std::size_t r0 = d.rows, r1 = 0, c0 = d.cols, c1 = 0; + bool edge_in = false, near_edge_in = false, beside_nodata_in = false, near_nodata_in = false; + for (std::size_t r = 0; r < d.rows; ++r) + for (std::size_t c = 0; c < d.cols; ++c) { + const std::size_t i = d.at(r, c); + REQUIRE((out.mask[i] == 0 || out.mask[i] == 1)); + if (invalid(d, i, nodata)) { + REQUIRE(out.mask[i] == 0); // NoData is never in, seed or not + continue; + } + if (seed[i]) + REQUIRE(out.mask[i] == 1); // every valid seed is in + if (!out.mask[i]) + continue; + ++count; + r0 = std::min(r0, r); + r1 = std::max(r1, r); + c0 = std::min(c0, c); + c1 = std::max(c1, c); + edge_in = edge_in || on_edge(d, r, c); + near_edge_in = near_edge_in || r <= 1 || c <= 1 || r + 2 >= d.rows || c + 2 >= d.cols; + for (std::size_t rr = r > 2 ? r - 2 : 0; rr <= std::min(r + 2, d.rows - 1); ++rr) + for (std::size_t cc = c > 2 ? c - 2 : 0; cc <= std::min(c + 2, d.cols - 1); ++cc) + if (invalid(d, d.at(rr, cc), nodata)) { + near_nodata_in = true; + const bool adjacent = (rr + 1 >= r && rr <= r + 1 && cc + 1 >= c && cc <= c + 1); + beside_nodata_in = beside_nodata_in || adjacent; + } + } + CHECK(out.nodes_in == count); + if (count > 0) { + CHECK(out.row_min == r0); + CHECK(out.row_max == r1); + CHECK(out.col_min == c0); + CHECK(out.col_max == c1); + } + if (edge_in) + CHECK(out.touches_edge); + if (!near_edge_in) + CHECK_FALSE(out.touches_edge); + if (beside_nodata_in) + CHECK(out.touches_nodata); + if (!near_nodata_in) + CHECK_FALSE(out.touches_nodata); + + // Every in-node has an 8-path through in-nodes to a seed. + Mask reached(n, 0); + std::deque queue; + for (std::size_t i = 0; i < n; ++i) + if (seed[i] && out.mask[i]) { + reached[i] = 1; + queue.push_back(i); + } + while (!queue.empty()) { + const std::size_t i = queue.front(); + queue.pop_front(); + each_neighbour(d, i / d.cols, i % d.cols, [&](std::size_t rr, std::size_t cc) { + const std::size_t j = d.at(rr, cc); + if (out.mask[j] && !reached[j]) { + reached[j] = 1; + queue.push_back(j); + } + }); + } + CHECK(reached == out.mask); +} + +Mask random_seeds(std::size_t n, double density, std::mt19937& rng) { + std::bernoulli_distribution pick(density); + Mask seed(n, 0); + for (auto& s : seed) + s = pick(rng) ? 1 : 0; + return seed; +} + +Mask blob(const Dem& d, std::mt19937& rng) { + std::uniform_int_distribution rr(0, d.rows - 1), cc(0, d.cols - 1); + std::size_t a = rr(rng), b = rr(rng), e = cc(rng), f = cc(rng); + Mask seed = none(d); + for (std::size_t r = std::min(a, b); r <= std::max(a, b); ++r) + for (std::size_t c = std::min(e, f); c <= std::max(e, f); ++c) + seed[d.at(r, c)] = 1; + return seed; +} + +} // namespace + +TEST_CASE("the flood equals lowest-neighbour descent on DEMs with no pit and no tie", + "[hydrology][upstream][oracle]") { + std::mt19937 rng{22001}; + std::uniform_int_distribution size(3, 30); + int cases = 0, nontrivial = 0; + for (int trial = 0; trial < 300; ++trial) { + const Dem d = no_pit_dem(size(rng), size(rng), rng); + const std::size_t n = d.rows * d.cols; + Mask seed = none(d); + switch (trial % 3) { + case 0: seed[std::uniform_int_distribution(0, n - 1)(rng)] = 1; break; + case 1: seed = blob(d, rng); break; + default: seed = random_seeds(n, 0.05, rng); break; + } + const Mask expected = descent_oracle(d, seed); + const UpstreamOutcome out = flood(d, seed); + INFO("trial " << trial << ", " << d.rows << " x " << d.cols); + REQUIRE(out.mask == expected); + check_invariants(d, seed, out); + ++cases; + const auto seeds = static_cast(std::count(seed.begin(), seed.end(), 1)); + nontrivial += out.nodes_in > seeds ? 1 : 0; + } + // The oracle must be able to disagree: most catchments are larger than + // their seed set, so a flood that returned the seeds alone would fail. + CHECK(cases == 300); + CHECK(nontrivial > 150); +} + +TEST_CASE("a closed depression is wholly in when it spills towards the seed, wholly out otherwise", + "[hydrology][upstream][depression]") { + // Every edge node is a high wall (100) except two outlets at 0: the seed + // (7, 0) on the west edge, and a plain outlet (3, 14) on the east. The + // interior is the lower of two funnels, one towards each. + Dem d{15, 15, std::vector(225)}; + for (std::size_t r = 0; r < 15; ++r) + for (std::size_t c = 0; c < 15; ++c) { + const float dr7 = std::abs(static_cast(r) - 7.0F); + const float dr3 = std::abs(static_cast(r) - 3.0F); + const float west = 10.0F + 2.0F * dr7 + static_cast(c); + const float east = 10.0F + 2.0F * dr3 + static_cast(14 - c) + 0.5F; + d.z[d.at(r, c)] = on_edge(d, r, c) ? 100.0F : std::min(west, east); + } + d.z[d.at(7, 0)] = 0.0F; + d.z[d.at(3, 14)] = 0.0F; + // Bowl A, 2 x 2 at rows 10-11, cols 4-5, in the west funnel (west 20-23, + // east above 33): a pit well below its rim. + for (std::size_t r : {10U, 11U}) + for (std::size_t c : {4U, 5U}) + d.z[d.at(r, c)] = 1.0F; + // Bowl B, 2 x 3 at rows 2-3, cols 10-12, in the east funnel. + for (std::size_t r : {2U, 3U}) + for (std::size_t c : {10U, 11U, 12U}) + d.z[d.at(r, c)] = 1.0F; + + Mask seed = none(d); + seed[d.at(7, 0)] = 1; + const UpstreamOutcome out = flood(d, seed); + check_invariants(d, seed, out); + for (std::size_t r : {10U, 11U}) + for (std::size_t c : {4U, 5U}) + CHECK(out.mask[d.at(r, c)] == 1); + for (std::size_t r : {2U, 3U}) + for (std::size_t c : {10U, 11U, 12U}) + CHECK(out.mask[d.at(r, c)] == 0); + CHECK(out.mask[d.at(3, 14)] == 0); // the other outlet + CHECK(out.touches_edge); // the seed is on the west edge +} + +TEST_CASE("a flat lake with shelves: a shelf that spills into it is in, one that spills away is out", + "[hydrology][upstream][flat]") { + // 24 x 24. The lake is the flat block rows 10-13, cols 10-13, at 50. Around + // it, by Chebyshev distance k from the block: a bowl rising to a rim at + // k = 3 (70, 80, 90), then a slope falling to the edges (90 - 5 (k - 3), + // 55 on the edge). The lake's outlet is a channel at 40 along row 11 from + // the lake to the east edge, so the lake is not itself a depression. + Dem d{24, 24, std::vector(576)}; + const auto k_of = [](std::size_t r, std::size_t c) { + const auto axis = [](std::size_t i) { + return i < 10 ? 10 - i : (i > 13 ? i - 13 : std::size_t{0}); + }; + return std::max(axis(r), axis(c)); + }; + for (std::size_t r = 0; r < 24; ++r) + for (std::size_t c = 0; c < 24; ++c) { + const std::size_t k = k_of(r, c); + d.z[d.at(r, c)] = k == 0 ? 50.0F + : k <= 3 ? 60.0F + 10.0F * static_cast(k) + : 90.0F - 5.0F * static_cast(k - 3); + } + for (std::size_t c = 14; c < 24; ++c) + d.z[d.at(11, c)] = 40.0F; + // Shelf N: a flat at 60 outside the rim (k = 4, row 6, cols 8-15), closed + // by the slope (80 at row 5) except for a channel at 55 through the rim + // (col 12, rows 7-9) down to the lake: it spills into the lake. + for (std::size_t c = 8; c <= 15; ++c) + d.z[d.at(6, c)] = 60.0F; + for (std::size_t r = 7; r <= 9; ++r) + d.z[d.at(r, 12)] = 55.0F; + // Shelf S: a flat at 60 outside the rim on the south (k = 5, row 18, cols + // 8-15), whose lowest way out is outwards (k = 6, 75), not inwards (k = 4, + // 85): it spills away. + for (std::size_t c = 8; c <= 15; ++c) + d.z[d.at(18, c)] = 60.0F; + + Mask seed = none(d); + for (std::size_t r = 10; r <= 13; ++r) + for (std::size_t c = 10; c <= 13; ++c) + seed[d.at(r, c)] = 1; + const UpstreamOutcome out = flood(d, seed); + check_invariants(d, seed, out); + + for (std::size_t c = 8; c <= 15; ++c) { + CHECK(out.mask[d.at(6, c)] == 1); + CHECK(out.mask[d.at(18, c)] == 0); + } + for (std::size_t r = 7; r <= 9; ++r) + CHECK(out.mask[d.at(r, 12)] == 1); + // The outlet channel is downstream of the lake: out. + for (std::size_t c = 14; c < 24; ++c) + CHECK(out.mask[d.at(11, c)] == 0); + // The bowl inside the rim drains to the lake, away from the channel's mouth. + for (std::size_t r = 0; r < 24; ++r) + for (std::size_t c = 0; c <= 12; ++c) + if (k_of(r, c) <= 2) + CHECK(out.mask[d.at(r, c)] == 1); + // Far from every edge and with no NoData: no flag. + CHECK_FALSE(out.touches_edge); + CHECK_FALSE(out.touches_nodata); +} + +TEST_CASE("a V-valley: the catchment is its side slopes up to the ridge lines, exactly", + "[hydrology][upstream][valley]") { + // 21 x 15. Valley along row 10, falling west at b = 0.37 per column, sides + // rising at a = 1 per row to ridges at rows 4 and 16, outer slopes falling + // at 2.9 per row to the north and south edges. Column 0 is a wall (1000) + // except the seed (10, 0) at -1, the valley's outlet. + constexpr std::size_t rows = 21, cols = 15, rv = 10, w = 6; + Dem d{rows, cols, std::vector(rows * cols)}; + for (std::size_t r = 0; r < rows; ++r) + for (std::size_t c = 0; c < cols; ++c) { + const auto off = static_cast(r > rv ? r - rv : rv - r); + const float side = off <= static_cast(w) + ? off + : static_cast(w) - 2.9F * (off - static_cast(w)); + d.z[d.at(r, c)] = c == 0 ? 1000.0F : 0.37F * static_cast(c) + side; + } + d.z[d.at(rv, 0)] = -1.0F; + + Mask seed = none(d); + seed[d.at(rv, 0)] = 1; + const UpstreamOutcome out = flood(d, seed); + check_invariants(d, seed, out); + + Mask expected = none(d); + expected[d.at(rv, 0)] = 1; + for (std::size_t r = rv - (w - 1); r <= rv + (w - 1); ++r) + for (std::size_t c = 1; c + 1 < cols; ++c) + expected[d.at(r, c)] = 1; + CHECK(out.mask == expected); + CHECK(out.nodes_in == 11 * 13 + 1); + CHECK(out.row_min == 5); + CHECK(out.row_max == 15); + CHECK(out.col_min == 0); + CHECK(out.col_max == 13); + CHECK(out.touches_edge); +} + +TEST_CASE("NoData is never in, a NoData seed is ignored, and a seed beside NoData sets the flag", + "[hydrology][upstream][nodata]") { + std::mt19937 rng{22002}; + Dem d = no_pit_dem(12, 12, rng); + d.z[d.at(5, 5)] = kNaN; + d.z[d.at(8, 3)] = kSentinel; + + SECTION("NaN and the sentinel, as seeds, are not in") { + Mask seed = none(d); + seed[d.at(5, 5)] = 1; + seed[d.at(8, 3)] = 1; + const UpstreamOutcome out = flood(d, seed, kSentinel); + CHECK(out.nodes_in == 0); + CHECK(std::all_of(out.mask.begin(), out.mask.end(), [](auto m) { return m == 0; })); + CHECK_FALSE(out.touches_nodata); + CHECK_FALSE(out.touches_edge); + } + SECTION("without a sentinel, -9999 is a value like any other") { + Mask seed = none(d); + seed[d.at(8, 3)] = 1; + const UpstreamOutcome out = flood(d, seed); + CHECK(out.mask[d.at(8, 3)] == 1); + } + SECTION("a seed beside NoData is in and sets touches_nodata") { + Mask seed = none(d); + seed[d.at(5, 6)] = 1; + const UpstreamOutcome out = flood(d, seed, kSentinel); + CHECK(out.mask[d.at(5, 6)] == 1); + CHECK(out.mask[d.at(5, 5)] == 0); + CHECK(out.touches_nodata); + check_invariants(d, seed, out, kSentinel); + } + SECTION("every valid node seeded: every valid node in, no NoData in") { + Mask seed(d.rows * d.cols, 1); + const UpstreamOutcome out = flood(d, seed, kSentinel); + CHECK(out.nodes_in == d.rows * d.cols - 2); + check_invariants(d, seed, out, kSentinel); + } +} + +TEST_CASE("a seed on the window's edge is allowed and sets touches_edge; the bounds are the in-nodes'", + "[hydrology][upstream][edge]") { + std::mt19937 rng{22003}; + const Dem d = no_pit_dem(9, 11, rng); + Mask seed = none(d); + seed[d.at(0, 4)] = 1; + const UpstreamOutcome out = flood(d, seed); + CHECK(out.mask[d.at(0, 4)] == 1); + CHECK(out.touches_edge); + CHECK(out.row_min == 0); + check_invariants(d, seed, out); + // An edge node is an outlet: nothing drains through it, so on this family + // the seed on the edge is its own whole catchment. + CHECK(out.mask == descent_oracle(d, seed)); +} + +TEST_CASE("no seed: an empty catchment and no flag", "[hydrology][upstream][empty]") { + std::mt19937 rng{22004}; + const Dem d = no_pit_dem(7, 6, rng); + const UpstreamOutcome out = flood(d, none(d)); + CHECK(out.nodes_in == 0); + CHECK(out.mask == none(d)); + CHECK_FALSE(out.touches_edge); + CHECK_FALSE(out.touches_nodata); +} + +TEST_CASE("a seed mask of the wrong size is std::invalid_argument", "[hydrology][upstream][refusal]") { + std::mt19937 rng{22005}; + const Dem d = no_pit_dem(5, 6, rng); + const Raster raster = d.raster(); + const Mask short_by_one(29, 0), long_by_one(31, 0), empty; + CHECK_THROWS_AS(upstream(raster, std::span{short_by_one}), + std::invalid_argument); + CHECK_THROWS_AS(upstream(raster, std::span{long_by_one}), + std::invalid_argument); + CHECK_THROWS_AS(upstream(raster, std::span{empty}), std::invalid_argument); +} + +TEST_CASE("one-row, one-column and single-node rasters: every node is an outlet, the seeds are the catchment", + "[hydrology][upstream][degenerate]") { + constexpr std::array, 4> shapes{{{1, 1}, {1, 7}, {7, 1}, {2, 2}}}; + for (const auto& [rows, cols] : shapes) { + Dem d{rows, cols, std::vector(rows * cols)}; + for (std::size_t i = 0; i < d.z.size(); ++i) + d.z[i] = static_cast(i % 3); // ties and a rise, deliberately + Mask seed = none(d); + seed[d.z.size() / 2] = 1; + const UpstreamOutcome out = flood(d, seed); + INFO(rows << " x " << cols); + CHECK(out.mask == seed); + CHECK(out.touches_edge); + } +} + +TEST_CASE("invariants on random DEMs with pits, flats and NoData", "[hydrology][upstream][property]") { + std::mt19937 rng{22006}; + std::uniform_int_distribution size(1, 24); + std::uniform_int_distribution level(0, 4); // few levels: flats, ties, pits + std::bernoulli_distribution hole(0.04); + for (int trial = 0; trial < 400; ++trial) { + Dem d{size(rng), size(rng), {}}; + d.z.resize(d.rows * d.cols); + const bool with_nodata = trial % 2 == 1; + for (auto& v : d.z) + v = with_nodata && hole(rng) ? (trial % 4 == 1 ? kNaN : kSentinel) + : static_cast(level(rng)); + const std::optional nodata = + with_nodata ? std::optional{kSentinel} : std::nullopt; + const std::size_t n = d.rows * d.cols; + const Mask s1 = random_seeds(n, 0.03, rng); + const Mask s2 = blob(d, rng); + Mask both(n, 0); + for (std::size_t i = 0; i < n; ++i) + both[i] = (s1[i] || s2[i]) ? 1 : 0; + + INFO("trial " << trial << ", " << d.rows << " x " << d.cols); + const UpstreamOutcome a = flood(d, s1, nodata); + const UpstreamOutcome b = flood(d, s2, nodata); + const UpstreamOutcome ab = flood(d, both, nodata); + check_invariants(d, s1, a, nodata); + check_invariants(d, s2, b, nodata); + check_invariants(d, both, ab, nodata); + + // Deterministic: the same input twice gives the same mask and counts. + const UpstreamOutcome again = flood(d, s1, nodata); + CHECK(again.mask == a.mask); + CHECK(again.nodes_in == a.nodes_in); + + // The flood's order does not depend on the labels, so catchments add: + // the catchment of a union of seed sets is the union of catchments. + Mask united(n, 0); + for (std::size_t i = 0; i < n; ++i) + united[i] = (a.mask[i] || b.mask[i]) ? 1 : 0; + CHECK(ab.mask == united); + } +} + +TEST_CASE("an owning Raster and a RasterView over the same buffer give the same flood", + "[hydrology][upstream][concept]") { + std::mt19937 rng{22007}; + const Dem d = no_pit_dem(13, 17, rng); + std::vector z64(d.z.begin(), d.z.end()); + const RasterView view{grid(d.rows, d.cols), z64.data(), std::nullopt}; + const Mask seed = blob(d, rng); + const UpstreamOutcome from_view = upstream(view, std::span{seed}); + const UpstreamOutcome from_raster = flood(d, seed); + CHECK(from_view.mask == from_raster.mask); + CHECK(from_view.nodes_in == from_raster.nodes_in); +} diff --git a/tests/python/catchment_fixtures.py b/tests/python/catchment_fixtures.py new file mode 100644 index 00000000..f8d06a0e --- /dev/null +++ b/tests/python/catchment_fixtures.py @@ -0,0 +1,139 @@ +"""Terrains, lakes and an in-memory repository for increment 22's suites. + +`docs/increments/22-auto-catchment.md`, "The red suites" (PR 1). Imports only +what earlier increments shipped, so a missing 22 module fails the tests that +use it, not the collection of this helper. + +Every terrain is on a 100 m lattice, so the window's 2000 m margin +(`WINDOW_MARGIN_M`) is 20 cells and a catchment can outgrow the first window +on a grid of a few hundred nodes. + +THE BOWL. An elliptical basin: z = 1000 (1 - |1 - b|) with +b = ((c - cc) / ax)^2 + ((r - rc) / ay)^2, rising from the centre to a rim at +b = 1 and falling outside it to the data's edges. A small tilt breaks the +mirror ties. The lake is a 7 x 7-node box round the centre, and its outlet is +a channel down column cc from the centre to the south edge, strictly +descending and below the lake: the lake is not a closed depression, as a real +lake with a river out of it is not. What drains into the lake is the part of +the basin that reaches it before the channel. +""" + +from __future__ import annotations + +from collections import deque +from collections.abc import Mapping + +import numpy as np +import numpy.typing as npt +import shapely +from shapely.geometry import Polygon, box + +from mosaic_fixtures import whole +from tin_engine.io.models import DemTile +from tin_engine.io.repository import TileFootprint + +X0 = 500_000.0 +Y0 = 6_600_000.0 +D = 100.0 +EPSG = "EPSG:25833" +MARGIN_CELLS = 20 # WINDOW_MARGIN_M / D + + +def lat(col: float, row: float) -> tuple[float, float]: + """A lattice position (col, row) as a point in EPSG:25833.""" + return (X0 + D * col, Y0 - D * row) + + +def bowl( + rows: int = 400, + cols: int = 200, + rc: int = 200, + cc: int = 100, + ax: float = 8.0, + ay: float = 55.0, +) -> npt.NDArray[np.float32]: + r, c = np.indices((rows, cols)).astype(np.float64) + b = ((c - cc) / ax) ** 2 + ((r - rc) / ay) ** 2 + z = 1000.0 * (1.0 - np.abs(1.0 - b)) + 0.0131 * c + 0.0077 * r + channel = np.arange(rows) >= rc + z[channel, cc] = -1.0 - 0.5 * (np.arange(rows)[channel] - rc) + return z.astype(np.float32) + + +def tile_of(z: npt.NDArray[np.floating], nodata: float | None = None) -> DemTile: + rows, cols = z.shape + return whole(rows, cols, array=z, x_min=X0, y_max=Y0, dx=D, dy=D, nodata=nodata) + + +def lake_box(rc: int = 200, cc: int = 100, half: float = 3.5) -> Polygon: + """The lake: a box round (rc, cc), its sides half a cell from any node.""" + x0, y0 = lat(cc - half, rc + half) + x1, y1 = lat(cc + half, rc - half) + return box(x0, y0, x1, y1) + + +class MemoryRepository: + """A `DemRepository` over tiles held in memory; records every load.""" + + def __init__(self, tiles: Mapping[str, DemTile]) -> None: + self.tiles = dict(tiles) + self.loads: list[str] = [] + + def footprints(self) -> tuple[TileFootprint, ...]: + return tuple( + TileFootprint(name=name, meta=tile.meta, dtype=tile.array.dtype) + for name, tile in sorted(self.tiles.items()) + ) + + def load(self, name: str) -> DemTile: + self.loads.append(name) + return self.tiles[name] + + +def node_grid(tile: DemTile) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]: + m = tile.meta + r, c = np.indices((m.rows, m.cols)).astype(np.float64) + return m.x_min + c * m.delta_x, m.y_max - r * m.delta_y + + +def seeds_in(tile: DemTile, lake: shapely.Geometry) -> npt.NDArray[np.uint8]: + """The design's seed mask: DEM nodes inside the lake, `shapely.contains_xy`.""" + x, y = node_grid(tile) + return shapely.contains_xy(lake, x, y).astype(np.uint8) + + +def full_flood(tile: DemTile, seed: npt.NDArray[np.uint8]) -> npt.NDArray[np.uint8]: + """One flood over the whole raster: the reference the window loop must equal.""" + from tin_engine._core import upstream + from tin_engine.raster import to_core + + return np.asarray(upstream(to_core(tile), seed).mask, dtype=np.uint8) + + +def filled(mask: npt.ArrayLike) -> npt.NDArray[np.uint8]: + """`mask` with its holes filled: every out-node the padding cannot reach + through 4-connected out-nodes becomes in (the (8, 4) pairing).""" + m = np.pad(np.asarray(mask) != 0, 1) + outside = np.zeros_like(m) + outside[0, 0] = True + queue = deque([(0, 0)]) + rows, cols = m.shape + while queue: + r, c = queue.popleft() + for rr, cc in ((r - 1, c), (r + 1, c), (r, c - 1), (r, c + 1)): + if 0 <= rr < rows and 0 <= cc < cols and not m[rr, cc] and not outside[rr, cc]: + outside[rr, cc] = True + queue.append((rr, cc)) + return (~outside[1:-1, 1:-1]).astype(np.uint8) + + +def placed( + mask: npt.ArrayLike, window_x_min: float, window_y_max: float, shape: tuple[int, int] +) -> npt.NDArray[np.uint8]: + """A window's mask on the whole raster's lattice (origin X0, Y0, spacing D).""" + m = np.asarray(mask, dtype=np.uint8) + r0 = round((Y0 - window_y_max) / D) + c0 = round((window_x_min - X0) / D) + out = np.zeros(shape, dtype=np.uint8) + out[r0 : r0 + m.shape[0], c0 : c0 + m.shape[1]] = m + return out diff --git a/tests/python/test_catchment.py b/tests/python/test_catchment.py new file mode 100644 index 00000000..e2aa8223 --- /dev/null +++ b/tests/python/test_catchment.py @@ -0,0 +1,352 @@ +"""`catchment.delineate`: seed, window loop and fine outline (increment 22, PR 1). + +`docs/increments/22-auto-catchment.md`, "Data flow", "The seed", "The window", +"The fine outline" and "The red suites" (PR 1, `test_catchment.py`). The DEM +is an in-memory repository (`catchment_fixtures.MemoryRepository`, as +`mosaic_fixtures` builds tiles); no file is read. + +Interface assumed (the design names `CatchmentRequest` (frozen Pydantic), +`delineate(request, repository) -> Catchment` (frozen), `CatchmentError`, +`WINDOW_MARGIN_M`; the fields are chosen here and stated in the handback): + +- `CatchmentRequest(seed=(x, y), seed_crs="EPSG:4326", lakes=None, + lakes_crs=None)`. `lakes` is a tuple of shapely polygons or multipolygons + in `lakes_crs` (the CLI reads them from `--lakes`; paths stop there). +- `Catchment`: `fine` (a shapely `Polygon` in the DEM's CRS, no interiors), + `crs` (text pyproj reads as the DEM's CRS), `nodes` (in-nodes), + `seed_nodes`, `fine_area` (m^2), `rings_dropped`, `holes_filled`, + `windows` (one entry per flood), and the last window's `mask` (in-nodes + non-zero) with its `meta` (`RasterMeta`). +- `CatchmentError` is a `ValueError`. + +THE WINDOW'S REFERENCE is one `_core.upstream` over the whole raster with the +same seed mask: the loop must end equal to it, node for node. + +THE FLAGS: the C++ suite pins `touches_edge` and `touches_nodata` only where +the design's words and any repair agree (see its header). The refusals here +are pinned by behaviour: a catchment cut by the data's edge, or by NoData, is +refused, whatever mechanism detects it. +""" + +from __future__ import annotations + +import asyncio +from typing import Any + +import numpy as np +import pytest +import shapely +from pydantic import ValidationError +from pyproj import Transformer +from shapely.geometry import MultiPolygon, Point, Polygon, box + +import tin_engine.mosaic as mosaic +from catchment_fixtures import ( + EPSG, + MARGIN_CELLS, + X0, + Y0, + D, + MemoryRepository, + bowl, + filled, + full_flood, + lake_box, + lat, + placed, + seeds_in, + tile_of, +) +from mosaic_fixtures import quadrants +from test_outline import square_area +from tin_engine.crs import parse_crs + + +@pytest.fixture(scope="module") +def api() -> Any: + import tin_engine.catchment as catchment + + return catchment + + +def repository_of(z: np.ndarray, nodata: float | None = None) -> MemoryRepository: + rows, cols = z.shape + return MemoryRepository( + quadrants(tile_of(z, nodata), row_cut=rows // 2, col_cut=cols // 2, overlap=1) + ) + + +def to_4326(x: float, y: float) -> tuple[float, float]: + lon, lat_ = Transformer.from_crs(EPSG, "EPSG:4326", always_xy=True).transform(x, y) + return float(lon), float(lat_) + + +def moved(geometry: Polygon, crs: str) -> Polygon: + t = Transformer.from_crs(EPSG, crs, always_xy=True) + return shapely.transform(geometry, lambda xy: np.column_stack(t.transform(xy[:, 0], xy[:, 1]))) + + +CENTRE = lat(100, 200) # the bowl's centre node, inside the lake + + +def request(api: Any, **kwargs: Any) -> Any: + kwargs.setdefault("seed", CENTRE) + kwargs.setdefault("seed_crs", EPSG) + if "lakes" not in kwargs: + kwargs["lakes"], kwargs["lakes_crs"] = (lake_box(),), EPSG + return api.CatchmentRequest(**kwargs) + + +def on_whole(result: Any, shape: tuple[int, int]) -> np.ndarray: + m = result.meta + assert np.asarray(result.mask).shape == (m.rows, m.cols) + assert m.delta_x == D and m.delta_y == D + return placed(np.asarray(result.mask) != 0, m.x_min, m.y_max, shape) + + +# --------------------------------------------------------------------------- +# The request and the error +# --------------------------------------------------------------------------- + + +def test_the_margin_is_2000_metres(api: Any) -> None: + assert api.WINDOW_MARGIN_M == 2000 + assert MARGIN_CELLS * D == api.WINDOW_MARGIN_M + + +def test_catchment_error_is_a_value_error(api: Any) -> None: + assert issubclass(api.CatchmentError, ValueError) + + +def test_the_request_is_frozen_and_seed_crs_defaults_to_wgs84(api: Any) -> None: + req = api.CatchmentRequest(seed=(8.5425, 61.3512)) + assert req.seed_crs == "EPSG:4326" + assert req.lakes is None + with pytest.raises(ValidationError): + req.seed = (0.0, 0.0) + + +# --------------------------------------------------------------------------- +# The seed +# --------------------------------------------------------------------------- + + +def test_a_lake_polygon_seeds_the_nodes_inside_it(api: Any) -> None: + z = bowl() + result = api.delineate(request(api), repository_of(z)) + seed = seeds_in(tile_of(z), lake_box()) + assert seed.sum() == 49 # the fixture's own premise: 7 x 7 nodes + assert result.seed_nodes == 49 + whole_mask = on_whole(result, z.shape) + assert np.all(whole_mask[seed == 1] == 1) + + +def test_a_lake_from_geojson_seeds_the_same_nodes(api: Any) -> None: + doc = {"type": "Polygon", "coordinates": [list(lake_box().exterior.coords)]} + lake = shapely.geometry.shape(doc) + z = bowl() + a = api.delineate(request(api, lakes=(lake,), lakes_crs=EPSG), repository_of(z)) + b = api.delineate(request(api), repository_of(z)) + assert a.seed_nodes == b.seed_nodes == 49 + assert np.array_equal(on_whole(a, z.shape), on_whole(b, z.shape)) + + +def test_a_multipolygon_contributes_the_part_containing_the_point(api: Any) -> None: + z = bowl() + far = box(*lat(20, 380), *lat(30, 370)) # 100 nodes, far from the bowl + lakes = (MultiPolygon([lake_box(), far]),) + result = api.delineate(request(api, lakes=lakes, lakes_crs=EPSG), repository_of(z)) + assert result.seed_nodes == 49 + + +def test_the_point_in_no_lake_is_refused_naming_the_point(api: Any) -> None: + x, y = lat(60, 60) + with pytest.raises(api.CatchmentError, match=f"{int(x)}"): + api.delineate(request(api, seed=(x, y)), repository_of(bowl())) + + +def test_the_point_in_two_lakes_is_refused(api: Any) -> None: + lakes = (lake_box(), lake_box().buffer(D)) + with pytest.raises(api.CatchmentError, match=f"{int(CENTRE[0])}"): + api.delineate(request(api, lakes=lakes, lakes_crs=EPSG), repository_of(bowl())) + + +def test_without_lakes_the_seed_is_the_nearest_node(api: Any) -> None: + z = bowl() + x, y = CENTRE + result = api.delineate( + api.CatchmentRequest(seed=(x + 30.0, y - 20.0), seed_crs=EPSG), repository_of(z) + ) + assert result.seed_nodes == 1 + seed = np.zeros(z.shape, dtype=np.uint8) + seed[200, 100] = 1 + assert np.array_equal(on_whole(result, z.shape), full_flood(tile_of(z), seed)) + + +def test_a_seed_in_wgs84_lands_on_the_same_node_as_in_the_dem_crs(api: Any) -> None: + z = bowl() + x, y = CENTRE[0] + 30.0, CENTRE[1] - 20.0 + utm = api.delineate(api.CatchmentRequest(seed=(x, y), seed_crs=EPSG), repository_of(z)) + wgs = api.delineate(api.CatchmentRequest(seed=to_4326(x, y)), repository_of(z)) + assert wgs.seed_nodes == 1 + assert np.array_equal(on_whole(wgs, z.shape), on_whole(utm, z.shape)) + + +def test_lakes_and_seed_in_other_crss_are_moved_into_the_dem_crs(api: Any) -> None: + z = bowl() + here = api.delineate(request(api), repository_of(z)) + there = api.delineate( + api.CatchmentRequest( + seed=to_4326(*CENTRE), lakes=(moved(lake_box(), "EPSG:3035"),), lakes_crs="EPSG:3035" + ), + repository_of(z), + ) + assert there.seed_nodes == 49 + assert np.array_equal(on_whole(there, z.shape), on_whole(here, z.shape)) + + +# --------------------------------------------------------------------------- +# The window +# --------------------------------------------------------------------------- + + +def test_a_catchment_larger_than_the_first_window_grows_and_equals_one_flood(api: Any) -> None: + z = bowl() + tile = tile_of(z) + reference = full_flood(tile, seeds_in(tile, lake_box())) + rows = np.nonzero(reference.any(axis=1))[0] + # The fixture's premise: the lake's box plus the margin does not hold it. + assert rows.min() < 200 - 3 - MARGIN_CELLS + + result = api.delineate(request(api), repository_of(z)) + assert len(result.windows) >= 2 + assert np.array_equal(on_whole(result, z.shape), reference) + assert result.nodes == int(reference.sum()) + assert parse_crs(result.crs) == parse_crs(EPSG) + + +def test_the_fine_outline_is_the_filled_outer_ring_round_the_seed(api: Any) -> None: + z = bowl() + tile = tile_of(z) + reference = full_flood(tile, seeds_in(tile, lake_box())) + result = api.delineate(request(api), repository_of(z)) + fine = result.fine + assert isinstance(fine, Polygon) + assert fine.is_valid + assert len(fine.interiors) == 0 + assert fine.exterior.is_ccw + assert fine.contains(Point(CENTRE)) + assert fine.area == pytest.approx(square_area(filled(reference)) * D * D, rel=1e-12) + assert result.fine_area == pytest.approx(fine.area, rel=1e-12) + r, c = np.nonzero(reference) + x, y = X0 + D * c, Y0 - D * r + assert shapely.contains_xy(fine, x, y).all() + + +@pytest.mark.parametrize( + ("side", "kwargs"), [("north", {"rc": 40}), ("east", {"cc": 190, "ax": 15.0})] +) +def test_a_catchment_cut_by_the_data_edge_is_refused_naming_the_side( + api: Any, side: str, kwargs: dict[str, Any] +) -> None: + rc, cc = kwargs.get("rc", 200), kwargs.get("cc", 100) + z = bowl(**kwargs) + with pytest.raises(api.CatchmentError, match=rf"(?i)\b{side}\b"): + api.delineate( + request(api, seed=lat(cc, rc), lakes=(lake_box(rc, cc),), lakes_crs=EPSG), + repository_of(z), + ) + + +def test_a_catchment_cut_by_nodata_is_refused(api: Any) -> None: + z = bowl() + z[150:153, 99:102] = np.nan # inside the basin, upstream of the lake + with pytest.raises(api.CatchmentError, match=r"(?i)no ?data"): + api.delineate(request(api), repository_of(z)) + + +def test_a_catchment_cut_by_the_sentinel_is_refused(api: Any) -> None: + z = bowl() + z[150:153, 99:102] = -32767.0 + with pytest.raises(api.CatchmentError, match=r"(?i)no ?data"): + api.delineate(request(api), repository_of(z, nodata=-32767.0)) + + +def test_the_memory_cap_refuses_before_the_flood(api: Any, monkeypatch: pytest.MonkeyPatch) -> None: + # The first window is the lake's box (7 cells) plus 2 x 20 cells: about 48 + # x 48 nodes. With physical memory 10 x that, 15a's cap (4 bytes a node + # against half of it) passes and the flood's (4 + 2 bytes a node) refuses. + first = 48 * 48 + fake = 10 * first + monkeypatch.setattr(mosaic, "physical_memory", lambda: fake) + monkeypatch.setattr(api, "physical_memory", lambda: fake, raising=False) + with pytest.raises(api.CatchmentError, match=r"(?i)memory"): + api.delineate(request(api), repository_of(bowl())) + + +# --------------------------------------------------------------------------- +# Holes filled, other pieces dropped +# --------------------------------------------------------------------------- + + +def plus_with_a_hole() -> tuple[np.ndarray, Polygon]: + """A 70 x 70 terrain falling from the centre (35, 35) to every edge, and a + lake whose nodes are the plus round the centre (not the centre) and one + node (35, 45) joined to it by a corridor holding no node. + + The plus is at 100, the centre at 200, the rest 50 - distance: the + centre's diagonal neighbours (48.6) are flooded before the plus, so they + reach the centre first and it is out, a hole in a catchment of five + nodes in two pieces.""" + r, c = np.indices((70, 70)).astype(np.float64) + z = (50.0 - np.hypot(r - 35, c - 35)).astype(np.float32) + for rr, cc in ((34, 35), (36, 35), (35, 34), (35, 36), (35, 45)): + z[rr, cc] = 100.0 + z[35, 35] = 200.0 + + def diamond(radius: float) -> Polygon: + return Polygon( + [lat(35 + radius, 35), lat(35, 35 - radius), lat(35 - radius, 35), lat(35, 35 + radius)] + ) + + annulus = diamond(1.4).difference(diamond(0.5)) + corridor = box(*lat(36, 35.4), *lat(45, 35.3)) + square = box(*lat(44.6, 35.4), *lat(45.4, 34.6)) + lake = shapely.union_all([annulus, corridor, square]) + assert isinstance(lake, Polygon) and len(lake.interiors) == 1 + return z, lake + + +def test_holes_are_filled_and_other_pieces_dropped_and_counted(api: Any) -> None: + z, lake = plus_with_a_hole() + tile = tile_of(z) + seed = seeds_in(tile, lake) + assert {(int(a), int(b)) for a, b in zip(*np.nonzero(seed), strict=True)} == { + (34, 35), (36, 35), (35, 34), (35, 36), (35, 45) + } # fmt: skip + result = api.delineate( + api.CatchmentRequest(seed=lat(35, 34), seed_crs=EPSG, lakes=(lake,), lakes_crs=EPSG), + MemoryRepository({"t.tif": tile}), + ) + assert result.nodes == 5 + assert result.seed_nodes == 5 + assert result.rings_dropped == 1 + assert result.holes_filled == 1 + fine = result.fine + assert len(fine.interiors) == 0 + assert fine.contains(Point(lat(35, 35))) # the hole, filled + assert not fine.contains(Point(lat(45, 35))) # the other piece, dropped + assert fine.area == pytest.approx(4.5 * D * D, rel=1e-12) + + +# --------------------------------------------------------------------------- +# Async callers +# --------------------------------------------------------------------------- + + +async def test_delineate_runs_in_a_worker_thread(api: Any) -> None: + z = bowl() + here = api.delineate(request(api), repository_of(z)) + there = await asyncio.to_thread(api.delineate, request(api), repository_of(z)) + assert np.array_equal(np.asarray(there.mask), np.asarray(here.mask)) + assert there.fine.equals_exact(here.fine, 0.0) diff --git a/tests/python/test_cli_catchment.py b/tests/python/test_cli_catchment.py new file mode 100644 index 00000000..2b31c738 --- /dev/null +++ b/tests/python/test_cli_catchment.py @@ -0,0 +1,283 @@ +"""`rasputin catchment`, and its output meshed by `rasputin mesh --domain` (increment 22, PR 1). + +`docs/increments/22-auto-catchment.md`, "The command", "The seed" and "The +red suites" (PR 1, `test_cli_catchment.py`), and the Bygdin numbers under +"What the data says about Bygdin". + +Pinned here, from the design: + +- `rasputin catchment --dem PATH --seed X Y [--seed-crs CRS] [--lakes PATH + [--lakes-layer NAME]] --out FILE.geojson`. `--seed-crs` defaults to + EPSG:4326 (`LON LAT`). `--lakes` is a `.gpkg` (with `--lakes-layer` when it + holds more than one features table) or `.geojson`/`.json`, in the file's + own CRS. +- The output is a `FeatureCollection` of one `Feature`, the polygon in the + DEM's CRS, with a `crs` member naming it; `domain.read_domain` reads it. + In PR 1 the polygon is the fine outline, so its area is the marching-squares + area of the catchment's nodes, holes filled; coordinates survive the round + trip (`repr` precision), so the area read back is that, to 1e-12. +- stderr has a line per window, and the fine outline's vertex count and area. + The wording is not pinned; these tests look for the words `window`, `fine`, + `vertices` and `km`. +- Every refusal is a non-zero exit and no file. + +The Bygdin test runs against Ola's data in ../rasputin_data and is skipped +when it is absent: the area must be within 2 % of NVE's 305.54 km^2 (delfelt +1187), as the design's acceptance asks. +""" + +from __future__ import annotations + +import json +import re +from pathlib import Path +from typing import Any + +import numpy as np +import pytest +from shapely.geometry import Point +from typer.testing import CliRunner + +from catchment_fixtures import EPSG, D, bowl, filled, full_flood, lake_box, lat, seeds_in, tile_of +from gpkg_fixtures import DTM10, OLA_NORWAY, Layer, Row, write_gpkg +from mosaic_fixtures import quadrants +from test_catchment import moved, to_4326 +from test_cli_mesh import plain +from test_cli_mesh_mosaic import write_tiles +from test_outline import square_area +from tin_engine.cli import app +from tin_engine.crs import parse_crs, reprojector +from tin_engine.domain import read_domain + +runner = CliRunner(env={"NO_COLOR": "1", "TERM": "dumb"}) + +SEED = lat(100, 200) # the bowl's centre node, inside the lake +NVE_BYGDIN_KM2 = 305.54 # NVE delfelt 1187, fetched 2026-09-29 (see the design) +BYGDIN_SEED = ("8.5425", "61.3512") + + +def invoke(*args: str) -> tuple[int, str]: + result = runner.invoke(app, list(args)) + return result.exit_code, plain(result.output) + + +def utm_seed(x: float = SEED[0], y: float = SEED[1]) -> tuple[str, ...]: + return ("--seed", repr(x), repr(y), "--seed-crs", EPSG) + + +def lake_geojson(path: Path, lake: Any = None, crs: str = EPSG) -> Path: + lake = lake_box() if lake is None else lake + doc = { + "type": "FeatureCollection", + "crs": {"type": "name", "properties": {"name": crs}}, + "features": [ + { + "type": "Feature", + "properties": {"name": "the lake"}, + "geometry": {"type": "Polygon", "coordinates": [list(lake.exterior.coords)]}, + } + ], + } + path.write_text(json.dumps(doc)) + return path + + +@pytest.fixture +def dem_dir(tmp_path: Path) -> Path: + write_tiles(tmp_path / "tiles", quadrants(tile_of(bowl()), row_cut=200, col_cut=100, overlap=1)) + return tmp_path / "tiles" + + +@pytest.fixture +def lakes(tmp_path: Path) -> Path: + return lake_geojson(tmp_path / "lakes.geojson") + + +def run(dem: Path, out: Path, *args: str) -> str: + code, output = invoke("catchment", "--dem", str(dem), *args, "--out", str(out)) + assert code == 0, output + assert out.is_file() + return output + + +def refused(tmp_path: Path, *args: str, out: str = "c.geojson") -> str: + target = tmp_path / out + code, output = invoke("catchment", *args, "--out", str(target)) + assert code != 0, output + # Refused by the command, not by Typer for want of one. + assert "No such command" not in output, output + assert not target.exists() + return output + + +def expected_area() -> float: + tile = tile_of(bowl()) + return square_area(filled(full_flood(tile, seeds_in(tile, lake_box())))) * D * D + + +# --------------------------------------------------------------------------- +# The round trip +# --------------------------------------------------------------------------- + + +def test_the_output_is_one_feature_in_the_dem_crs( + tmp_path: Path, dem_dir: Path, lakes: Path +) -> None: + out = tmp_path / "c.geojson" + run(dem_dir, out, *utm_seed(), "--lakes", str(lakes)) + doc = json.loads(out.read_text()) + assert doc["type"] == "FeatureCollection" + (feature,) = doc["features"] + assert feature["type"] == "Feature" + assert feature["geometry"]["type"] == "Polygon" + assert parse_crs(doc["crs"]["properties"]["name"]) == parse_crs(EPSG) + + +def test_read_domain_reads_it_back_with_the_fine_area( + tmp_path: Path, dem_dir: Path, lakes: Path +) -> None: + out = tmp_path / "c.geojson" + run(dem_dir, out, *utm_seed(), "--lakes", str(lakes)) + domain = read_domain(out) + assert parse_crs(domain.crs) == parse_crs(EPSG) + assert domain.polygon.is_valid + assert len(domain.polygon.interiors) == 0 + assert domain.polygon.area == pytest.approx(expected_area(), rel=1e-12) + + +def test_mesh_accepts_it_as_a_domain(tmp_path: Path, dem_dir: Path, lakes: Path) -> None: + out = tmp_path / "c.geojson" + run(dem_dir, out, *utm_seed(), "--lakes", str(lakes)) + vtk = tmp_path / "m.vtk" + code, output = invoke( + "mesh", "--dem", str(dem_dir), "--domain", str(out), "--tolerance", "5", "--out", str(vtk) + ) + assert code == 0, output + assert vtk.is_file() and vtk.stat().st_size > 0 + + +def test_stderr_reports_the_windows_and_the_fine_outline( + tmp_path: Path, dem_dir: Path, lakes: Path +) -> None: + output = run(dem_dir, tmp_path / "c.geojson", *utm_seed(), "--lakes", str(lakes)) + assert re.search(r"(?i)\bwindow", output), output + assert re.search(r"(?i)\bfine\b.*\bvertices\b|\bvertices\b.*\bfine\b", output), output + assert re.search(r"(?i)\bfine\b.*km|km.*\bfine\b", output), output + + +# --------------------------------------------------------------------------- +# The seed +# --------------------------------------------------------------------------- + + +def test_the_seed_crs_defaults_to_wgs84(tmp_path: Path, dem_dir: Path, lakes: Path) -> None: + utm, wgs = tmp_path / "utm.geojson", tmp_path / "wgs.geojson" + run(dem_dir, utm, *utm_seed(), "--lakes", str(lakes)) + lon, lat_ = to_4326(*SEED) + run(dem_dir, wgs, "--seed", repr(lon), repr(lat_), "--lakes", str(lakes)) + assert read_domain(wgs).polygon.equals_exact(read_domain(utm).polygon, 0.0) + + +def test_a_lake_file_in_another_crs(tmp_path: Path, dem_dir: Path, lakes: Path) -> None: + here, there = tmp_path / "here.geojson", tmp_path / "there.geojson" + run(dem_dir, here, *utm_seed(), "--lakes", str(lakes)) + other = lake_geojson(tmp_path / "l3035.geojson", moved(lake_box(), "EPSG:3035"), "EPSG:3035") + run(dem_dir, there, *utm_seed(), "--lakes", str(other)) + assert read_domain(there).polygon.equals_exact(read_domain(here).polygon, 0.0) + + +def gpkg_lakes(path: Path) -> Path: + """Two features tables in EPSG:3035: `lakes` holding the lake, `other` a + box far from it. With two tables, `--lakes-layer` is needed.""" + lake = Row(1, moved(lake_box(), "EPSG:3035"), {"Code_18": "512"}) + far = Row(1, moved(lake_box(380, 20), "EPSG:3035"), {"Code_18": "512"}) + layers = [Layer("lakes", 3035, [lake], rtree=False), Layer("other", 3035, [far], rtree=False)] + return write_gpkg(path, layers) + + +def test_a_geopackage_layer_seeds_like_the_geojson( + tmp_path: Path, dem_dir: Path, lakes: Path +) -> None: + here, there = tmp_path / "here.geojson", tmp_path / "there.geojson" + run(dem_dir, here, *utm_seed(), "--lakes", str(lakes)) + gpkg = gpkg_lakes(tmp_path / "lakes.gpkg") + run(dem_dir, there, *utm_seed(), "--lakes", str(gpkg), "--lakes-layer", "lakes") + assert read_domain(there).polygon.equals_exact(read_domain(here).polygon, 0.0) + + +def test_a_geopackage_of_two_tables_without_a_layer_is_refused( + tmp_path: Path, dem_dir: Path +) -> None: + gpkg = gpkg_lakes(tmp_path / "lakes.gpkg") + refused(tmp_path, "--dem", str(dem_dir), *utm_seed(), "--lakes", str(gpkg)) + + +def test_without_lakes_the_seed_is_a_pour_point(tmp_path: Path, dem_dir: Path) -> None: + out = tmp_path / "c.geojson" + run(dem_dir, out, *utm_seed(SEED[0] + 30.0, SEED[1] - 20.0)) + tile = tile_of(bowl()) + seed = np.zeros((400, 200), dtype=np.uint8) + seed[200, 100] = 1 + area = square_area(filled(full_flood(tile, seed))) * D * D + assert read_domain(out).polygon.area == pytest.approx(area, rel=1e-12) + + +# --------------------------------------------------------------------------- +# Refusals: a non-zero exit and no file +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize("out", ["c.vtk", "c.wkt", "c"]) +def test_an_output_that_is_not_geojson_is_refused( + tmp_path: Path, dem_dir: Path, lakes: Path, out: str +) -> None: + refused(tmp_path, "--dem", str(dem_dir), *utm_seed(), "--lakes", str(lakes), out=out) + + +def test_lakes_layer_without_lakes_is_refused(tmp_path: Path, dem_dir: Path) -> None: + output = refused(tmp_path, "--dem", str(dem_dir), *utm_seed(), "--lakes-layer", "lakes") + assert "--lakes" in output + + +def test_a_seed_in_no_lake_is_refused(tmp_path: Path, dem_dir: Path, lakes: Path) -> None: + refused(tmp_path, "--dem", str(dem_dir), *utm_seed(*lat(60, 60)), "--lakes", str(lakes)) + + +def test_a_catchment_cut_by_the_data_edge_writes_nothing(tmp_path: Path) -> None: + write_tiles( + tmp_path / "cut", quadrants(tile_of(bowl(rc=40)), row_cut=200, col_cut=100, overlap=1) + ) + lake = lake_geojson(tmp_path / "l.geojson", lake_box(40, 100)) + output = refused( + tmp_path, "--dem", str(tmp_path / "cut"), *utm_seed(*lat(100, 40)), "--lakes", str(lake) + ) + assert re.search(r"(?i)\bnorth\b", output), output + + +# --------------------------------------------------------------------------- +# Bygdin, on Ola's data +# --------------------------------------------------------------------------- + + +@pytest.mark.skipif( + not (DTM10.is_dir() and OLA_NORWAY.is_file()), reason="Ola's DTM10 and CORINE are not here" +) +def test_bygdin_is_within_two_percent_of_nve(tmp_path: Path) -> None: + out = tmp_path / "bygdin.geojson" + output = run( + DTM10, + out, + "--seed", + *BYGDIN_SEED, + "--lakes", + str(OLA_NORWAY), + "--lakes-layer", + "corine2018", + ) + domain = read_domain(out) + assert parse_crs(domain.crs) == parse_crs(EPSG) + km2 = domain.polygon.area / 1e6 + assert abs(km2 / NVE_BYGDIN_KM2 - 1.0) <= 0.02, (km2, output) + # The seed is mid-lake, so the outline contains it. + ((x, y),) = reprojector("EPSG:4326", EPSG)([[float(v) for v in BYGDIN_SEED]]) + assert domain.polygon.contains(Point(x, y)) diff --git a/tests/python/test_core_upstream.py b/tests/python/test_core_upstream.py new file mode 100644 index 00000000..37ede905 --- /dev/null +++ b/tests/python/test_core_upstream.py @@ -0,0 +1,122 @@ +"""`_core.upstream`, the flood's binding (increment 22, PR 1). + +`docs/increments/22-auto-catchment.md`, "Where it runs" and "Flow and +membership: one flood". The flood's behaviour is the C++ suite's +(`tests/cpp/unit/test_hydrology_upstream.cpp`); this file pins what crosses +the binding. + +Interface assumed (names chosen here, stated in the handback; they mirror the +C++ `UpstreamOutcome`): + +- `_core.upstream(view: RasterView, seed: NDArray[uint8]) -> UpstreamOutcome`, + `seed` of the raster's shape `(rows, cols)`, non-zero for a seed. Any other + shape is a `ValueError`. Releases the GIL. +- `UpstreamOutcome.mask`: a `uint8` array of shape `(rows, cols)`, 1 for an + in-node and 0 otherwise, owned by the outcome (it outlives the view). + `nodes_in`, `row_min`, `row_max`, `col_min`, `col_max` (inclusive), + `touches_edge`, `touches_nodata`. +""" + +from __future__ import annotations + +import gc +from typing import Any + +import numpy as np +import pytest + +from tin_engine._core import raster_view + +ROWS, COLS, RV, W = 21, 15, 10, 6 + + +def valley(dtype: type = np.float32) -> np.ndarray: + """The C++ suite's V-valley: falling west at 0.37 a column, sides rising at + 1 a row to ridges at rows 4 and 16, outer slopes at 2.9; column 0 a wall + at 1000 except the seed (10, 0) at -1.""" + r, c = np.indices((ROWS, COLS)) + off = np.abs(r - RV).astype(dtype) + side = np.where(off <= W, off, dtype(W) - dtype(2.9) * (off - dtype(W))) + z = (dtype(0.37) * c.astype(dtype) + side).astype(dtype) + z[:, 0] = 1000 + z[RV, 0] = -1 + z.flags.writeable = False + return z + + +def view_of(z: np.ndarray, nodata: float | None = None) -> Any: + return raster_view( + z, x_min=500_000.0, y_max=6_600_000.0, delta_x=10.0, delta_y=5.0, nodata=nodata + ) + + +def seed_at(*nodes: tuple[int, int]) -> np.ndarray: + seed = np.zeros((ROWS, COLS), dtype=np.uint8) + for r, c in nodes: + seed[r, c] = 1 + return seed + + +@pytest.fixture(scope="module") +def upstream() -> Any: + from tin_engine._core import upstream as the_upstream + + return the_upstream + + +def expected_valley() -> np.ndarray: + e = np.zeros((ROWS, COLS), dtype=np.uint8) + e[RV, 0] = 1 + e[RV - 5 : RV + 6, 1 : COLS - 1] = 1 + return e + + +@pytest.mark.parametrize("dtype", [np.float32, np.float64]) +def test_the_valley_through_the_binding(upstream: Any, dtype: type) -> None: + out = upstream(view_of(valley(dtype)), seed_at((RV, 0))) + assert isinstance(out.mask, np.ndarray) + assert out.mask.dtype == np.uint8 + assert out.mask.shape == (ROWS, COLS) + assert np.array_equal(out.mask, expected_valley()) + assert out.nodes_in == 11 * 13 + 1 + assert (out.row_min, out.row_max, out.col_min, out.col_max) == (5, 15, 0, 13) + assert out.touches_edge is True + assert out.touches_nodata is False + + +def test_the_mask_outlives_the_view(upstream: Any) -> None: + view = view_of(valley()) + out = upstream(view, seed_at((RV, 0))) + del view + gc.collect() + assert np.array_equal(out.mask, expected_valley()) + + +def test_a_bool_seed_is_accepted_as_well(upstream: Any) -> None: + out = upstream(view_of(valley()), seed_at((RV, 0)).astype(bool)) + assert np.array_equal(out.mask, expected_valley()) + + +@pytest.mark.parametrize( + "shape", [(ROWS, COLS - 1), (ROWS + 1, COLS), (COLS, ROWS), (ROWS * COLS,), (0, 0)] +) +def test_a_seed_of_another_shape_is_a_value_error(upstream: Any, shape: tuple[int, ...]) -> None: + with pytest.raises(ValueError, match=r"(?i)shape|size"): + upstream(view_of(valley()), np.zeros(shape, dtype=np.uint8)) + + +def test_the_sentinel_is_nodata_and_never_in(upstream: Any) -> None: + z = np.array(valley(), copy=True) + z[RV, 5] = -9999.0 + z.flags.writeable = False + out = upstream(view_of(z, nodata=-9999.0), seed_at((RV, 5), (RV, 6))) + assert out.mask[RV, 5] == 0 + assert out.mask[RV, 6] == 1 + assert out.touches_nodata is True + + +def test_no_seed_no_catchment(upstream: Any) -> None: + out = upstream(view_of(valley()), np.zeros((ROWS, COLS), dtype=np.uint8)) + assert out.nodes_in == 0 + assert not out.mask.any() + assert out.touches_edge is False and out.touches_nodata is False diff --git a/tests/python/test_outline.py b/tests/python/test_outline.py new file mode 100644 index 00000000..0901b2bf --- /dev/null +++ b/tests/python/test_outline.py @@ -0,0 +1,253 @@ +"""The fine outline: marching squares on a catchment's node mask (increment 22, PR 1). + +`docs/increments/22-auto-catchment.md`, "The fine outline" and "The red +suites" (PR 1, `test_outline.py`). No `_core`: the tracer is numpy only. + +Interface assumed (the design fixes the module, the function and the +guarantee; the representation below is chosen here and stated in the +handback): + +- `tin_engine.outline.trace(mask) -> list[numpy.ndarray]`. `mask` is a 2-D + array, bool or uint8, non-zero for an in-node. Each ring is an `(k, 2)` + float array of `(row, col)` in the mask's own indices (so a vertex on the + padding side of row 0 has row -0.5). Every vertex is the midpoint of a + lattice edge between an in-node and an out-node: one coordinate is an + integer and the other a half. A ring may or may not repeat its first vertex + at the end; this suite accepts both and counts vertices without the repeat. +- In the world frame (x = col, y = -row: x east, y north, as the tracer's + `x = x_min + col * dx`, `y = y_max - row * dy` gives), rings at even + nesting depth (outer) are counter-clockwise and rings at odd depth (holes) + clockwise. + +THE AREA ORACLE is independent of any tracing: marching squares with the +(8, 4) saddle rule gives each square of four nodes a fixed in-area by its +corners (0 in: 0; 1: 1/8; 2 adjacent: 1/2; 2 diagonal, joined: 3/4; 3: 7/8; +4: 1). The signed sum of the rings' areas (outer positive, holes negative) +must equal the sum over every square of the padded mask. +""" + +from __future__ import annotations + +import itertools +from collections.abc import Callable +from typing import Any + +import numpy as np +import numpy.typing as npt +import pytest +import shapely +from shapely.geometry import LinearRing, Point, Polygon + +from importscan import first_party_imports + +Ring = npt.NDArray[np.float64] +Trace = Callable[[Any], list[Any]] + + +@pytest.fixture(scope="module") +def trace() -> Trace: + from tin_engine.outline import trace as the_trace + + return the_trace + + +def open_ring(ring: npt.ArrayLike) -> Ring: + """The ring's vertices without a repeated closing vertex.""" + a = np.asarray(ring, dtype=np.float64) + assert a.ndim == 2 and a.shape[1] == 2, a.shape + if len(a) > 1 and np.array_equal(a[0], a[-1]): + a = a[:-1] + return a + + +def world(ring: Ring) -> list[tuple[float, float]]: + """(row, col) -> (x, y) = (col, -row).""" + return [(float(c), float(-r)) for r, c in ring] + + +def signed_area(ring: Ring) -> float: + """Shoelace in the world frame: positive counter-clockwise.""" + xy = np.asarray(world(ring)) + x, y = xy[:, 0], xy[:, 1] + return 0.5 * float(np.dot(x, np.roll(y, -1)) - np.dot(np.roll(x, -1), y)) + + +def square_area(mask: npt.ArrayLike) -> float: + """The oracle: the marching-squares in-area of `mask`, square by square.""" + m = np.pad(np.asarray(mask) != 0, 1).astype(int) + a, b, c, d = m[:-1, :-1], m[:-1, 1:], m[1:, :-1], m[1:, 1:] # tl, tr, bl, br + n = a + b + c + d + diagonal = (n == 2) & (a == d) + area = np.select( + [n == 1, (n == 2) & ~diagonal, diagonal, n == 3, n == 4], [1 / 8, 1 / 2, 3 / 4, 7 / 8, 1.0] + ) + return float(area.sum()) + + +def depth(ring: Ring, others: list[Ring]) -> int: + """How many other rings enclose this one (tested at its first vertex).""" + p = Point(world(ring)[0]) + return sum(Polygon(world(o)).contains(p) for o in others) + + +def check_guarantee(mask: npt.NDArray[np.uint8], rings: list[Ring]) -> None: + """The design's guarantee, every clause, on one mask.""" + rows, cols = mask.shape + for ring in rings: + assert len(ring) >= 4 + frac = np.abs(ring - np.round(ring)) + halves = np.isclose(frac, 0.5) + # One coordinate an integer, the other a half, on every vertex. + assert np.all(halves.sum(axis=1) == 1), ring + assert np.all(np.isclose(frac, 0.0) | halves) + assert ring[:, 0].min() >= -0.5 and ring[:, 0].max() <= rows - 0.5 + assert ring[:, 1].min() >= -0.5 and ring[:, 1].max() <= cols - 0.5 + lr = LinearRing(world(ring)) + assert lr.is_valid and lr.is_simple, shapely.validation.explain_validity(Polygon(lr)) + # No two rings share a point. + for a, b in itertools.combinations(rings, 2): + assert LinearRing(world(a)).disjoint(LinearRing(world(b))) + # Outer counter-clockwise, holes clockwise, by nesting depth. + for i, ring in enumerate(rings): + d = depth(ring, rings[:i] + rings[i + 1 :]) + assert (signed_area(ring) > 0) == (d % 2 == 0), (i, d, signed_area(ring)) + # Area: the signed sum is the square-by-square oracle. + assert sum(map(signed_area, rings)) == pytest.approx(square_area(mask), abs=1e-12) + # In-nodes inside an odd number of rings, out-nodes an even number, each at + # least a quarter of the cell's diagonal from every ring. + polygons = [Polygon(world(r)) for r in rings] + boundaries = [LinearRing(world(r)) for r in rings] + quarter = np.sqrt(2.0) / 4 + for r, c in itertools.product(range(rows), range(cols)): + p = Point(float(c), float(-r)) + inside = sum(poly.contains(p) for poly in polygons) + assert inside % 2 == (1 if mask[r, c] else 0), (r, c, inside) + for b in boundaries: + assert b.distance(p) >= quarter - 1e-12, (r, c) + + +def traced(trace: Trace, mask: npt.ArrayLike) -> list[Ring]: + m = np.asarray(mask, dtype=np.uint8) + rings = [open_ring(r) for r in trace(m)] + check_guarantee(m, rings) + return rings + + +def test_the_tracer_imports_no_core() -> None: + import tin_engine.outline as outline + + assert not any(name.startswith("tin_engine._core") for name in first_party_imports(outline)) + + +def test_an_empty_mask_has_no_ring(trace: Trace) -> None: + assert list(trace(np.zeros((4, 5), dtype=np.uint8))) == [] + + +def test_a_single_node_is_one_diamond_of_half_a_cell(trace: Trace) -> None: + mask = np.zeros((5, 6), dtype=np.uint8) + mask[2, 3] = 1 + (ring,) = traced(trace, mask) + assert len(ring) == 4 + assert {tuple(v) for v in ring.tolist()} == {(1.5, 3.0), (2.0, 3.5), (2.5, 3.0), (2.0, 2.5)} + assert signed_area(ring) == pytest.approx(0.5) + + +def test_a_two_by_two_block_is_an_octagon_of_three_and_a_half(trace: Trace) -> None: + mask = np.zeros((5, 5), dtype=np.uint8) + mask[1:3, 2:4] = 1 + (ring,) = traced(trace, mask) + assert len(ring) == 8 + assert signed_area(ring) == pytest.approx(3.5) + + +def test_an_l_is_one_ring(trace: Trace) -> None: + mask = np.zeros((6, 6), dtype=np.uint8) + mask[1:5, 1] = 1 + mask[4, 1:4] = 1 + (ring,) = traced(trace, mask) + assert signed_area(ring) > 0 + + +def test_a_diagonal_pair_is_one_ring_by_the_saddle_rule(trace: Trace) -> None: + mask = np.zeros((4, 4), dtype=np.uint8) + mask[1, 1] = mask[2, 2] = 1 + (ring,) = traced(trace, mask) + assert signed_area(ring) == pytest.approx(1 / 8 * 6 + 3 / 4) + + +def test_the_anti_diagonal_pair_is_one_ring_too(trace: Trace) -> None: + mask = np.zeros((4, 4), dtype=np.uint8) + mask[1, 2] = mask[2, 1] = 1 + assert len(traced(trace, mask)) == 1 + + +def test_a_pinch_is_one_simple_ring(trace: Trace) -> None: + # Two 2 x 2 blocks meeting at one diagonal. + mask = np.zeros((6, 6), dtype=np.uint8) + mask[1:3, 1:3] = 1 + mask[3:5, 3:5] = 1 + assert len(traced(trace, mask)) == 1 + + +def test_a_ring_of_nodes_is_an_outer_ring_and_a_clockwise_hole(trace: Trace) -> None: + mask = np.zeros((5, 5), dtype=np.uint8) + mask[1:4, 1:4] = 1 + mask[2, 2] = 0 + rings = traced(trace, mask) + assert len(rings) == 2 + outer, hole = sorted(rings, key=signed_area, reverse=True) + assert signed_area(outer) > 0 + assert signed_area(hole) == pytest.approx(-0.5) # the diamond around the out-node + + +def test_a_diagonal_ring_holds_a_hole_the_out_node_cannot_leave(trace: Trace) -> None: + # A plus of four in-nodes round an out-node: the four are 8-connected, so + # the centre is enclosed (the (8, 4) pairing), a clockwise hole of half a cell. + mask = np.zeros((5, 5), dtype=np.uint8) + for r, c in ((1, 2), (2, 1), (2, 3), (3, 2)): + mask[r, c] = 1 + rings = traced(trace, mask) + assert len(rings) == 2 + assert sorted(map(signed_area, rings)) == pytest.approx([-0.5, 4.5]) + + +def test_disjoint_pieces_are_disjoint_rings(trace: Trace) -> None: + mask = np.zeros((7, 9), dtype=np.uint8) + mask[1:3, 1:3] = 1 + mask[4, 6] = 1 + mask[1, 7] = 1 + assert len(traced(trace, mask)) == 3 + + +@pytest.mark.parametrize("dtype", [np.uint8, np.bool_]) +def test_a_mask_filling_the_window_closes_through_the_padding(trace: Trace, dtype: Any) -> None: + mask = np.ones((3, 4), dtype=dtype) + rings = [open_ring(r) for r in trace(mask)] + check_guarantee(mask.astype(np.uint8), rings) + (ring,) = rings + assert ring[:, 0].min() == -0.5 and ring[:, 0].max() == 2.5 + assert ring[:, 1].min() == -0.5 and ring[:, 1].max() == 3.5 + + +def test_a_single_row_and_a_single_column(trace: Trace) -> None: + traced(trace, np.ones((1, 5))) + traced(trace, np.ones((5, 1))) + traced(trace, np.ones((1, 1))) + + +@pytest.mark.parametrize("seed", range(40)) +def test_random_masks_meet_the_guarantee(trace: Trace, seed: int) -> None: + rng = np.random.default_rng(22_100 + seed) + rows, cols = (int(v) for v in rng.integers(1, 13, size=2)) + density = [0.2, 0.5, 0.8][seed % 3] + mask = (rng.random((rows, cols)) < density).astype(np.uint8) + traced(trace, mask) + + +def test_a_checkerboard_is_one_ring_with_holes(trace: Trace) -> None: + # Every in-node touches its diagonal neighbours: one 8-connected piece. + r, c = np.indices((7, 7)) + mask = ((r + c) % 2 == 0).astype(np.uint8) + rings = traced(trace, mask) + outer = [ring for ring in rings if signed_area(ring) > 0] + assert len(outer) == 1 From 06776de4e78929ba59d18eaeff878f20be50d2a9 Mon Sep 17 00:00:00 2001 From: Ola Skavhaug Date: Tue, 29 Sep 2026 01:33:31 +0200 Subject: [PATCH 07/14] 22 design: flags, loop termination, closed lakes, lake reading Co-Authored-By: Claude Opus 5.5 --- docs/increments/22-auto-catchment.md | 159 +++++++++++++++++++++++---- 1 file changed, 135 insertions(+), 24 deletions(-) diff --git a/docs/increments/22-auto-catchment.md b/docs/increments/22-auto-catchment.md index 5228c918..f18792ce 100644 --- a/docs/increments/22-auto-catchment.md +++ b/docs/increments/22-auto-catchment.md @@ -215,7 +215,7 @@ cli.py catchment v catchment.py delineate(request, repository) -> Catchment [no paths] 1. seed point -> DEM CRS (crs.reprojector); lake polygon containing it - -> DEM CRS (seed.py, shapely/pyproj) + -> DEM CRS (lakes read by cli.py via feature_input.read_lakes) 2. window = seed bounds grown by the margin 3. loop: plan_mosaic + assemble (15a) -> DemTile, read-only seed mask: DEM nodes inside the lake (shapely.contains_xy) @@ -281,12 +281,36 @@ Why this and not the alternatives: x 168680 and the dam at x 168087, so a sliver of river below the dam may be included. The acceptance run measures it (nodes in ours and not in NVE's). -The source is any polygon layer: a GeoPackage table (`--lakes-layer`, read -with 16b's `io/geopackage.py`, box query on the seed point, in the layer's -own CRS) or a GeoJSON `FeatureCollection`. A multipolygon contributes the -part containing the point. No class filter: the polygon under the point is -the one meant. Refusals: the point in no polygon, or in two (overlapping -input); both name the point in the source's CRS. +The source is any polygon layer: a GeoPackage table (`--lakes-layer`) or a +GeoJSON `FeatureCollection`. A multipolygon contributes the part containing +the point. No class filter: the polygon under the point is the one meant. +Refusals: the point in no polygon, or in two (overlapping input); both name +the point in the source's CRS. + +**Who reads the lakes** (revised 2026-09-29, after the PR 1 red step). +`feature_input.py`, the module that already reads 16b's sources, gains +`read_lakes(path, layer, point, point_crs) -> (tuple of shapely geometries, +crs text)`. It reuses 16b's reading by splitting the file-reading half of +`_Tally.source` into a public `read_source(path, layer, attribute, box_for)`, +which `_Tally.source` then calls, so 16b's behaviour is unchanged: + +- `.gpkg`: `io/repository.open_geopackage`, `io/geopackage.layer_info` + (`--lakes-layer`, or the only features table), then + `io/geopackage.query_features` with the box `box_for(layer CRS)`, here the + seed point moved into the layer's CRS (a point box; the R-tree widening + makes it a superset). The attribute column is the layer's primary key, + since lakes need no class. +- `.geojson` / `.json`: `json` and `shapely.geometry.shape` per feature, the + CRS from the `crs` member or WGS 84, as 16b does. +- `.gml`: comes free with `read_source`; not advertised. + +`read_lakes` returns every polygon or multipolygon the query yields, in the +source's CRS; lines and points are skipped. Which one contains the point, and +the refusals above, are `catchment.py`'s: `CatchmentRequest` carries +`lakes` (a tuple of shapely geometries) and `lakes_crs`, and never a path. +`cli.py` calls `read_lakes` and builds the request. This matches the red +suite as committed: `test_catchment.py` passes geometries, and +`test_cli_catchment.py` passes files through the CLI. **Without `--lakes`**, the seed is the one DEM node nearest the point, a pour point with no snapping. It is there because it costs five lines and gives the @@ -318,10 +342,36 @@ per node (unreached, out, in), which is also the returned mask, plus the queue is one more byte per node, owned by numpy. Returned: the mask (a numpy `uint8` array the outcome owns), the number of -nodes in, the bounding rows and columns of the in-nodes, and two flags: -`touches_edge` (an in-node on the window's edge) and `touches_nodata` (an -in-node with a NoData neighbour). A seed mask whose shape is not the raster's -is a `ValueError` in the binding. +nodes in, the bounding rows and columns of the in-nodes, and two flags. A seed +mask whose shape is not the raster's is a `ValueError` in the binding. + +**The flags** (revised 2026-09-29, after the PR 1 red step; the first wording, +"an in-node on the window's edge", could only ever fire for a seed, because +every edge node and every node beside NoData is an outlet and an outlet is in +only if it is a seed). The flags say where the catchment may continue beyond +what the window knows. An outlet's own drainage is unknown: the window assumed +it leaves, and in the full DEM it may instead run into the catchment. So: + +- `touches_edge` is true when some in-node is an edge outlet or is an + 8-neighbour of one. On a grid that is exactly: an in-node in the first or + last two rows or columns (row <= 1, row >= rows - 2, the same for columns). +- `touches_nodata` is true when some in-node is an outlet beside NoData or is + an 8-neighbour of one, which puts it within two nodes of NoData. + +An outlet can be both kinds; then both flags may be set. No comparison of +heights is made: an edge node lower than its in-neighbour might still, in the +full DEM, fill and spill back, so any contact counts. Conservative by design: +a flag never misses a truncated catchment, and a catchment that merely comes +within one node of an edge it does not cross is flagged too. The window loop +then grows past it, so the price is a growth step, not a wrong answer. + +Checked against the red suite as committed (1e1b3bb): the C++ +`check_invariants` asserts the flag when an in-node is on the edge or beside +NoData, and its absence when no in-node is within one node of the edge or two +of NoData; this rule sets the flag exactly on the first and never on the +second. The named cases (a seed on the edge, the V-valley, the flat lake far +from every edge, no seed, the degenerate rasters, the binding's seed beside +NoData) agree. **No test changes are needed.** ### The window @@ -332,22 +382,82 @@ re-plan pattern as 15b's `_domain_plan`: the margin, `WINDOW_MARGIN_M = 2000` metres. Default (main session / @architect, 2026-09-29), for Ola to confirm; no flag tonight. 2. Plan and assemble that box (15a), flood it. -3. If the in-nodes' bounds grown by the margin fit inside the window, stop. - Otherwise the new box is the old box joined with the in-nodes' bounds grown - by twice the margin, the margin doubles, and the loop repeats from 2. The - box only grows and the margin doubles, so the total work is within about - twice the last flood's. -4. Refuse, with the side named, when the catchment reaches a side the data - cannot extend (the planned window is smaller than the box asked for on - that side), or touches NoData: the catchment is truncated, and a truncated - catchment is a wrong one. Default (main session / @architect, - 2026-09-29), for Ola to confirm; alternative: write it with a warning - under an `--allow-truncated` flag. +3. and 4., revised 2026-09-29 after the PR 1 red step (the first wording + could loop for ever near the data's edge, where "the bounds plus the + margin fit" stays false while the window cannot grow). All comparisons + are in node indices on the plan's lattice, and all are closed: a box that + ends exactly on the window's last node line fits. + + Let E be the data's node rectangle on the chosen lattice (the union of its + tiles' extents, which `plan_mosaic` clamps every window to), and W the + window just flooded, already inside E. + + a. If `touches_nodata`, refuse: the catchment is truncated by missing + data, which no growth can fix. + b. The need N is the in-nodes' bounds grown by the margin, clamped to E. + A side of W *can grow* when N reaches past W on that side; since N is + clamped to E, that also means W is not yet at E there. + c. If some side can grow: the next window is W joined with N, the margin + doubles, and the loop repeats from 2. + d. Otherwise decide by the flag. If `touches_edge`, refuse, naming the + sides (from the bounds: an in-node in the first or last two rows or + columns), each of which is then at E: the catchment is cut by the + data's edge. If not, accept. + + **It terminates.** Step c runs only when N reaches past W on some side, + and the next window contains N, so it has at least one more row or column + than W; every window lies inside E, which is finite. So step c runs at + most (rows of E + columns of E) times, and in practice a handful, since + the margin doubles. The memory cap (step 5) may refuse earlier. Every + exit is an accept or a refusal from a or d. + + **Accepting means**: the catchment comes no closer than two nodes to any + edge of the window, and either the window already holds the in-nodes' + bounds plus the margin, or it is at the data's edge on the sides where it + does not. The margin is a heuristic against the known limit below; the + flag is the rule. + + Refusing a truncated catchment stays the default (main session / + @architect, 2026-09-29), for Ola to confirm; alternative: write it with a + warning under an `--allow-truncated` flag. 5. Memory: before each flood, refuse if nodes x (itemsize + 2) exceeds half the physical memory (`mosaic.physical_memory`, dtype-aware as 15a's cap is), naming the window's size. 15a's own cap on the array still applies inside `plan_mosaic`. +### A lake in a closed depression + +Added 2026-09-29, after the PR 1 red step. **A lake seed need not have an +outlet, and nothing is refused.** Default (main session / @architect, +2026-09-29), for Ola to confirm. + +Within a window every depression drains: Priority-Flood fills it to its +lowest rim and routes it out over that rim, so every lake has an outlet as +far as the flood is concerned, a real one or the spill point of its +depression. The catchment is still "every node whose flooding chain passes +through the lake". For the nodes of the depression itself that is what the +design wants when the lake polygon covers the depression, as it does for a +lake whose DEM surface is its water line: all those nodes are seeds. + +Where it falls short: when the filled depression is larger than the lake +polygon (a lake drawn smaller than its DEM basin, or a dry pan beside it +below the spill level), the nodes of the depression outside the polygon are +flooded first in, first out from the spill point, so by distance in steps. +Those reached through lake nodes are in; those reached round the lake from +the spill point are out, although water there would run into the lake. The +error is bounded by the part of the depression outside the polygon and +reached before the lake; on Bygdin, where CORINE's polygon covers the DEM's +water surface, the acceptance's node comparison with NVE's polygon is where +it would show. + +The alternative, for later: seed the whole depression. A first flood +computes the filled levels (one value per node, 4 bytes at float32, exact +because every level is some node's z); every node 8-connected to a seed +through nodes whose filled level is above their own z joins the seed; a +second flood labels as now. Twice the flood's time and 4 more bytes per +node, and a behaviour the red suite would need to pin (a seed inside a pit). +Not tonight. + **Known limit.** Not touching the window's edge does not prove the catchment complete. A closed depression that straddles the window's edge drains out through the edge in the window, but in the full DEM it may fill and spill into @@ -536,15 +646,16 @@ request and the result are frozen; the CLI is the only place with paths. | `src_python/tin_engine/outline.py` | the tracer | 70 | | `src_python/tin_engine/catchment.py` | request, seed, window loop, result | 170 | | `src_python/tin_engine/dem_input.py` | repository helper split out | 10 | +| `src_python/tin_engine/feature_input.py` | `read_source` split out of `_Tally.source`, `read_lakes` | 30 | | `src_python/tin_engine/cli.py` | `catchment` command, report, writer | 100 | | `project_structure.md` | the two C++ modules and two Python modules | docs | -About 830 lines, over the 700 ceiling (CLAUDE.md §2), so two PRs on this +About 860 lines, over the 700 ceiling (CLAUDE.md §2), so two PRs on this branch. ### The PR split -- **PR 1, the fine catchment** (about 480 lines): `upstream.hpp` and its +- **PR 1, the fine catchment** (about 510 lines): `upstream.hpp` and its binding, `outline.py`, `catchment.py`, the repository helper, and the `catchment` command writing the fine outline. It already answers "what is Bygdin's catchment" and can be compared against NVE. Red, green, review. From 58f69041b9db89cdcbd7c4c70e599152b84bcabd Mon Sep 17 00:00:00 2001 From: Ola Skavhaug Date: Tue, 29 Sep 2026 01:54:26 +0200 Subject: [PATCH 08/14] green: 22 PR1 the fine catchment (upstream flood, outline, window, rasputin catchment) Makes the PR 1 red suite (1e1b3bb) pass, per docs/increments/22-auto-catchment.md as amended in 06776de. No test file changed except removing the red-step if(EXISTS hydrology/upstream.hpp) guard in tests/cpp/CMakeLists.txt. - include/terrain/hydrology/upstream.hpp: Priority-Flood labelling of one catchment, keys (level, push counter); flags as "The flags" words them. - bindings/core.cpp, _core.pyi: _core.upstream and UpstreamOutcome (read-only (rows, cols) uint8 mask anchored to the outcome), GIL released. - outline.py: marching squares with the (8, 4) saddle rule, boundary squares only. - catchment.py: CatchmentRequest, delineate, the window loop, the fine outline. - feature_input.py: read_source split out of _Tally.source; read_lakes. - dem_input.py: repository_for, shared by open_dem and catchment. - cli.py: rasputin catchment, writing the fine outline as GeoJSON. Design deviation, for review: step c of "The window" doubles the margin and tests growth against the doubled margin, so every step grows and the loop only stops at the data's edge (Bygdin: 8 min, then refused by 15a's coverage check). Here growth is tested with the base 2000 m margin and the step takes the doubled one. Bygdin: 3 windows, 304.9096 km2 against NVE's 305.54 (-0.21 %), 6.3 s end to end. Co-Authored-By: Claude Opus 5.5 --- bindings/core.cpp | 59 +++++ include/terrain/hydrology/upstream.hpp | 134 ++++++++++ src_python/tin_engine/_core.pyi | 27 ++ src_python/tin_engine/catchment.py | 327 +++++++++++++++++++++++++ src_python/tin_engine/cli.py | 124 +++++++++- src_python/tin_engine/dem_input.py | 19 +- src_python/tin_engine/feature_input.py | 156 ++++++++---- src_python/tin_engine/outline.py | 83 +++++++ tests/cpp/CMakeLists.txt | 18 +- 9 files changed, 874 insertions(+), 73 deletions(-) create mode 100644 include/terrain/hydrology/upstream.hpp create mode 100644 src_python/tin_engine/catchment.py create mode 100644 src_python/tin_engine/outline.py diff --git a/bindings/core.cpp b/bindings/core.cpp index 37fb5742..8e5b51d9 100644 --- a/bindings/core.cpp +++ b/bindings/core.cpp @@ -12,6 +12,7 @@ #include #include #include +#include #include #include #include @@ -291,6 +292,13 @@ template return BoundRasterView{array, terrain::raster::RasterView{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_(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 auto cols = static_cast(b.cols); + return readonly_view(self, b.outcome.mask.data(), + {static_cast(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; + 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(s.shape(0)) != g.rows() + || static_cast(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 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"); } diff --git a/include/terrain/hydrology/upstream.hpp b/include/terrain/hydrology/upstream.hpp new file mode 100644 index 00000000..144b8ab0 --- /dev/null +++ b/include/terrain/hydrology/upstream.hpp @@ -0,0 +1,134 @@ +#pragma once + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace terrain::hydrology { + +// The catchment of a seed set: every DEM node whose drainage path over the +// filled surface passes through a seed (docs/increments/22-auto-catchment.md, +// "Flow and membership: one flood"; Barnes, Lehman and Mulla 2014a, +// Priority-Flood, labelling one catchment). +struct UpstreamOutcome { + std::vector mask; // row-major, 1 in, 0 otherwise (out, NoData) + std::size_t nodes_in{}; + std::size_t row_min{}, row_max{}, col_min{}, col_max{}; // inclusive; nodes_in > 0 + // Where the catchment may continue beyond what the window knows: an + // in-node that is an edge outlet or an 8-neighbour of one, and likewise + // for an outlet beside NoData ("The flags"). + bool touches_edge{}; + bool touches_nodata{}; +}; + +namespace detail { + +enum : std::uint8_t { kUnreached = 0, kIn = 1, kOut = 2, kNoData = 3 }; + +// Calls f(j) for each 8-neighbour j of node i, in a fixed order: the order is +// part of the result, since first in, first out breaks ties. +template +void each_neighbour(std::size_t i, std::size_t rows, std::size_t cols, F&& f) { + const std::size_t r = i / cols, c = i % cols; + const std::size_t r0 = r > 0 ? r - 1 : 0, r1 = r + 1 < rows ? r + 1 : r; + const std::size_t c0 = c > 0 ? c - 1 : 0, c1 = c + 1 < cols ? c + 1 : c; + for (std::size_t rr = r0; rr <= r1; ++rr) + for (std::size_t cc = c0; cc <= c1; ++cc) + if (rr != r || cc != c) + f(rr * cols + cc); +} + +} // namespace detail + +// `seed` is row-major, one byte per node, non-zero for a seed; a NoData seed +// is ignored. Keys are (level, push counter): equal levels pop first in, +// first out, so the result depends on nothing but the input. +template +[[nodiscard]] UpstreamOutcome upstream(const R& z, std::span seed) { + using namespace detail; + const raster::RasterGeometry& g = z.geometry(); + const std::size_t rows = g.rows(), cols = g.cols(), n = g.size(); + if (seed.size() != n) + throw std::invalid_argument("upstream: the seed mask's size is not the raster's"); + + UpstreamOutcome out; + std::vector& state = out.mask; + state.assign(n, kUnreached); + for (std::size_t r = 0; r < rows; ++r) + for (std::size_t c = 0; c < cols; ++c) + if (z.is_nodata(raster::CellIndex{r, c})) + state[r * cols + c] = kNoData; + + using Entry = std::tuple; + std::priority_queue, std::greater<>> queue; + std::uint64_t counter = 0; + const auto level_of = [&](std::size_t i) { + return static_cast(z.value_at(raster::CellIndex{i / cols, i % cols})); + }; + + // Outlets: every valid node on the window's edge or beside NoData. + std::vector beside_nodata; + for (std::size_t i = 0; i < n; ++i) { + if (state[i] == kNoData) + continue; + const std::size_t r = i / cols, c = i % cols; + bool nodata_near = false; + each_neighbour(i, rows, cols, [&](std::size_t j) { nodata_near |= state[j] == kNoData; }); + if (nodata_near) + beside_nodata.push_back(i); + if (nodata_near || r == 0 || c == 0 || r + 1 == rows || c + 1 == cols) { + state[i] = seed[i] != 0 ? kIn : kOut; + queue.emplace(level_of(i), counter++, i); + } + } + + while (!queue.empty()) { + const auto [level, order, i] = queue.top(); + queue.pop(); + const std::uint8_t label = state[i]; + each_neighbour(i, rows, cols, [&](std::size_t j) { + if (state[j] != kUnreached) + return; + state[j] = seed[j] != 0 ? kIn : label; + const double zj = level_of(j); + queue.emplace(zj > level ? zj : level, counter++, j); + }); + } + + out.row_min = rows; + out.col_min = cols; + for (std::size_t i = 0; i < n; ++i) { + if (state[i] != kIn) + continue; + const std::size_t r = i / cols, c = i % cols; + ++out.nodes_in; + out.row_min = r < out.row_min ? r : out.row_min; + out.row_max = r > out.row_max ? r : out.row_max; + out.col_min = c < out.col_min ? c : out.col_min; + out.col_max = c > out.col_max ? c : out.col_max; + } + if (out.nodes_in == 0) { + out.row_min = out.col_min = 0; + } else { + out.touches_edge = out.row_min <= 1 || out.col_min <= 1 || out.row_max + 2 >= rows + || out.col_max + 2 >= cols; + } + for (const std::size_t o : beside_nodata) { + bool near = state[o] == kIn; + each_neighbour(o, rows, cols, [&](std::size_t j) { near |= state[j] == kIn; }); + out.touches_nodata = out.touches_nodata || near; + } + for (auto& s : state) + s = s == kIn ? 1 : 0; + return out; +} + +} // namespace terrain::hydrology diff --git a/src_python/tin_engine/_core.pyi b/src_python/tin_engine/_core.pyi index 922758fb..d2663579 100644 --- a/src_python/tin_engine/_core.pyi +++ b/src_python/tin_engine/_core.pyi @@ -444,3 +444,30 @@ def refine( depend on ``threads``. ``min_angle_deg`` > 0 first improves the start mesh's angles with DEM nodes; 0 is off. ``constraint_feet`` inserts a node's foot on a nearby constraint segment instead of the node.""" + +@final +class UpstreamOutcome: + """What :func:`upstream` returned. Bounds are inclusive and meaningful + when ``nodes_in`` > 0.""" + + @property + def mask(self) -> npt.NDArray[np.uint8]: + """Read-only ``(rows, cols)``: 1 for a node in the catchment, 0 otherwise.""" + @property + def nodes_in(self) -> int: ... + @property + def row_min(self) -> int: ... + @property + def row_max(self) -> int: ... + @property + def col_min(self) -> int: ... + @property + def col_max(self) -> int: ... + @property + def touches_edge(self) -> bool: ... + @property + def touches_nodata(self) -> bool: ... + +def upstream(view: RasterView, seed: npt.ArrayLike) -> UpstreamOutcome: + """Every node draining into a seed of the ``(rows, cols)`` mask + (Priority-Flood); another shape is a ``ValueError``. Releases the GIL.""" diff --git a/src_python/tin_engine/catchment.py b/src_python/tin_engine/catchment.py new file mode 100644 index 00000000..deada36d --- /dev/null +++ b/src_python/tin_engine/catchment.py @@ -0,0 +1,327 @@ +"""The catchment of a lake, from the DEM (increment 22, PR 1: the fine outline). + +`docs/increments/22-auto-catchment.md`, "Data flow", "The seed", "The +window" and "The fine outline". :func:`delineate` moves the seed point and the +lake into the DEM's CRS, plans and assembles a window round them (15a), floods +it (`_core.upstream`), grows the window until the catchment lies clearly +inside it or is cut by the data's edge or NoData, and traces the outline of +the in-nodes, holes filled and other pieces dropped. + +No paths: the DEM comes through a `DemRepository`, and the lakes arrive as +shapely geometries the CLI read. Blocking (the flood releases the GIL); an +async caller runs :func:`delineate` in `asyncio.to_thread`. +""" + +from __future__ import annotations + +import math +import time +from dataclasses import dataclass +from typing import Any, Self + +import numpy as np +import numpy.typing as npt +import shapely +from pydantic import BaseModel, ConfigDict, model_validator +from shapely.geometry import Point, Polygon + +from tin_engine._core import UpstreamOutcome, upstream +from tin_engine.crs import crs_label, reprojector +from tin_engine.io.models import RasterMeta +from tin_engine.io.repository import DemRepository +from tin_engine.mosaic import ( + Bounds, + MosaicError, + MosaicPlan, + assemble, + physical_memory, + plan_mosaic, +) +from tin_engine.outline import trace +from tin_engine.raster import to_core + +#: Metres round the seed's bounds, and round the catchment's when the window +#: grows; doubled at each growth step. +WINDOW_MARGIN_M = 2000 +_SIDES = ("north", "west", "south", "east") + + +class CatchmentError(ValueError): + """A catchment that cannot be delineated, in words for the person asking.""" + + +class CatchmentRequest(BaseModel): + """The seed point in `seed_crs`, and the lakes (shapely polygons or + multipolygons in `lakes_crs`) of which the one containing it is the seed. + Without lakes the seed is the DEM node nearest the point.""" + + model_config = ConfigDict(frozen=True, arbitrary_types_allowed=True) + + seed: tuple[float, float] + seed_crs: str = "EPSG:4326" + lakes: tuple[Any, ...] | None = None + lakes_crs: str | None = None + + @model_validator(mode="after") + def _lakes_have_a_crs(self) -> Self: + if self.lakes is not None and self.lakes_crs is None: + raise ValueError("lakes need their CRS") + return self + + +@dataclass(frozen=True, slots=True) +class Window: + """One flood: the window's node extent, its shape, the flood's seconds, + and the sides it grew on afterwards (none: the catchment is contained).""" + + bounds: Bounds + rows: int + cols: int + seconds: float + grown: tuple[str, ...] + + +@dataclass(frozen=True, slots=True) +class Catchment: + """The fine outline (holes filled) in the DEM's CRS `crs`, the seed point + there, the lake's area (None without one), node counts, the outline's + area and what was left out of it, every window, and the last window's + mask of in-nodes with its `meta`.""" + + fine: Polygon + crs: str + seed: tuple[float, float] + lake_area: float | None + nodes: int + seed_nodes: int + fine_area: float + rings_dropped: int + dropped_nodes: int + holes_filled: int + holes_area: float + windows: tuple[Window, ...] + mask: npt.NDArray[np.uint8] + meta: RasterMeta + trace_seconds: float + + +def delineate(request: CatchmentRequest, repository: DemRepository) -> Catchment: + """The catchment of the request's seed over the repository's DEM, or a + :class:`CatchmentError` (truncated by the data's edge or NoData, no lake + or two under the point, over the memory cap).""" + footprints = repository.footprints() + epsgs = sorted({f.meta.epsg for f in footprints}) + if len(epsgs) != 1: + raise CatchmentError(f"the tiles are in {len(epsgs)} CRSs, EPSG:{epsgs}; need one") + dem_crs = f"EPSG:{epsgs[0]}" + ((x, y),) = reprojector(request.seed_crs, dem_crs)([request.seed]) + if not (math.isfinite(x) and math.isfinite(y)): + raise CatchmentError(f"the seed {request.seed} has no image in {dem_crs}") + lake = _lake(request, dem_crs) + x0, y0, x1, y1 = lake.bounds if lake is not None else (x, y, x, y) + margin = float(WINDOW_MARGIN_M) + plan = _plan( + footprints, + Bounds(x_min=x0 - margin, y_min=y0 - margin, x_max=x1 + margin, y_max=y1 + margin), + ) + windows: list[Window] = [] + while True: + m = plan.meta + itemsize = np.result_type(*(t.dtype for t in plan.tiles)).itemsize + need, cap = m.rows * m.cols * (itemsize + 2), physical_memory() // 2 + if need > cap: + raise CatchmentError( + f"the window is {m.rows} x {m.cols} nodes, {need} bytes to flood, over the cap " + f"of half the physical memory ({cap} bytes)" + ) + tile = assemble(plan, repository.load).tile + seed = _seed_mask(m, lake, (float(x), float(y))) + t0 = time.perf_counter() + out = upstream(to_core(tile), seed) + seconds = time.perf_counter() - t0 + if out.touches_nodata: + raise CatchmentError("the catchment reaches NoData in the DEM: it is truncated") + if out.nodes_in == 0: + raise CatchmentError("no seed node has data: the seed lies on NoData") + # Grow while the in-nodes' bounds plus the base margin reach past the + # window; the step itself takes the margin doubled at each growth. + grown = _grown(plan, _plan(footprints, _joined(m, out, float(WINDOW_MARGIN_M)))) + windows.append(Window(Bounds(**_bounds_of(m)), m.rows, m.cols, seconds, grown)) + if not grown: + break + margin *= 2 + plan = _plan(footprints, _joined(m, out, margin)) + if out.touches_edge: + sides = _edge_sides(m, out) + raise CatchmentError(f"the catchment is cut by the data's {' and '.join(sides)} edge") + return _outline(dem_crs, lake, (float(x), float(y)), out, m, seed, tuple(windows)) + + +def _lake(request: CatchmentRequest, dem_crs: str) -> Polygon | None: + """The one lake part containing the seed point, in the DEM's CRS.""" + if request.lakes is None or request.lakes_crs is None: + return None + ((px, py),) = reprojector(request.seed_crs, request.lakes_crs)([request.seed]) + point = Point(px, py) + found = [p for g in request.lakes for p in shapely.get_parts(g) if p.contains(point)] + where = f"({px!r}, {py!r}) in {request.lakes_crs}" + if not found: + raise CatchmentError(f"the seed point {where} is in no lake") + if len(found) > 1: + raise CatchmentError(f"the seed point {where} is in {len(found)} lakes; give one") + move = reprojector(request.lakes_crs, dem_crs) + lake = shapely.transform(found[0], move) + if not np.isfinite(shapely.get_coordinates(lake)).all(): + raise CatchmentError(f"the lake at {where} has a vertex with no image in {dem_crs}") + assert isinstance(lake, Polygon) + return lake + + +def _plan(footprints: Any, bounds: Bounds) -> MosaicPlan: + try: + return plan_mosaic(footprints, bounds) + except MosaicError as exc: + raise CatchmentError(str(exc)) from exc + + +def _seed_mask(m: RasterMeta, lake: Polygon | None, point: tuple[float, float]) -> Any: + """The lake's nodes (`shapely.contains_xy`, over its bounding box only), + or the node nearest the point.""" + seed = np.zeros((m.rows, m.cols), dtype=np.uint8) + if lake is None: + r = round((m.y_max - point[1]) / m.delta_y) + c = round((point[0] - m.x_min) / m.delta_x) + if not (0 <= r < m.rows and 0 <= c < m.cols): + raise CatchmentError(f"the seed point {point} is outside the DEM") + seed[r, c] = 1 + return seed + x0, y0, x1, y1 = lake.bounds + r0, r1 = ( + max(0, math.floor((m.y_max - y1) / m.delta_y)), + min(m.rows - 1, math.ceil((m.y_max - y0) / m.delta_y)), + ) + c0, c1 = ( + max(0, math.floor((x0 - m.x_min) / m.delta_x)), + min(m.cols - 1, math.ceil((x1 - m.x_min) / m.delta_x)), + ) + if r0 > r1 or c0 > c1: + return seed + r, c = np.indices((r1 - r0 + 1, c1 - c0 + 1)) + x, y = m.x_min + (c + c0) * m.delta_x, m.y_max - (r + r0) * m.delta_y + seed[r0 : r1 + 1, c0 : c1 + 1] = shapely.contains_xy(lake, x, y) + return seed + + +def _extent(m: RasterMeta) -> tuple[float, float, float, float]: + return ( + m.x_min, + m.y_max - (m.rows - 1) * m.delta_y, + m.x_min + (m.cols - 1) * m.delta_x, + m.y_max, + ) + + +def _bounds_of(m: RasterMeta) -> dict[str, float]: + return dict(zip(("x_min", "y_min", "x_max", "y_max"), _extent(m), strict=True)) + + +def _joined(m: RasterMeta, out: UpstreamOutcome, margin: float) -> Bounds: + """The window joined with the in-nodes' bounds grown by `margin`.""" + x0, y0, x1, y1 = _extent(m) + return Bounds( + x_min=min(m.x_min + out.col_min * m.delta_x - margin, x0), + y_min=min(m.y_max - out.row_max * m.delta_y - margin, y0), + x_max=max(m.x_min + out.col_max * m.delta_x + margin, x1), + y_max=max(m.y_max - out.row_min * m.delta_y + margin, y1), + ) + + +def _grown(old: MosaicPlan, new: MosaicPlan) -> tuple[str, ...]: + """The sides on which `new`'s window reaches past `old`'s.""" + a, b = old.window, new.window + past = ( + b.row0 < a.row0, + b.col0 < a.col0, + b.row0 + b.rows > a.row0 + a.rows, + b.col0 + b.cols > a.col0 + a.cols, + ) + return tuple(side for side, p in zip(_SIDES, past, strict=True) if p) + + +def _edge_sides(m: RasterMeta, out: UpstreamOutcome) -> tuple[str, ...]: + """The window's sides with an in-node in their first two node lines.""" + near = ( + out.row_min <= 1, + out.col_min <= 1, + out.row_max + 2 >= m.rows, + out.col_max + 2 >= m.cols, + ) + return tuple(side for side, p in zip(_SIDES, near, strict=True) if p) + + +def _outline( + dem_crs: str, + lake: Polygon | None, + point: tuple[float, float], + out: UpstreamOutcome, + m: RasterMeta, + seed: Any, + windows: tuple[Window, ...], +) -> Catchment: + """The outer ring round the seed point (the pour node without a lake), + holes filled and every other ring dropped, counted.""" + mask = np.asarray(out.mask) + t0 = time.perf_counter() + lattice = trace(mask) + rings = [ + np.column_stack([m.x_min + r[:, 1] * m.delta_x, m.y_max - r[:, 0] * m.delta_y]) + for r in lattice + ] + areas = [_signed_area(r) for r in rings] + if lake is None: + r, c = np.argwhere(seed)[0] + point = (m.x_min + c * m.delta_x, m.y_max - r * m.delta_y) + inside = Point(point) + around = [k for k, a in enumerate(areas) if a > 0 and Polygon(rings[k]).contains(inside)] + if not around: + raise CatchmentError(f"no ring of the catchment contains the seed point {point}") + chosen = max(around, key=lambda k: areas[k]) + fine = Polygon(rings[chosen]) + shapely.prepare(fine) + within = [k != chosen and fine.contains(Point(rings[k][0])) for k in range(len(rings))] + dropped = [k for k, a in enumerate(areas) if a > 0 and k != chosen and not within[k]] + holes = [k for k, a in enumerate(areas) if a < 0 and within[k]] + return Catchment( + fine=fine, + crs=crs_label(dem_crs), + seed=point, + lake_area=None if lake is None else lake.area, + nodes=int(out.nodes_in), + seed_nodes=int(np.count_nonzero(seed)), + fine_area=fine.area, + rings_dropped=len(dropped), + dropped_nodes=_nodes_in(mask, [lattice[k] for k in dropped]), + holes_filled=len(holes), + holes_area=-sum(areas[k] for k in holes), + windows=windows, + mask=mask, + meta=m, + trace_seconds=time.perf_counter() - t0, + ) + + +def _signed_area(ring: npt.NDArray[np.float64]) -> float: + """Shoelace: positive counter-clockwise.""" + x, y = ring[:, 0], ring[:, 1] + return 0.5 * float(np.dot(x, np.roll(y, -1)) - np.dot(np.roll(x, -1), y)) + + +def _nodes_in(mask: Any, rings: list[npt.NDArray[np.float64]]) -> int: + """In-nodes inside the dropped outer rings (in lattice units), for the report.""" + count = 0 + for ring in rings: + # Only the nodes in the ring's box: its vertices are on half-cells. + (r0, c0), (r1, c1) = np.ceil(ring.min(axis=0)), np.floor(ring.max(axis=0)) + r, c = np.nonzero(mask[int(r0) : int(r1) + 1, int(c0) : int(c1) + 1]) + count += int(shapely.contains_xy(Polygon(ring), r + r0, c + c0).sum()) + return count diff --git a/src_python/tin_engine/cli.py b/src_python/tin_engine/cli.py index 0dbb0d38..352f007d 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 @@ -68,9 +69,10 @@ sample, triangulate, ) +from tin_engine.catchment import CatchmentRequest, delineate from tin_engine.chains import start_chains from tin_engine.crs import crs_label, parse_crs, transform_description -from tin_engine.dem_input import DemInput, DemRequest, open_dem +from tin_engine.dem_input import DemInput, DemRequest, open_dem, repository_for from tin_engine.domain import DomainError, DomainPolygon, read_domain from tin_engine.elevation import Trimmed, trim from tin_engine.feature_input import ( @@ -80,6 +82,7 @@ FeatureSet, FeatureSource, open_features, + read_lakes, ) from tin_engine.features import DEFAULT_VOCABULARY from tin_engine.grid_domain import default_stride, refine_start_stride, subsample @@ -1280,3 +1283,122 @@ def _off_node(xy: npt.NDArray[np.float64], meta: RasterMeta) -> int: meta.y_max - row * meta.delta_y == xy[:, 1] ) return int(np.count_nonzero(~node)) + + +#: What the suffix of ``catchment --out`` may be: GeoJSON, which ``--domain`` reads. +CATCHMENT_SUFFIXES = (".geojson", ".json") + + +@app.command() +def catchment( + out: Annotated[ + Path, typer.Option("--out", help="Where to write the polygon: .geojson or .json.") + ], + dem: Annotated[ + list[Path], + typer.Option("--dem", help="A GeoTIFF DEM, several (repeat --dem), or one directory."), + ], + seed: Annotated[ + tuple[float, float], + typer.Option("--seed", metavar="X Y", help="A point in the lake, in --seed-crs."), + ], + seed_crs: Annotated[ + str, typer.Option("--seed-crs", help="The seed's CRS, anything pyproj reads (LON LAT).") + ] = "EPSG:4326", + lakes: Annotated[ + Path | None, + typer.Option( + "--lakes", + help="Lake polygons (.gpkg, .geojson or .json); the one under the seed is the " + "seed. Without it, the seed is the DEM node nearest the point.", + ), + ] = None, + lakes_layer: Annotated[ + str | None, typer.Option("--lakes-layer", help="The GeoPackage's features table.") + ] = None, + out_parent: Annotated[ + Path | None, + typer.Option("--out-parent", help="Refuse any output path resolving outside this."), + ] = None, +) -> None: + """Write the catchment of a lake, from the DEM, as a GeoJSON polygon in the + DEM's CRS that ``mesh --domain`` reads (increment 22, PR 1: the fine + outline, drawn between DEM nodes, holes filled). A catchment cut by the + data's edge or by NoData is refused, and nothing is written.""" + if lakes is None and lakes_layer is not None: + raise typer.BadParameter("applies only with --lakes", param_hint="--lakes-layer") + if out.suffix.lower() not in CATCHMENT_SUFFIXES: + raise typer.BadParameter(f"use {' or '.join(CATCHMENT_SUFFIXES)}", param_hint="--out") + target = _destination(out, out_parent, out.stem) + try: + found = None if lakes is None else read_lakes(lakes, lakes_layer, seed, seed_crs) + except FeatureError as exc: + raise typer.BadParameter(str(exc), param_hint="--lakes") from exc + try: + repository, _ = repository_for(tuple(dem)) + request = CatchmentRequest( + seed=seed, + seed_crs=seed_crs, + lakes=None if found is None else found[0], + lakes_crs=None if found is None else found[1], + ) + result = delineate(request, repository) + except OSError as exc: + raise typer.BadParameter(f"cannot read {exc.filename}: {exc}", param_hint="--dem") from exc + except ValueError as exc: + raise typer.BadParameter(_words(exc), param_hint="--dem") from exc + for k, w in enumerate(result.windows, 1): + b = w.bounds + typer.echo( + f"window {k}: x {b.x_min:.0f}-{b.x_max:.0f}, y {b.y_min:.0f}-{b.y_max:.0f}, " + f"{w.rows} x {w.cols} nodes, flood {w.seconds:.2f} s, " + + (f"grown: touches {', '.join(w.grown)}" if w.grown else "contained"), + err=True, + ) + cell = result.meta.delta_x * result.meta.delta_y + if result.lake_area is None: + typer.echo( + f"seed: the pour node at {result.seed} (a pour point must lie on the flow line; " + "it is not snapped)", + err=True, + ) + else: + typer.echo( + f"seed: lake of {result.lake_area / 1e6:.6f} km2, {result.seed_nodes} nodes", err=True + ) + typer.echo( + f"catchment: {result.nodes} nodes, {result.nodes * cell / 1e6:.6f} km2 of node area", + err=True, + ) + vertices = len(result.fine.exterior.coords) - 1 + typer.echo( + f"fine outline: {vertices} vertices, {result.fine_area / 1e6:.6f} km2, " + f"{result.rings_dropped} rings dropped ({result.dropped_nodes} nodes), " + f"{result.holes_filled} holes filled ({result.holes_area / 1e6:.6f} km2), " + f"traced in {result.trace_seconds:.2f} s", + err=True, + ) + properties = { + "seed": list(seed), + "seed_crs": seed_crs, + "nodes": result.nodes, + "fine_vertices": vertices, + "fine_area_m2": result.fine_area, + "windows": [[w.rows, w.cols] for w in result.windows], + } + doc = { + "type": "FeatureCollection", + "crs": {"type": "name", "properties": {"name": result.crs}}, + "features": [ + { + "type": "Feature", + "properties": properties, + "geometry": { + "type": "Polygon", + "coordinates": [[list(xy) for xy in result.fine.exterior.coords]], + }, + } + ], + } + target.write_text(json.dumps(doc), encoding="utf-8") + typer.echo(f"{target}") diff --git a/src_python/tin_engine/dem_input.py b/src_python/tin_engine/dem_input.py index 01895384..3dd5d16c 100644 --- a/src_python/tin_engine/dem_input.py +++ b/src_python/tin_engine/dem_input.py @@ -68,16 +68,21 @@ class DemInput: seams: tuple[Seam, ...] = () +def repository_for( + sources: tuple[Path, ...], nodata: float | None = None +) -> tuple[TiffDemRepository, str]: + """The repository over one directory of tiles or over the files given, + and the run's name: the directory's name, or the first file's stem.""" + first = sources[0] + if first.is_dir(): + return TiffDemRepository.from_directory(first, nodata=nodata), first.resolve().name + return TiffDemRepository(sources, nodata=nodata), first.stem + + def open_dem(request: DemRequest) -> DemInput: """List, plan and assemble. Every refusal is a `ValueError` (`GeoTiffError`, `MosaicError`) or, for a file that cannot be opened, an `OSError`.""" - first = request.sources[0] - if first.is_dir(): - repository = TiffDemRepository.from_directory(first, nodata=request.nodata) - label = first.resolve().name - else: - repository = TiffDemRepository(request.sources, nodata=request.nodata) - label = first.stem + repository, label = repository_for(request.sources, request.nodata) footprints = repository.footprints() if request.domain is None or not footprints: plan, domain, grown = plan_mosaic(footprints, request.bounds, None), None, None diff --git a/src_python/tin_engine/feature_input.py b/src_python/tin_engine/feature_input.py index c72cc34f..0ab54970 100644 --- a/src_python/tin_engine/feature_input.py +++ b/src_python/tin_engine/feature_input.py @@ -28,6 +28,7 @@ import time from collections.abc import Callable, Iterable, Mapping from contextlib import closing +from dataclasses import dataclass from pathlib import Path from typing import Any, Literal @@ -259,57 +260,22 @@ def __init__(self, domain: DomainPolygon, dem: CRS, vocabulary: EdgeVocabulary) def source(self, source: FeatureSource) -> tuple[str, str | None]: """Read one source and take in its features; its CRS text and layer.""" - path, attribute = source.path, source.class_map.attribute - suffix = path.suffix.lower() - if suffix not in SUFFIXES: - raise FeatureError(f"{path.name}: unknown suffix {suffix or '(none)'}; use {SUFFIXES}") - if source.layer is not None and suffix != ".gpkg": - raise FeatureError(f"{path.name}: a layer applies to a GeoPackage only") + read = read_source( + source.path, + source.layer, + source.class_map.attribute, + lambda own: source_region(self.domain, self.dem, own).bounds, + source.crs, + ) + if read.scanned: + self.scanned.append(str(read.layer)) try: - if suffix == ".gpkg": - with closing(open_geopackage(path)) as conn: - layer = layer_info(conn, source.layer) - if layer.rtree is None: - self.scanned.append(layer.table) - own = self._own(source, layer.crs) - box = source_region(self.domain, self.dem, own).bounds - scale = 1.0 if parse_crs(own).is_geographic else 100_000.0 - found = query_features(conn, layer, box, attribute, scale) - self._take(source, own, ((r.fid, r.geometry, r.value) for r in found)) - return own, layer.table - if suffix == ".gml": - with path.open("rb") as stream: - doc = read_gml(stream, attribute) - if doc.crs is None and source.crs is None: - raise FeatureError(f"{path.name}: no geometry names its srsName; give its CRS") - own = self._own(source, doc.crs or str(source.crs)) - self._take(source, own, ((f.fid, f.geometry, f.value) for f in doc.features)) - return own, None - doc = json.loads(path.read_text()) - try: # a malformed document's structure raises any of these - member = doc.get("crs") - text = member["properties"]["name"] if member else GEOJSON_DEFAULT_CRS - rows = [ - (f.get("id", k), f["geometry"] and shape(f["geometry"]), f["properties"] or {}) - for k, f in enumerate(doc["features"]) - ] - values = [(k, g, p.get(attribute)) for k, g, p in rows] - except (KeyError, TypeError, AttributeError, shapely.errors.ShapelyError) as exc: - raise FeatureError(f"{path.name}: not a GeoJSON FeatureCollection ({exc})") from exc - own = self._own(source, text) - self._take(source, own, values) - return own, None + self._take(source, read.crs, read.rows) except FeatureError: raise - except (OSError, ValueError, KeyError, sqlite3.Error) as exc: - raise FeatureError(f"{path.name}: {exc}") from exc - - def _own(self, source: FeatureSource, text: str) -> str: - """The file's CRS text, checked against the one given, if any.""" - if source.crs is not None and parse_crs(source.crs) != parse_crs(text): - raise FeatureError(f"{source.path.name} is in {text} but the given CRS is {source.crs}") - parse_crs(text) - return text + except (ValueError, KeyError) as exc: + raise FeatureError(f"{source.path.name}: {exc}") from exc + return read.crs, read.layer def _take(self, source: FeatureSource, own: str, rows: Iterable[tuple[Any, Any, Any]]) -> None: name, cmap = source.path.name, source.class_map @@ -360,6 +326,100 @@ def _take(self, source: FeatureSource, own: str, rows: Iterable[tuple[Any, Any, self.seconds += time.perf_counter() - t0 +@dataclass(frozen=True, slots=True) +class SourceRows: + """What :func:`read_source` read: the file's CRS text, the GeoPackage + layer (None otherwise) and whether it was scanned without an R-tree, and + ``(fid, geometry, attribute value)`` per feature, in file order.""" + + crs: str + layer: str | None + scanned: bool + rows: list[tuple[Any, Any, Any]] + + +def read_source( + path: Path, + layer: str | None, + attribute: str | None, + box_for: Callable[[str], tuple[float, float, float, float]], + crs: str | None = None, +) -> SourceRows: + """One source's features, unclipped, in its own CRS, by suffix (R1). + + ``attribute`` None reads a GeoPackage's primary key. ``box_for(own CRS)`` + is the box a GeoPackage's R-tree is queried with. ``crs``, when given, + must be the file's; a GML file naming none takes it. Every refusal is a + :class:`FeatureError` naming the file.""" + suffix = path.suffix.lower() + if suffix not in SUFFIXES: + raise FeatureError(f"{path.name}: unknown suffix {suffix or '(none)'}; use {SUFFIXES}") + if layer is not None and suffix != ".gpkg": + raise FeatureError(f"{path.name}: a layer applies to a GeoPackage only") + try: + if suffix == ".gpkg": + with closing(open_geopackage(path)) as conn: + info = layer_info(conn, layer) + own = _own(path, crs, info.crs) + scale = 1.0 if parse_crs(own).is_geographic else 100_000.0 + found = query_features(conn, info, box_for(own), attribute or info.pk, scale) + rows = [(r.fid, r.geometry, r.value) for r in found] + return SourceRows(own, info.table, info.rtree is None, rows) + if suffix == ".gml": + with path.open("rb") as stream: + doc = read_gml(stream, attribute or "") + if doc.crs is None and crs is None: + raise FeatureError(f"{path.name}: no geometry names its srsName; give its CRS") + own = _own(path, crs, doc.crs or str(crs)) + return SourceRows( + own, None, False, [(f.fid, f.geometry, f.value) for f in doc.features] + ) + doc = json.loads(path.read_text()) + try: # a malformed document's structure raises any of these + member = doc.get("crs") + text = member["properties"]["name"] if member else GEOJSON_DEFAULT_CRS + rows = [ + (f.get("id", k), f["geometry"] and shape(f["geometry"]), f["properties"] or {}) + for k, f in enumerate(doc["features"]) + ] + values = [(k, g, p.get(attribute)) for k, g, p in rows] + except (KeyError, TypeError, AttributeError, shapely.errors.ShapelyError) as exc: + raise FeatureError(f"{path.name}: not a GeoJSON FeatureCollection ({exc})") from exc + return SourceRows(_own(path, crs, text), None, False, values) + except FeatureError: + raise + except (OSError, ValueError, KeyError, sqlite3.Error) as exc: + raise FeatureError(f"{path.name}: {exc}") from exc + + +def _own(path: Path, given: str | None, text: str) -> str: + """The file's CRS text, checked against the one given, if any.""" + if given is not None and parse_crs(given) != parse_crs(text): + raise FeatureError(f"{path.name} is in {text} but the given CRS is {given}") + parse_crs(text) + return text + + +def read_lakes( + path: Path, layer: str | None, point: tuple[float, float], point_crs: str +) -> tuple[tuple[BaseGeometry, ...], str]: + """Increment 22: every polygon or multipolygon of a lake source near + ``point`` (in ``point_crs``), in the source's own CRS, with its CRS text. + A GeoPackage is queried with the point's box; lines and points are skipped.""" + + def box_for(own: str) -> tuple[float, float, float, float]: + ((x, y),) = reprojector(point_crs, own)([point]) + return (float(x), float(y), float(x), float(y)) + + read = read_source(path, layer, None, box_for) + lakes = tuple( + shapely.force_2d(g) + for _, g, _ in read.rows + if g is not None and not g.is_empty and g.geom_type in ("Polygon", "MultiPolygon") + ) + return lakes, read.crs + + def _finite(geometry: BaseGeometry) -> bool: return bool(np.isfinite(shapely.get_coordinates(geometry)).all()) diff --git a/src_python/tin_engine/outline.py b/src_python/tin_engine/outline.py new file mode 100644 index 00000000..22ef0621 --- /dev/null +++ b/src_python/tin_engine/outline.py @@ -0,0 +1,83 @@ +"""The fine outline of a node mask: marching squares (increment 22, PR 1). + +`docs/increments/22-auto-catchment.md`, "The fine outline". Every vertex is +the midpoint of a lattice edge between an in-node and an out-node; in each +square of four nodes a segment joins two such midpoints with the in-nodes on +its left in the world frame (x = col east, y = -row north). The one ambiguous +square, two in-nodes on a diagonal, keeps them joined: the catchment is +8-connected and the outside 4-connected (Kong and Rosenfeld 1989). The mask is +padded by one row and column of out-nodes, so every ring closes. + +Each boundary midpoint starts exactly one segment and ends exactly one, so +the segments close into simple rings that share no point; outer rings are +counter-clockwise in the world frame and holes clockwise. + +Numpy only: no `_core`, no paths. +""" + +from __future__ import annotations + +from typing import Any + +import numpy as np +import numpy.typing as npt + +Ring = npt.NDArray[np.float64] + +#: A square's corners counter-clockwise in the world frame, as (row, col) +#: offsets from its north-west node: south-west, south-east, north-east, +#: north-west. Edge k runs from corner k to corner k + 1. +_CORNERS = ((1, 0), (1, 1), (0, 1), (0, 0)) + + +def trace(mask: npt.ArrayLike) -> list[Ring]: + """Every ring of `mask`'s in-nodes (non-zero), as `(k, 2)` float arrays of + `(row, col)` in the mask's own indices, each vertex a half on one axis. + Not closed: the first vertex is not repeated.""" + padded = np.pad(np.asarray(mask) != 0, 1) + rows, cols = padded.shape + full = [padded[r : rows - 1 + r, c : cols - 1 + c] for r, c in _CORNERS] + count: npt.NDArray[np.uint8] = np.sum([f.astype(np.uint8) for f in full], axis=0) + # Only squares with in- and out-corners carry a segment. + i, j = np.nonzero((count > 0) & (count < 4)) + corner = [f[i, j] for f in full] + # A midpoint's key is its doubled padded (row, col), flattened: unique. + width = 2 * cols + mid = [ + (2 * i + _CORNERS[k][0] + _CORNERS[(k + 1) % 4][0]) * width + + 2 * j + + _CORNERS[k][1] + + _CORNERS[(k + 1) % 4][1] + for k in range(4) + ] + starts: list[Any] = [] + ends: list[Any] = [] + for k in range(4): + # Edge k leaves the in-nodes (corner k in, k + 1 out): the segment + # ends at the next edge counter-clockwise that enters them again. + leaving = corner[k] & ~corner[(k + 1) % 4] + after = corner[(k + 2) % 4] + end = np.where( + after, + mid[(k + 1) % 4], + np.where(corner[(k + 3) % 4], mid[(k + 2) % 4], mid[(k + 3) % 4]), + ) + starts.append(mid[k][leaving]) + ends.append(end[leaving]) + start, stop = np.concatenate(starts), np.concatenate(ends) + order = np.argsort(start, kind="stable") + following = order[np.searchsorted(start, stop, sorter=order)] + seen = np.zeros(len(start), dtype=bool) + rings: list[Ring] = [] + for first in order: + if seen[first]: + continue + loop = [] + s = int(first) + while not seen[s]: + seen[s] = True + loop.append(start[s]) + s = int(following[s]) + keys = np.asarray(loop, dtype=np.int64) + rings.append(np.column_stack([keys // width / 2 - 1, keys % width / 2 - 1])) + return rings diff --git a/tests/cpp/CMakeLists.txt b/tests/cpp/CMakeLists.txt index f8c91a8d..57e27d1a 100644 --- a/tests/cpp/CMakeLists.txt +++ b/tests/cpp/CMakeLists.txt @@ -316,20 +316,4 @@ add_terrain_backend_test(prop_noding_verifier_sweep property/prop_noding_verifie # (docs/increments/22-auto-catchment.md, "Flow and membership: one flood" and # "The red suites"). Not named invariant-critical; the descent oracle and the # invariants on random DEMs are the defence. It starts no threads. -# -# RED UNTIL hydrology/upstream.hpp EXISTS (increment 12's scaffold, 39dba00): -# while the header is absent the target is EXCLUDE_FROM_ALL and not registered -# with ctest, so `cmake --build build` and `ctest` stay green; see the red with -# `cmake --build build --target test_hydrology_upstream`. Once the header lands, -# re-running `cmake -S . -B build` registers the suite with no edit here. -# REMOVE AT GREEN: replace the if/else with the plain add_terrain_test line. -if(EXISTS "${PROJECT_SOURCE_DIR}/include/terrain/hydrology/upstream.hpp") - add_terrain_test(test_hydrology_upstream unit/test_hydrology_upstream.cpp) -else() - message(STATUS "test_hydrology_upstream: include/terrain/hydrology/upstream.hpp absent; " - "suite built only on request (increment 22 red step)") - add_executable(test_hydrology_upstream EXCLUDE_FROM_ALL unit/test_hydrology_upstream.cpp) - target_link_libraries(test_hydrology_upstream PRIVATE - terrain_headers terrain_test_support Catch2::Catch2WithMain) - target_compile_options(test_hydrology_upstream PRIVATE -Wall -Wextra -Wpedantic -Werror) -endif() +add_terrain_test(test_hydrology_upstream unit/test_hydrology_upstream.cpp) From 0e1c29a223f3fd1898f4d36820e6147e685e0a27 Mon Sep 17 00:00:00 2001 From: Ola Skavhaug Date: Tue, 29 Sep 2026 01:55:20 +0200 Subject: [PATCH 09/14] 22 design: the window loop checks the base margin Co-Authored-By: Claude Opus 5.5 --- docs/increments/22-auto-catchment.md | 47 ++++++++++++++++++++++------ 1 file changed, 37 insertions(+), 10 deletions(-) diff --git a/docs/increments/22-auto-catchment.md b/docs/increments/22-auto-catchment.md index f18792ce..299bcc33 100644 --- a/docs/increments/22-auto-catchment.md +++ b/docs/increments/22-auto-catchment.md @@ -394,21 +394,38 @@ re-plan pattern as 15b's `_domain_plan`: a. If `touches_nodata`, refuse: the catchment is truncated by missing data, which no growth can fix. - b. The need N is the in-nodes' bounds grown by the margin, clamped to E. - A side of W *can grow* when N reaches past W on that side; since N is - clamped to E, that also means W is not yet at E there. - c. If some side can grow: the next window is W joined with N, the margin - doubles, and the loop repeats from 2. + b. The need N is the in-nodes' bounds grown by the **base** margin (2000 + m, never doubled), clamped to E. A side of W *can grow* when N reaches + past W on that side; since N is clamped to E, that also means W is not + yet at E there. + c. If some side can grow: the step margin doubles (4000 m at the first + growth, then 8000 m, ...), the next window is W joined with the + in-nodes' bounds grown by the step margin, clamped to E, and the loop + repeats from 2. Only the step doubles; the test in b always uses the + base margin. + + Revised 2026-09-29 at the PR 1 green step (58f6904): the previous + wording let the test in b use the doubled margin too, so the need grew + as fast as the window, every step grew, and the loop only stopped at the + data's edge. @developer's first Bygdin run took 8 minutes and was then + refused by 15a's coverage check. With the base margin in b, Bygdin gives + 304.91 km² (-0.21 % against NVE) in 6.3 s over 3 windows. d. Otherwise decide by the flag. If `touches_edge`, refuse, naming the sides (from the bounds: an in-node in the first or last two rows or columns), each of which is then at E: the catchment is cut by the data's edge. If not, accept. - **It terminates.** Step c runs only when N reaches past W on some side, - and the next window contains N, so it has at least one more row or column - than W; every window lies inside E, which is finite. So step c runs at - most (rows of E + columns of E) times, and in practice a handful, since - the margin doubles. The memory cap (step 5) may refuse earlier. Every + **It terminates.** Step c runs only when N reaches past W on some side. + The next window contains the in-nodes' bounds grown by the step margin, + which is at least the base margin, so it contains N (both clamped to E) + and has at least one more row or column than W; every window lies inside + E, which is finite. So step c runs at most (rows of E + columns of E) + times. **And it stops early**: after a growth step the window holds the + last catchment's bounds plus the base margin, so a further step happens + only if the catchment itself grew past that in the new window. On a + catchment far from the data's edge the number of windows is one more + than the number of times the catchment outgrew its window's base margin, + in practice two or three. The memory cap (step 5) may refuse earlier. Every exit is an accept or a refusal from a or d. **Accepting means**: the catchment comes no closer than two nodes to any @@ -417,6 +434,16 @@ re-plan pattern as 15b's `_domain_plan`: does not. The margin is a heuristic against the known limit below; the flag is the rule. + **What a test should pin** (for @tester; described only): on a + synthetic DEM whose data extends more than four base margins beyond the + full catchment on every side, and whose catchment is larger than the + first window (the bowl fixture, placed in a larger raster), the loop + accepts in at most 3 windows, and the final window reaches the data's + edge on no side. A second, cheaper pin: a catchment that lies inside the + first window with the base margin to spare is accepted in exactly one + window. The first of these fails on the pre-58f6904 loop, which grows + until it meets the data's edge. + Refusing a truncated catchment stays the default (main session / @architect, 2026-09-29), for Ola to confirm; alternative: write it with a warning under an `--allow-truncated` flag. From 608e3665df87e619f8693c762d5d1fa75a61131e Mon Sep 17 00:00:00 2001 From: Ola Skavhaug Date: Tue, 29 Sep 2026 01:57:00 +0200 Subject: [PATCH 10/14] 22 PR1 tests: the window loop stops early The two pins the amended design asks for (0e1c29a, "What a test should pin"): - test_the_loop_stops_early_far_from_the_data_edge: the bowl, with more than four base margins of data beyond the catchment on every side and a catchment larger than the first window, is accepted in 2 or 3 windows, the last reaching the data's edge on no side, equal to one flood. - test_a_catchment_inside_the_first_window_takes_one_window: a lake on a peak, 49 nodes, accepted in one window. Both pass at 58f6904. Against the old rule (the check on the doubled margin, the step on the undoubled one: a scratchpad copy of catchment.py with those two lines changed, src_python untouched) the first fails: 6 windows, 49x49 up to the whole 400x200 raster. The second passes there too; it is the cheaper pin the design describes, guarding the other direction. Co-Authored-By: Claude Opus 5.5 --- tests/python/test_catchment.py | 39 ++++++++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) diff --git a/tests/python/test_catchment.py b/tests/python/test_catchment.py index e2aa8223..df8697e3 100644 --- a/tests/python/test_catchment.py +++ b/tests/python/test_catchment.py @@ -284,6 +284,45 @@ def test_the_memory_cap_refuses_before_the_flood(api: Any, monkeypatch: pytest.M api.delineate(request(api), repository_of(bowl())) +def test_the_loop_stops_early_far_from_the_data_edge(api: Any) -> None: + """The design's pin after 0e1c29a: grow-or-stop is decided on the base + margin, and only the step doubles. The loop that checked the doubled + margin grew on every step until it met the data's edge.""" + z = bowl() + tile = tile_of(z) + reference = full_flood(tile, seeds_in(tile, lake_box())) + r, c = np.nonzero(reference) + rows, cols = z.shape + # The fixture's premises: more than four base margins of data beyond the + # catchment on every side, and a catchment larger than the first window. + assert min(r.min(), c.min(), rows - 1 - r.max(), cols - 1 - c.max()) > 4 * MARGIN_CELLS + assert r.min() < 200 - 3 - MARGIN_CELLS + + result = api.delineate(request(api), repository_of(z)) + assert 2 <= len(result.windows) <= 3 + m = result.meta + assert m.x_min > X0 + assert m.y_max < Y0 + assert m.x_min + (m.cols - 1) * D < X0 + (cols - 1) * D + assert m.y_max - (m.rows - 1) * D > Y0 - (rows - 1) * D + assert np.array_equal(on_whole(result, z.shape), reference) + + +def test_a_catchment_inside_the_first_window_takes_one_window(api: Any) -> None: + """A peak with the lake on top: nothing drains into the lake, so the + catchment is the lake's 49 nodes, inside the first window (the lake's box + plus the base margin) with the margin to spare.""" + r, c = np.indices((100, 100)).astype(np.float64) + z = (100.0 - np.hypot(r - 50, c - 50)).astype(np.float32) + lake = lake_box(50, 50) + result = api.delineate( + api.CatchmentRequest(seed=lat(50, 50), seed_crs=EPSG, lakes=(lake,), lakes_crs=EPSG), + repository_of(z), + ) + assert result.nodes == result.seed_nodes == 49 + assert len(result.windows) == 1 + + # --------------------------------------------------------------------------- # Holes filled, other pieces dropped # --------------------------------------------------------------------------- From 0bbdbd0c7509f0e544d6b189d900c96de40aa1ab Mon Sep 17 00:00:00 2001 From: Ola Skavhaug Date: Tue, 29 Sep 2026 07:55:08 +0200 Subject: [PATCH 11/14] 22 PR1 docs: as built at PR 1 Co-Authored-By: Claude Opus 5.5 --- docs/increments/22-auto-catchment.md | 56 ++++++++++++++++------------ project_structure.md | 35 +++++++++++++---- 2 files changed, 60 insertions(+), 31 deletions(-) diff --git a/docs/increments/22-auto-catchment.md b/docs/increments/22-auto-catchment.md index 299bcc33..a64dfddc 100644 --- a/docs/increments/22-auto-catchment.md +++ b/docs/increments/22-auto-catchment.md @@ -1,7 +1,9 @@ # Increment 22 — auto-catchment: the catchment of a lake, from the DEM (Bygdin first) -Status: **designed** (`@architect`, 2026-09-29), on branch -`increment22-autocatchment` off master `b4847d7`. Nothing is built. Ola was +Status: **PR 1 built, under review** (2026-09-29). Designed by `@architect` +on branch `increment22-autocatchment` off master `b4847d7`; PR 1 built at +`58f6904` (green), under review; PR 2 (the outline reduction) designed, in a +stacked PR ("As built" below covers PR 1 only). Ola was asleep while this was written; every choice he would normally make is marked "Default (main session / @architect, 2026-09-29), for Ola to confirm", with the alternative, so the loop can run tonight. @@ -132,9 +134,9 @@ design leans on a detail beyond the abstract, it says so. flooded it. - **O'Callaghan and Mark 1984**, "The extraction of drainage networks from digital elevation data", *Computer Vision, Graphics, and Image Processing*, - doi:10.1016/S0734-189X(84)80047-X. Crossref lists it as volume 27(2), page - 247; it is commonly cited as 28(3):323-344. Cite it by DOI until someone - resolves which. D8: each node drains to its steepest neighbour, slope + 28(3):323-344, doi:10.1016/S0734-189X(84)80011-0 (checked on Crossref, + 2026-09-29, at review; the DOI first cited here, 80047-X, is a one-page + item in 27(2):247). D8: each node drains to its steepest neighbour, slope measured with the diagonal's length. **Departure, and why:** the flood drains each node to its *lowest* filled neighbour, not its steepest, so the diagonal distance plays no part. The two differ only where a diagonal @@ -664,21 +666,27 @@ request and the result are frozen; the CLI is the only place with paths. ### New and changed files -| File | What | Production lines (estimate) | -|---|---|---| -| `include/terrain/hydrology/upstream.hpp` | the flood, `UpstreamOutcome` | 110 | -| `include/terrain/vector_simplify/area_collapse.hpp` | the reduction, `ReduceOutcome`, edge grid | 260 | -| `bindings/core.cpp` | `upstream`, `reduce_ring`, two outcome classes | 80 | -| `src_python/tin_engine/_core.pyi` | their stubs | 30 | -| `src_python/tin_engine/outline.py` | the tracer | 70 | -| `src_python/tin_engine/catchment.py` | request, seed, window loop, result | 170 | -| `src_python/tin_engine/dem_input.py` | repository helper split out | 10 | -| `src_python/tin_engine/feature_input.py` | `read_source` split out of `_Tally.source`, `read_lakes` | 30 | -| `src_python/tin_engine/cli.py` | `catchment` command, report, writer | 100 | -| `project_structure.md` | the two C++ modules and two Python modules | docs | - -About 860 lines, over the 700 ceiling (CLAUDE.md §2), so two PRs on this -branch. +| File | What | Estimate | As built, PR 1 | +|---|---|---|---| +| `include/terrain/hydrology/upstream.hpp` | the flood, `UpstreamOutcome` | 110 | 106 | +| `include/terrain/vector_simplify/area_collapse.hpp` | the reduction, `ReduceOutcome`, edge grid | 260 | | +| `bindings/core.cpp` | `upstream`, `reduce_ring`, two outcome classes | 80 | 48 | +| `src_python/tin_engine/_core.pyi` | their stubs | 30 | 19 | +| `src_python/tin_engine/outline.py` | the tracer | 70 | 50 | +| `src_python/tin_engine/catchment.py` | request, seed, window loop, result | 170 | 244 | +| `src_python/tin_engine/dem_input.py` | repository helper split out | 10 | 1 | +| `src_python/tin_engine/feature_input.py` | `read_source` split out of `_Tally.source`, `read_lakes` | 30 | 38 net | +| `src_python/tin_engine/cli.py` | `catchment` command, report, writer | 100 | 113 | +| **Total** | | about 860 | **619 net** (674 added, 55 removed) | +| `project_structure.md` | the two C++ modules and two Python modules | docs | | + +About 860 lines estimated, over the 700 ceiling (CLAUDE.md §2), so two PRs. +"As built" is `@reviewer`'s count at review (2026-09-29), by CLAUDE.md §2's +rule (blank lines, comments and docstrings not counted), over PR 1 = +`master..608e366`. The totals are `@reviewer`'s; the per-file numbers are +theirs too and were not recounted here. The PR 1 column, counted at +`58f6904`, sums to 619. PR 1 is under the 700 ceiling. PR 2's figures are the +estimates until it is built. ### The PR split @@ -690,7 +698,9 @@ branch. binding, and the command reducing by default with `--outline-tolerance`. Red, green, review, then the Bygdin acceptance run. -Both go on `increment22-autocatchment`, PR 2 stacked on PR 1. Neither touches +Both go on `increment22-autocatchment`, PR 2 stacked on PR 1. As built: PR 1 +is `master..608e366` (red `1e1b3bb`, green `58f6904`, the window test +`608e366`); PR 2 is designed, in a stacked PR. Neither touches refine or mesh code, so the 1 m benchmark and scaling sweep (README, rule 2) do not apply; the Bygdin run below is this increment's acceptance. @@ -701,7 +711,7 @@ Lean: no throwaway implementations, no mutation round. Default (main session invariants on random inputs are the defence. None is named invariant-critical for mutation testing tonight. -**PR 1, C++ (Catch2), `tests/cpp/unit/hydrology_upstream.cpp`:** +**PR 1, C++ (Catch2), `tests/cpp/unit/test_hydrology_upstream.cpp`:** - *Brute-force oracle*. For random DEMs with no interior pit and no ties, built as z = 10 x (steps to the nearest edge) + a distinct fraction per @@ -747,7 +757,7 @@ for mutation testing tonight. **PR 2, C++ and Python:** -- `tests/cpp/unit/area_collapse.cpp`: status for a negative, NaN or infinite +- `tests/cpp/unit/test_area_collapse.cpp`: status for a negative, NaN or infinite tolerance, a clockwise ring, fewer than 4 vertices; tolerance 0 removes only collinear vertices; a traced rectangle of nodes reduces to at most 8 vertices at a tolerance of one cell, with its area unchanged; a keep-point diff --git a/project_structure.md b/project_structure.md index 7f5e37ff..26091d90 100644 --- a/project_structure.md +++ b/project_structure.md @@ -60,13 +60,18 @@ include/terrain/ # public C++ headers, header-only where possible refine.hpp # RefineOptions, RefineOutcome, the round loop (14), # Delaunay insertion (14b), the quality-start call # (20), constraint feet (20b) + hydrology/ + upstream.hpp # upstream(z, seed) -> UpstreamOutcome: one + # Priority-Flood labelling the nodes that drain + # into the seed set, plus edge/NoData flags (22) src/ # C++ implementation, one directory per module # (only predicates/ and cdt/ exist; rest planned) predicates/ # exact orient2d/incircle; namespace terrain::pred parallel_util/ # (none: header-only, include/terrain/parallel_util/) vector_simplify/ # Visvalingam-Whyatt, Douglas-Peucker, topology checks - hydrology/ # pit fill, flow direction, accumulation, catchments, streams + # (no hydrology/ here: it is header-only, + # include/terrain/hydrology/) # (no noding/ here, and none planned: the noder is # header-only, its driver a template on the kernel) cdt/ # thin wrapper over vendored Detria @@ -102,7 +107,17 @@ src_python/tin_engine/ # public Python API (distribution name: rasputin) # read, mapped to masks by a ClassMap, pre-clipped, # moved to the DEM's CRS and clipped to the domain # as linework -> FeatureSet (16b); opens GeoJSON - # and .gml itself; never imports _core + # and .gml itself; never imports _core. + # read_source: one file's raw rows (16b's reader, + # shared); read_lakes: the polygons of --lakes + # near the seed point, in the file's CRS (22) + outline.py # trace(mask): marching-squares rings between in- + # and out-nodes, (8, 4) saddle rule; numpy only, + # never imports _core (22) + catchment.py # CatchmentRequest -> delineate(request, repo) -> + # Catchment: seed, the window loop over 15a's + # plan, _core.upstream and the fine ring; takes a + # DemRepository and no path (22) chains.py # start_chains: the domain's rings, then every # feature line, as (indices, role, mask) (16b); # never imports _core @@ -117,7 +132,9 @@ 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 + # one exception: the catchment GeoJSON writer is + # in cli.py (22); see "Python API surface" __init__.py ply.py # arrays -> PLY bytes; takes no path and opens nothing vtk_legacy.py # arrays + EdgeVocabulary -> legacy .vtk bytes, for @@ -159,9 +176,9 @@ docs/increments/ # per-increment design records; see its README ``` raster ────────────────┐ ▼ -predicates ───────────┬─→ vector_simplify hydrology (depends on raster, - │ │ parallel_util) - └─→ noding ←──────────┘ (catchment outline, river polylines) +predicates ───────────┬─→ vector_simplify hydrology (raster only) + │ + └─→ noding │ ▼ cdt (thin wrapper over vendored Detria, @@ -324,7 +341,7 @@ Visvalingam-Whyatt for area-preserving polygon simplification; Douglas-Peucker f ### `hydrology` -Implements `auto_catchments.md`: pit filling (priority-flood with epsilon, plus optional Lindsay hybrid), D8 flow direction, parallel flow accumulation (Barnes-Lehman-Mulla), catchment BFS, stream extraction with Strahler ordering, mask polygonization. May depend on RichDEM (MIT) as either reference or external library — decision deferred to first implementation pass. +Header-only; depends on `raster` only. `upstream.hpp` (increment 22): `upstream(z, seed)`, a template on the `RasterSource` concept. One Priority-Flood (Barnes, Lehman and Mulla 2014) from the window's edge and from the nodes beside NoData, ties first in, first out, labels a node *in* when it is a seed or was flooded from an *in* node. So each node drains to the node that flooded it, its lowest filled neighbour, and there is no separate fill, no flat resolution and no D8 or accumulation pass. Returns the mask, the in-nodes' count and bounds, and `touches_edge` / `touches_nodata` (an in-node that is, or neighbours, an outlet of that kind). Serial; the binding releases the GIL. What `auto_catchments.md` also sketches (epsilon filling, D8, accumulation, streams and Strahler order, RichDEM) is not built; the outline tracer is Python (`outline.py`). See `docs/increments/22-auto-catchment.md`. ### `noding` @@ -387,7 +404,7 @@ CMake-based, building the header-only core plus one Python extension. C++20 requ - **Detria** (header-only) is vendored at a pinned SHA in `lib/detria/`. CMake does **not** search for a system copy: version skew in a geometry kernel across machines is a reproducibility hazard and vendoring a header costs nothing. -- **External deps under consideration:** RichDEM (optional, MIT), Eigen (if linear algebra needs grow beyond what we want to hand-roll). +- **External deps under consideration:** RichDEM (optional, MIT; increment 22 does not use it), Eigen (if linear algebra needs grow beyond what we want to hand-roll). - **No CGAL and no GDAL** in the new core: prohibited by `CLAUDE.md` §2. **No Boost.Geometry** either, which is a scope choice and not a prohibition — reading this line as one is what put an unauthored ban in §2 for six @@ -402,6 +419,8 @@ There is no `setup.py`. It was removed with the foundation reset, along with the The public API is the `tin_engine` package, calling into `tin_engine._core`. `tin_engine.viz` is the renderer that turns a triangulation into an SVG a person can look at; it consumes the `typing.Protocol`s in `viz/protocols.py` and **never imports `_core`**, so it is testable with no compiled extension in the process. `cli.py` is the single composition root that joins the two -- the same shape as the rule below that exactly one module adapts decoded raster data into `_core`. The `rasputin draw` command drives it: it looks a fixture up in `viz.fixtures.GALLERY`, maps that fixture's **string** chain roles onto `_core.ChainRole` -- `viz/` may not name the enum, so the mapping is the composition root's -- runs `build_pslg` and `triangulate`, and hands the fixture itself to `build_scene` as the `PslgLike`, which is what lets a fixture the validator rejects still be drawn. It is also the only place a path exists, and it resolves and refuses one before writing. +**The catchment GeoJSON writer is in `cli.py`** (increment 22, `rasputin catchment`): it builds a one-feature `FeatureCollection` with a `crs` member and writes it, an exception to "all file encoding lives in `io/`". Recommended (@architect, 2026-09-29): move it, as a function from a polygon, its CRS text and its properties to bytes that takes no path and opens nothing, into `io/geojson.py`, the shape `io/ply.py` and `io/vtk_legacy.py` already have; `cli.py` keeps only the write. It is small, and a GUI or API worker writing a catchment would otherwise have to go through the CLI. Until it moves, this paragraph is the record of the exception. + The pre-migration `rasputin.*` modules (`mesh.py`, `geometry.py`, `reader.py`, `tin_repository.py` and friends) are archived under `legacy/rasputin/` rather than kept in place, so this is a re-implementation against the new backend rather than a rewiring of stable modules. Porting proceeds one entry point at a time; shapes worth preserving should be read out of `legacy/` before being reintroduced. ## What gets deleted, eventually From d3265fce8efbaacec6e607262847e4ce89e356a1 Mon Sep 17 00:00:00 2001 From: Ola Skavhaug Date: Tue, 29 Sep 2026 07:55:21 +0200 Subject: [PATCH 12/14] ROADMAP: increment 22 PR 1 built Co-Authored-By: Claude Opus 5.5 --- ROADMAP.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ROADMAP.md b/ROADMAP.md index c70a0225..ae6b9821 100644 --- a/ROADMAP.md +++ b/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 | planned: ruled by Ola 2026-09-28 (16b's Q2) as a separate increment, next after 16b | `docs/increments/16b-terrain-polygons.md` (Q2, R6) | | 16d | A record of the input polylines in the mesh. Today a constraint edge keeps only its type bits; the source feature, its attributes (a river's name or order, a road's class) and which edges belong to one polyline are lost, and where edges from two features merge only the union of their bits survives. Proposed: a feature table (one row per input feature: source ID and attributes) and, per constraint edge, the list of features it came from (a list, since a merged edge belongs to several). The noder already tracks each piece's input chain; the open part is the output: a per-edge list in `.vtk`/`.ply`, or a side file | proposed by Ola 2026-09-28; to be designed by `@architect` once Ola has said what it is for; after 16b | - | | 20c | Soft quality criterion: a penalty that each Steiner node or constraint split must pay for in angle gained, instead of 20's hard 25°; applied at the start and during DEM refinement; may split constraint segments when that improves the mesh. Ola's rulings on 20's C1-C3 | to design after 16b (`@architect` measures cost against 20 first) | `docs/increments/20-start-quality.md` (Ola's rulings) | -| — | Auto-catchment: the watershed upstream of a coordinate, computed from the DEM and handed to `--domain`, so a catchment no longer has to be supplied as a file (Ola, 2026-09-27: "not far into the future"). The textbook route is depression handling (Priority-Flood, Barnes, Lehman and Mulla 2014), D8 flow directions (O'Callaghan and Mark 1984) and accumulation, the pour point snapped to the strongest flow nearby, the upstream cells traced and their outline turned into a polygon; the literature check is `@architect`'s. Open for its design: whether it runs in the C++ core (a 10 m tile is 25 M cells); how a stair-stepped cell outline becomes a domain polygon, which meets input coarsening; **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) | **designed as increment 22** (`@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) | **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 From a612269138f29c64a24921f45afed5308087cf78 Mon Sep 17 00:00:00 2001 From: Ola Skavhaug Date: Tue, 29 Sep 2026 09:12:39 +0200 Subject: [PATCH 13/14] Retrospective agenda: role bleed as the main focus (Ola), with the GNN analogy Co-Authored-By: Claude Opus 5.5 --- docs/retrospectives/next.md | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/docs/retrospectives/next.md b/docs/retrospectives/next.md index e59d0c00..a525c502 100644 --- a/docs/retrospectives/next.md +++ b/docs/retrospectives/next.md @@ -4,7 +4,21 @@ When: after auto-catchment is ready to use (Ola, 2026-09-28). Run by `@orchestrator`. Items are added as they come up; the retrospective itself gets its own dated file here, and this file is then emptied. -## Agents taking on each other's work (Ola, 2026-09-28) +## Agents taking on each other's work (Ola, 2026-09-28): the main focus + +Ola, 2026-09-29: "when we do our retrospective, we should have a specific +focus on the role bleed issue." Names for it: role drift or role bleed in +practice; "disobey role specification" (failure mode 1.2) in the MAST +taxonomy of multi-agent failures (Cemri et al. 2025, arXiv 2503.13657), +which traces most such failures to weak role definitions and missing checks. +Ola's analogy (2026-09-28): over-smoothing in GNNs, where repeated message +passing makes every node look alike. Here every hand-off carries the whole +context, and each persona picks up a bit of the others' jobs until the roles +blur. The GNN remedies map across: skip connections (restate each persona's +role and limits in every brief) and a bounded reach for messages (hard limits +on what each persona may touch). Inside transformers the same effect is +called rank collapse or over-smoothing, which skip connections also counter +(Dong, Cordonnier and Loukas 2021; from memory, not checked). Ola: "One concern I have is that the agents 'leak' responsibilities to one another." All personas are the same model with different briefs, and each From 24797d097827a3081f84e401f795d59e160c08ab Mon Sep 17 00:00:00 2001 From: Ola Skavhaug Date: Tue, 29 Sep 2026 09:48:33 +0200 Subject: [PATCH 14/14] 22 PR1: GCC -Werror fixes (enum in ?:, dangling reference) upstream.hpp: cast kIn to std::uint8_t so both arms of the conditional have the same type (GCC -Wextra: enumerated and non-enumerated type). core.cpp: hold the raster geometry by value rather than as a reference bound through std::visit (GCC -Wdangling-reference false positive). Co-Authored-By: Claude Opus 5.5 --- bindings/core.cpp | 4 ++-- include/terrain/hydrology/upstream.hpp | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/bindings/core.cpp b/bindings/core.cpp index 8e5b51d9..0382a3ae 100644 --- a/bindings/core.cpp +++ b/bindings/core.cpp @@ -984,8 +984,8 @@ whether it may continue past the window's edge or past NoData. [](const BoundRasterView& raster, const py::object& seed) { using U8 = py::array_t; const auto s = U8::ensure(seed); - const auto& g = std::visit([](const auto& v) -> const auto& { return v.geometry(); }, - raster.view); + const auto g = std::visit([](const auto& v) -> const auto& { return v.geometry(); }, + raster.view); if (!s || s.ndim() != 2 || static_cast(s.shape(0)) != g.rows() || static_cast(s.shape(1)) != g.cols()) throw py::value_error(std::format( diff --git a/include/terrain/hydrology/upstream.hpp b/include/terrain/hydrology/upstream.hpp index 144b8ab0..39853cfc 100644 --- a/include/terrain/hydrology/upstream.hpp +++ b/include/terrain/hydrology/upstream.hpp @@ -97,7 +97,7 @@ template each_neighbour(i, rows, cols, [&](std::size_t j) { if (state[j] != kUnreached) return; - state[j] = seed[j] != 0 ? kIn : label; + state[j] = seed[j] != 0 ? static_cast(kIn) : label; const double zj = level_of(j); queue.emplace(zj > level ? zj : level, counter++, j); });