Seed the cost surface only from stream cells (#41)

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: flooded
 Title: Portable Floodplain Delineation from DEM and Stream Network
-Version: 0.4.0
+Version: 0.4.1
 Authors@R: c(
     person("Allan", "Irvine", , "[email protected]", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0002-3495-2128")),

NEWS.md

@@ -1,3 +1,35 @@
+# flooded 0.4.1
+
+- Fix `fl_cost_distance()` seeding on every zero-friction cell rather than only stream cells (#41).
+  Seeds are encoded by setting stream cells to `0` and calling `terra::costDist(target = 0)`, which
+  matches *every* zero cell — so any cell whose friction was already exactly zero acted as a free
+  cost source. Friction exactly equal to `0` is now floored to `1e-6` before seeding. Flat ground
+  stays cheap to cross (0.1 accumulated over a 100 km path at 10 m, against a default
+  `cost_threshold` of 2500); it simply stops being a source. Negative friction is deliberately not
+  floored, so `terra::costDist()`'s own rejection of a negative cost surface is left intact.
+- The fix strictly *removes* spurious reach from the cost mask; it never adds any. Measured on the
+  two DEMs this package ships, and the answer differs by dataset — check your own rather than
+  assuming:
+
+  | DEM | exact-zero slope cells | cost-mask change | delineation change |
+  |---|---|---|---|
+  | bundled `dem.tif` / `slope.tif`, 10 m | 0 of 45,726 | none | none |
+  | `pars_dem.tif` (MRDEM-30, 30 m, 20.9 Mcell) | 80 of 10.7 M | -2,289 cells (214 ha), 0 added | none |
+
+  So MRDEM-30 *does* contain exact zeros, the cost mask *does* move — and on both shipped datasets
+  the delineation does not, because the slope, distance and flood criteria plus morphological
+  cleanup absorb every affected cell. `fl_valley_confine()` returns the same 53,635 cells on the
+  bundled tile and the same 521,028 cells on the Parsnip Watershed Group, with zero cells differing
+  in either direction. The shipped vignette artifacts are therefore still current.
+- Do not read that as a general guarantee. Where cost is the binding criterion — flatter terrain, a
+  laxer `slope_threshold`, a larger `flood_factor` — results will move. Exposure is highest on
+  integer-metre DEMs, hydro-flattened lake surfaces and void-filled plateaus. Check with
+  `sum(terra::values(slope) == 0, na.rm = TRUE)`.
+- The effect is largest under `fl_valley_attribute()`, where cost is what separates one watercourse's
+  floodplain from another's: a flat patch inside a group's corridor would have spread that group's
+  mask across ground its own streams never reach.
+- Fix a stray one-space indent in `fl_valley_poly()`.
+
 # flooded 0.4.0
 
 - New `fl_valley_attribute()` — attribute a finished `fl_valley_confine()` delineation to the stream groups that produced it, so a floodplain can be filtered and queried per watercourse or reach rather than only per network (#40). Returns one `sf` row per group; rows overlap where ground is genuinely shared between watercourses, which near a confluence is most of it. The delineation is never recomputed, so changing the grouping key relabels the output without moving a boundary.

R/fl_cost_distance.R

@@ -20,6 +20,21 @@
 #'
 #' Cells that are `NA` in `friction` are impassable barriers.
 #'
+#' Seeds are encoded by setting stream cells to zero in the friction surface,
+#' which is only unambiguous if no other cell is zero. Friction rasters do
+#' contain exact zeros — integer-metre DEMs, hydro-flattened lake surfaces and
+#' void-filled plateaus all quantize to perfectly flat — so cells with friction
+#' exactly `0` are floored to `1e-6` before seeding. Flat ground therefore
+#' remains cheap to cross but is no longer a cost source: at 10 m resolution a
+#' 100 km path over floored ground accumulates 0.1, against a typical
+#' `cost_threshold` of 2500.
+#'
+#' Negative friction is not floored — [terra::costDist()] rejects a negative
+#' cost surface, and that error is left intact.
+#'
+#' If your friction is in units whose typical values approach `1e-6`, floor the
+#' raster yourself before calling.
+#'
 #' @examples
 #' dem <- terra::rast(system.file("testdata/dem.tif", package = "flooded"))
 #' slope <- terra::rast(system.file("testdata/slope.tif", package = "flooded"))
@@ -43,11 +58,24 @@ fl_cost_distance <- function(friction, streams) {
          call. = FALSE)
   }
 
-  # costDist(x, target) finds cells in x equal to `target` as seed points.
-
-  # Set stream cells to 0 in the friction raster so costDist treats them as
+  # costDist(x, target) seeds on every cell in x equal to `target`, so encoding
+  # stream cells as 0 is only correct if nothing else is 0. Friction rasters do
+  # contain exact zeros — integer-metre DEMs, hydro-flattened lake surfaces and
+  # void-filled plateaus all quantize to perfectly flat — and each one was
+  # silently acting as a free cost source (#41).
+  #
+  # Floor them to a negligible positive value so that zero means "stream cell"
+  # by construction. Percent slope runs 0-100+, so 1e-6 is eight orders below
+  # the signal: costDist accumulates friction x distance in map units, making a
+  # 100 km path over floored ground worth 0.1 against a default cost_threshold
+  # of 2500. Flat ground stays cheap to cross; it just stops being a source.
+  #
+  # `== 0` and not `<= 0`: terra rejects a negative cost surface outright, and
+  # flooring negatives would disable that guard, turning meaningless input into
+  # plausible-looking output. NA is left alone (NA == 0 is NA), so impassable
+  # cells stay impassable.
+  friction <- terra::ifel(friction == 0, 1e-6, friction)
 
-  # target cells (cost = 0 starting points).
   cost <- terra::ifel(!is.na(streams), 0, friction)
   out <- terra::costDist(cost, target = 0)
   names(out) <- "cost_distance"

R/fl_valley_poly.R

@@ -29,7 +29,7 @@ fl_valley_poly <- function(x, dissolve = TRUE) {
   x_mask <- terra::ifel(x == 1, 1L, NA)
 
   # Polygonize
- polys <- terra::as.polygons(x_mask, dissolve = dissolve)
+  polys <- terra::as.polygons(x_mask, dissolve = dissolve)
 
   # Name the column on the SpatVector, before it becomes an sf. Renaming by
   # position afterwards hits the geometry column when as.polygons() returns no

inst/notes/methodology.md

@@ -117,10 +117,48 @@ intersection, and waterbodies get no spatial filter, so those cells can satisfy
 group; `attr(x, "fl_fallback_cells")` reports how many. An unusually large count is also the signal
 that `max_width` / `cost_threshold` do not match the values the delineation was built with.
 
+### Only stream cells seed the cost surface
+
+`fl_cost_distance()` encodes seeds by setting stream cells to zero and calling
+`terra::costDist(target = 0)`, which seeds on *every* zero cell. Until 0.4.1 that included any cell
+whose friction was already exactly zero, so flat ground acted as a free cost source. Fixed by
+flooring `friction == 0` to `1e-6` before seeding (flooded#41).
+
+Flat ground stays cheap to *cross* — the floor accumulates 0.1 over a 100 km path at 10 m against a
+default `cost_threshold` of 2500 — it just stops being a *source*. Negative friction is deliberately
+not floored, so `costDist()`'s own rejection of a negative cost surface is left intact.
+
+Which DEMs contain exact zeros, and what changes when they do — measured on both datasets this
+package ships:
+
+| DEM | exact-zero slope cells | cost mask (`< 2500`) | valley cells |
+|---|---|---|---|
+| bundled `dem.tif` / `slope.tif`, 10 m | 0 of 45,726 (min 1.42e-14) | unchanged | 53,635 -> 53,635 |
+| `pars_dem.tif` (MRDEM-30, 30 m, 20.9 Mcell) | 80 of 10.7 M | -2,289 cells (214 ha), 0 added | 521,028 -> 521,028 |
+
+Two things worth separating. MRDEM-30 — the package default source — **does** produce exact zeros
+at watershed scale, so a small clip returning none is not evidence about the source. And the fix
+only ever *removes* cells from the cost mask, never adds, because it can only raise a cost that was
+spuriously zero.
+
+But a change in the cost mask is not a change in the delineation. Both shipped datasets come out
+bit-identical, because `fl_valley_confine()` intersects cost with slope, distance and flood and then
+runs morphological cleanup — enough to absorb all 2,289 Parsnip cells. That is a property of these
+two datasets, not a guarantee: wherever cost is the binding criterion (flat terrain, a lax
+`slope_threshold`, a large `flood_factor`) the delineation will move.
+
+The 1 m lidar run in `vignettes/stac-dem.Rmd` emits `[costDist] distance algorithm did not
+converge`, which is the shape of large zero-cost plateaus — suggestive, not confirmed, since that
+vignette is pre-baked and was not re-run.
+
+It matters most under attribution. Cost is what separates one watercourse's floodplain from
+another's, so a flat patch inside a group's corridor would spread that group's mask across ground
+its own streams never reach — cost failing to discriminate exactly where floodplains are, on flat
+ground.
+
 ### See also
 
 - `fl_valley_attribute()` docs — corridor cropping (`crop_margin`) is an approximation, not a bound
-- flooded#41 — `fl_cost_distance()` seeds every zero-friction cell; matters more per-group
 - flooded#44 — production-scale timing is unmeasured; the k=5 figure does not establish k=340
 - NewGraphEnvironment/floodplains#40 — driver-side half (config surface, key column on the gpkg)
 

planning/active/review-round1.md

@@ -0,0 +1,123 @@
+# Review round 1 — #41 zero-friction seeding fix
+
+Reviewed: staged diff (`R/fl_cost_distance.R`, `R/fl_valley_poly.R`,
+`tests/testthat/test-fl_cost_distance.R`, plus `man/`, `inst/notes/methodology.md`,
+`planning/` in the full staged set). Every claim below was probed, not reasoned.
+terra 1.9.34.
+
+The fix itself is correct. All findings are about artifacts and claims the fix
+orphans, plus one guard weaker than its comment says.
+
+## Findings
+
+- **[bug]** `inst/vignette-data/pars_valleys.tif` and the `floodplain` layer in
+  `inst/vignette-data/pars.gpkg` are cached outputs of the **buggy** code, and the
+  Parsnip DEM they came from does contain exact zeros — so the published
+  `vignettes/pars-floodplain.Rmd` now shows a floodplain the package no longer
+  produces. Measured on `inst/vignette-data/pars_dem.tif`
+  (`tan(terrain(dem,"slope",unit="degrees")*pi/180)*100`, the exact `slope = NULL`
+  path `fl_valley_confine()` takes):
+
+  ```
+  non-NA cells   10,688,931
+  exact zeros            80      <- 80 free cost sources under the old code
+  min nonzero     2.58e-15
+  ```
+
+  Re-running cost distance over that DEM + `pars.gpkg` streams, old vs new:
+
+  ```
+  cost-mask (cost < 2500) cells differing:  2,289 of 10,688,931
+  old-only (spurious under the bug):        2,289   (~214 ha at 30 m)
+  ```
+
+  The `vca` chunk at `vignettes/pars-floodplain.Rmd:218` is `eval = FALSE`, so the
+  vignette renders the committed raster/gpkg rather than recomputing — the stale
+  numbers publish to pkgdown and, per the appendix-port convention in `CLAUDE.md`,
+  flow into the fp_peace reporting appendix. Nothing in the diff flags this.
+  Either regenerate the cached bundle or state in the vignette which version
+  produced it. (Secondary, already acknowledged in `inst/notes/methodology.md:142`:
+  `vignettes/stac-dem.Rmd:311` is pre-baked and carries a
+  `[costDist] distance algorithm did not converge` warning plus a
+  "3,883,179 valley cells / 27.7 %" figure from a 1 m lidar run under the bug.)
+
+- **[bug]** `NEWS.md` (0.4.1 bullet 2) and `inst/notes/methodology.md:136-139`
+  generalize a one-AOI measurement into a claim about the package's default DEM
+  source: *"the bundled tile and a 30 m MRDEM-30 clip over the same AOI both
+  contain zero exact-zero slope cells"* / *"neither the bundled tile nor the
+  package's default DEM source is affected here"*. The package **ships** a second
+  MRDEM-30 clip — Parsnip — with 80 exact-zero cells and a 2,289-cell mask delta
+  (above). A reader on MRDEM-30 will read this as "I am not affected" and be
+  wrong. Scope the claim to the Bulkley AOI it was measured on, and add the
+  Parsnip row to the methodology table.
+
+- **[fragile]** `tests/testthat/test-fl_cost_distance.R:107-119` — the guard is
+  weaker than its comment. It claims to catch *"flooring to 1, or to the median
+  friction"*. Measured by patching the floor and re-running the file (namespace
+  **and** globalenv, see note below):
+
+  | floor  | `through_flat` | `over_slope` | line 118 |
+  |--------|----------------|--------------|----------|
+  | 1e-6   | 4101.22        | 4596.19      | pass (correct) |
+  | 1      | 4150.72        | 4596.19      | **pass** — over-correction undetected |
+  | 10     | 4596.19        | 4596.19      | fail ✓ |
+  | 100    | 7621.93        | 4596.19      | fail ✓ |
+
+  The 6×6 flat patch is 60 m across against a ~4.1 km path, so a floor of 1 barely
+  moves the total. Floor 1 is exactly the over-correction the comment describes:
+  1 % slope × 100 km = 100,000, forty times `cost_threshold`, i.e. flat ground more
+  than ~2.5 km from a stream would fail the cost mask. If the guard is meant to
+  hold, make the flat patch large enough to dominate the path, or assert an
+  absolute bound on the flat-crossing increment rather than a relative one.
+
+- **[fragile]** `tests/testthat/test-fl_cost_distance.R:133` — `expect_error(...,
+  "negative friction")` matches terra's own wording, not this package's;
+  `fl_cost_distance()` raises nothing here. Verified message is
+  `[costDist] negative friction values not allowed` (terra 1.9.34). A terra
+  rewording turns this green guard red and reads as a regression in `flooded`.
+  Low priority — worth a comment naming terra as the source of the string.
+
+## Verified — not findings
+
+Recorded so they are not re-litigated.
+
+- **Restore-the-bug.** Reverting `R/fl_cost_distance.R:77` and re-running the file:
+  `FAIL 3 | PASS 18`. Genuine detectors are lines **95** (zero cells ≠ seed cells:
+  37 vs 1), **104** (`0.0 <= 50.0`), **156** (17 zeros vs 1). The other four new
+  tests pass in both states and are labelled in-file as guards, which is accurate.
+  *Method note:* `testthat::test_file()` here resolves `fl_cost_distance` through
+  globalenv, not `asNamespace("flooded")` — patching only the namespace gives a
+  false green. Both bindings must be assigned.
+- **Checklist 5, `== 0` vs tolerance.** `terra::costDist(target=)` is exact
+  equality against the same doubles R compares, so nothing can seed without
+  satisfying `friction == 0`. The one theoretical gap — a nonzero double below
+  float32 normal min (~1.18e-38) flushing to 0 when the `ifel` intermediate spills
+  to FLT4S — is unreachable for percent slope: bundled min 1.42e-14, Parsnip min
+  2.58e-15. `-0.0 == 0` is TRUE, so signed zero is floored. NaN condition is NA →
+  cell stays NA.
+- **Checklist 6 + 9, type promotion.** `terra::ifel(friction == 0, 1e-6, friction)`
+  on an INT2S **file-backed** source returns FLT4S both in memory and under
+  `terraOptions(todisk=TRUE)` (default write datatype is FLT4S regardless of source
+  type) — verified; the floor is not rounded away. 1e-6 in float32 is
+  9.99999997e-07, not 0. `NA == 0` is NA and `ifel` keeps NA, confirmed by the
+  barrier test at line 159.
+- **Checklist 7, aliasing.** `friction <- terra::ifel(friction == 0, ...)`
+  evaluates the argument before rebinding; no recursion or aliasing. Ordering is
+  right — the floor runs *before* seed encoding; reversed it would erase the seeds.
+- **Convergence under the floor.** A 1500×1500 grid that is 87 % floored produces
+  no `did not converge` warning, exactly one zero cell, and a sane range. The
+  Parsnip run (10.7 M cells) also emitted none. Flooring does not make push-broom
+  convergence worse.
+- **Checklist 8, callers + inventory.** `terra::costDist` has exactly one call site
+  in `R/`. `fl_valley_confine():147` and `fl_group_cells()` in
+  `fl_valley_attribute.R:315` both pass single-layer percent slope (cropped in the
+  attribute path) and a rasterized stream seed layer — nothing they pass breaks
+  under the new behaviour. The bundled-data equality test at line 175 covers both
+  routes, and its premise assertions (lines 182, 188) are real.
+- **Full suite:** `FAIL 0 | WARN 0 | SKIP 0 | PASS 246`.
+- **`R/fl_valley_poly.R:32`** is a one-space indent fix, no behaviour change.
+- **Bookkeeping, not a defect:** `DESCRIPTION` (0.4.0 → 0.4.1) and `NEWS.md` are
+  modified but **unstaged**. Committing the staged set as-is lands the fix without
+  the version bump and NEWS entry. That matches this repo's "version bump as the
+  final commit" convention — noting only so it is deliberate rather than an
+  oversight.

planning/archive/2026-08-issue-41-fl-cost-distance-zero-friction-seeds/README.md

@@ -0,0 +1,15 @@
+# Issue #41 — fl_cost_distance() seeds every zero-friction cell
+
+`fl_cost_distance()` documented stream cells as its seed points and did not deliver that: seeds are
+encoded by setting stream cells to `0` and calling `terra::costDist(target = 0)`, which matches
+*every* zero cell, so any cell whose friction was already exactly zero acted as a free cost source.
+Fixed by flooring `friction == 0` to `1e-6` before seeding. Two things changed from the plan under
+measurement — the issue's proposed `<= 0` floor would have disabled terra's own rejection of a
+negative cost surface, and a first over-correction guard turned out too weak to catch a floor of 1.
+
+The fix strictly removes cells from the cost mask and never adds. On both DEMs this package ships
+the delineation is unmoved (bundled 53,635 cells; Parsnip 521,028, with 2,289 cost-mask cells
+absorbed by the other three criteria and cleanup) — so no shipped artifact needed regenerating. That
+is a property of these datasets, not a guarantee.
+
+Closed by PR #45, released as v0.4.1.