diff --git a/DESCRIPTION b/DESCRIPTION index e2acf1a..f77daa2 100644 --- a/DESCRIPTION +++ b/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", , "al@newgraphenvironment.com", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-3495-2128")), diff --git a/NEWS.md b/NEWS.md index 32be3b2..6f1ad13 100644 --- a/NEWS.md +++ b/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. diff --git a/R/fl_cost_distance.R b/R/fl_cost_distance.R index 6430320..7b7755f 100644 --- a/R/fl_cost_distance.R +++ b/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" diff --git a/R/fl_valley_poly.R b/R/fl_valley_poly.R index d85caf8..5b3ce05 100644 --- a/R/fl_valley_poly.R +++ b/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 diff --git a/inst/notes/methodology.md b/inst/notes/methodology.md index 59062e5..d21a3b3 100644 --- a/inst/notes/methodology.md +++ b/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) diff --git a/man/fl_cost_distance.Rd b/man/fl_cost_distance.Rd index 6f36aae..94c4bb9 100644 --- a/man/fl_cost_distance.Rd +++ b/man/fl_cost_distance.Rd @@ -29,6 +29,21 @@ weighted distance. The \code{friction} raster defines per-cell traversal cost and \code{streams} identifies seed cells (cost = 0). Cells that are \code{NA} in \code{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 \code{0} are floored to \code{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 +\code{cost_threshold} of 2500. + +Negative friction is not floored — \code{\link[terra:costDist]{terra::costDist()}} rejects a negative +cost surface, and that error is left intact. + +If your friction is in units whose typical values approach \code{1e-6}, floor the +raster yourself before calling. } \examples{ dem <- terra::rast(system.file("testdata/dem.tif", package = "flooded")) diff --git a/planning/active/review-round1.md b/planning/active/review-round1.md new file mode 100644 index 0000000..6b93ea0 --- /dev/null +++ b/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. diff --git a/planning/archive/2026-08-issue-41-fl-cost-distance-zero-friction-seeds/README.md b/planning/archive/2026-08-issue-41-fl-cost-distance-zero-friction-seeds/README.md new file mode 100644 index 0000000..31dea4b --- /dev/null +++ b/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. diff --git a/planning/archive/2026-08-issue-41-fl-cost-distance-zero-friction-seeds/findings.md b/planning/archive/2026-08-issue-41-fl-cost-distance-zero-friction-seeds/findings.md new file mode 100644 index 0000000..1c0d5b9 --- /dev/null +++ b/planning/archive/2026-08-issue-41-fl-cost-distance-zero-friction-seeds/findings.md @@ -0,0 +1,194 @@ +# Findings — fl_cost_distance() seeds every zero-friction cell (#41) + +## Issue context + +`fl_cost_distance()` treats every zero-friction cell as a stream seed, not just stream cells. + +```r +# R/fl_cost_distance.R:51-52 +cost <- terra::ifel(!is.na(streams), 0, friction) +out <- terra::costDist(cost, target = 0) +``` + +`costDist(target = 0)` finds *all* cells equal to 0. Stream cells are set to 0 deliberately — but so +is any cell whose friction (percent slope) is already exactly 0. The roxygen says otherwise. + +Filed separately from #40 because it changes VCA output on any DEM containing exact zeros — a +result-changing decision that deserves its own diff and regression check, not a side effect of an +attribution feature. Found during the plan review for #40. + +## Reproduction (measured 2026-08-27) + +Synthetic 50x50 grid, friction 10% everywhere, one stream cell at [45,45], a 6x6 exact-zero patch at +rows/cols 10-15: + +``` + current fixed +zero-cost cells 37 1 (stream cells: 1; 36 patch cells were seeds) +cost at flat-patch centre 0 4101 +cost one cell beyond it 50 4151 +``` + +## Why the bundled tile cannot catch it + +| raster | min | cells == 0 | +|---|---|---| +| `inst/testdata/slope.tif` | 1.422896e-14 | 0 | +| slope derived from `dem.tif` by the `slope = NULL` path in `fl_valley_confine()` | 1.422896e-14 | 0 | + +Both verified this session. `terrain()`-derived slope on a continuous float DEM essentially never +produces an exact zero. Production DEMs do: integer-metre DEMs, hydro-flattened lake surfaces, +void-filled plateaus. + +Corroboration already in the repo: the 1 m lidar run at `vignettes/stac-dem.Rmd:311` emits +`[costDist] distance algorithm did not converge`, which is the shape of large zero-cost plateaus. + +## Design: why flooring, not a sentinel target + +`terra::costDist(x, target)` takes a target *value* in `x`; there is no separate seed-raster +argument. So seeds must be encoded as a value in the friction raster. Two ways: + +**A — sentinel `target = -1`.** Leaves genuinely flat ground free to cross. +**B — floor non-positive friction to `EPS`, keep `target = 0`.** + +Probed both on the grid above; numerically identical (patch centre 4101 either way, exactly 1 zero +cell). **B chosen.** A *assumes* no friction cell equals −1; B *constructs* the invariant that +nothing but a stream cell equals 0. A is the same bug relocated to a rarer value. + +## Properties verified by probe, not assumed + +| property | probe result | +|---|---| +| `NA` friction survives the floor (barriers preserved) | `ifel(NA <= 0, ...)` -> `NA`. Confirmed | +| integer-typed friction does not round the floor back to 0 | `INT2S` raster promotes to float under `ifel`; min after floor 1e-06, zeros 0 | +| flat ground stays cheap to *cross*, not made a barrier | 4101 through the flat patch vs 4596 over equal-length 10 %-friction ground | +| a stream cell whose friction is `NA` still seeds at 0 | confirmed — seeding runs after flooring | +| all-flat friction raster | exactly 1 zero cell (the stream); corner cost 1.2e-04 | + +## Scale of `EPS = 1e-6` + +`costDist` accumulates friction × distance in map units — confirmed empirically: a 10x10 grid at +10 m res, friction 1e-6, corner-to-corner ≈ 127 m gave 1.202e-04 ≈ 127 × 1e-6. + +So at 10 m cells a **100 km** path over floored ground accrues **0.1** against a default +`cost_threshold` of 2500. Percent slope runs 0–100+; `EPS` is eight orders below. Negligible by +construction, and the probe above shows it does not read as a barrier. + +## No new argument + +`EPS` stays an internal constant. `fl_cost_distance()` is exported, so a user in unusual friction +units already has a one-line escape hatch — pre-floor their own raster before calling. A parameter +buys nothing they cannot already do. This keeps the change a pure bug fix (0.4.0 -> 0.4.1). + +## Relationship to #40 + +`fl_cost_distance()` is the package's only `costDist` call site — `R/fl_valley_confine.R:147` and +`R/fl_valley_attribute.R:315`. One fix covers both. + +Under `fl_valley_confine()` alone a spurious flat seed lowers cost slightly in ground that was +probably valley anyway. Under per-group attribution it is worse: a flat patch is a free source for +whichever group's corridor contains it, so that group's cost mask spreads across ground its own +streams never reach. Cost then stops discriminating between groups exactly where floodplains are — +on flat ground — and membership degenerates toward the distance buffer. + +#40's coverage argument ("a valley cell's global cost is reproduced by its own group's seeds") is +*strengthened* by this fix, not perturbed. + +## Revised during implementation: floor `== 0`, not `<= 0` + +The issue proposed `terra::ifel(friction <= 0, 1e-6, friction)`. Probing found that +`terra::costDist()` **rejects a negative cost surface outright**: + +``` +Error: [costDist] negative friction values not allowed +``` + +So flooring `<= 0` would have silently disabled a real guard, converting meaningless input into +plausible-looking output — the same failure direction the rest of this fix exists to close. Landed +as `friction == 0` instead, which eliminates exactly the set that `costDist(target = 0)` would +mistake for seeds and nothing more. Verified the guard still fires after the change. + +`NA == 0` is `NA`, so barriers are still untouched. + +## Does the default DEM source hit this? Yes — corrected after review + +My first measurement was one 30 m MRDEM-30 clip over the Bulkley test AOI: **0** exact-zero slope +cells. I wrote that into NEWS and `methodology.md` as "the default source is unaffected". A +code-review agent caught that the package's own *other* MRDEM-30 dataset contradicts it. Verified: + +| DEM shipped by this package | 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 | +| MRDEM-30 clip over the same AOI, 30 m | 0 of 45,726 (min 0.0041) | unchanged | n/a | +| `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 | + +Textbook "an inventory is only complete relative to a boundary — name the boundary". The +measurement was correct and the generalisation from it was not. A small clip returning zero is not +evidence about the source; scale is what surfaces the zeros. + +### But the reviewer's conclusion from it was wrong in the other direction + +It reported the shipped `pars_valleys.tif` as **stale** — produced under the bug, published by +pkgdown, feeding the fp_peace report. That inference came from the cost mask, an *intermediate*. +Measuring the actual output instead — the shipped raster **is** the old code's output, so it is its +own oracle — gives: + +``` +shipped(old): 521028 cells (48603.1 ha) +current(new): 521028 cells (48603.1 ha) +only in old : 0 only in new : 0 +``` + +Bit-identical. `fl_valley_confine()` intersects cost with slope, distance and flood and then runs +morphological cleanup, which absorbs all 2,289 cells. No shipped artifact needs regenerating — +which matters, because `data-raw/wsg_vignette_data.R` needs an fwapg database connection that is not +available here. + +Worth keeping as the general lesson: **a moved intermediate is not a moved output.** Both directions +of this finding cost a measurement to settle, and both were worth making. + +One property that does generalise: the fix can only *raise* a cost that was spuriously zero, so it +strictly removes cells from the cost mask and never adds (0 added, confirmed). Where cost is the +binding criterion the delineation will shrink; it can never grow. + +## Bug-restoration check + +Tests were run against the unfixed function before Phase 3. Three went red — and only the three +that target the bug: + +``` +test-fl_cost_distance.R:95 zero_cells != seed_cells (37 zeros vs 1 expected) +test-fl_cost_distance.R:104 patch centre 0.0 <= near-stream 50.0 +test-fl_cost_distance.R:156 INT2S output had 17 zeros, expected 1 +[ FAIL 3 | WARN 0 | SKIP 0 | PASS 18 ] +``` + +The other new tests pass in both states by design — they guard the *opposite* over-correction +(flooring high enough to make flat ground a barrier) and terra's negative-friction guard, neither of +which the bug touches. + +## Review round 1 — what landed + +| finding | verdict | action | +|---|---|---| +| shipped `pars_valleys.tif` is stale | **wrong** — cost mask moved, output did not (0 cells) | none; recorded above | +| NEWS/methodology overclaim the default source is unaffected | **right** | rewritten with the two-DEM table | +| over-correction guard cannot catch a floor of 1 | **right** | replaced with a negligibility ratio guard | +| negative-friction test pins terra's error string | **right**, minor | loosened to `"negative"` with a comment | + +The third is the one worth remembering. The test asserted "crossing flat ground costs less than +crossing sloped ground", which *any* floor below the sloped friction satisfies — so a floor of 1 +passed, while costing 1e5 over a 100 km path, 40x a default `cost_threshold`. The assertion was +correct and the property it encoded was too weak. Replaced with a ratio bound (flat traverse under +3e-5 of the sloped equivalent), then verified by restoration: floors of 1e-6 and 1e-4 pass, 1e-3 and +1 fail, with the patch confirmed to have taken effect on each run. + +Also noted by the reviewer and worth carrying: under `pkgload::load_all()` / `test_file()`, +`fl_cost_distance` resolves through `globalenv()` as well as the namespace, so a restoration probe +that patches only `asNamespace("flooded")` gives a **false green**. Patch both, and print a value +that proves the patch took before trusting the run. + +## Errors Encountered + +| Error | Resolution | +|-------|------------| diff --git a/planning/archive/2026-08-issue-41-fl-cost-distance-zero-friction-seeds/progress.md b/planning/archive/2026-08-issue-41-fl-cost-distance-zero-friction-seeds/progress.md new file mode 100644 index 0000000..5325987 --- /dev/null +++ b/planning/archive/2026-08-issue-41-fl-cost-distance-zero-friction-seeds/progress.md @@ -0,0 +1,24 @@ +# Progress — fl_cost_distance() seeds every zero-friction cell (#41) + +## Session 2026-08-27 + +- Plan-mode exploration — reproduced the bug on the issue's synthetic grid, probed and rejected the + sentinel-target alternative, verified NA / integer / flat-crossing / all-flat edge cases +- Phases approved by user +- Created branch `41-fl-cost-distance-seeds-every-zero-friction-c` off main +- Scaffolded PWF baseline with approved phases +- Phase 2 — 8 new tests in `test-fl_cost_distance.R`, all synthetic (the bundled tile cannot reach + this failure mode). Confirmed 3 red against the unfixed function before writing the fix +- Phase 2 revision — probing showed `costDist()` rejects negative friction, so the issue's proposed + `<= 0` floor would have disabled that guard. Landed `== 0` instead; task_plan corrected +- Phase 3 — floored `friction == 0` to `1e-6` before seeding; roxygen `@details` rewritten. Suite + 246 pass / 0 fail +- Phase 4 — measured MRDEM-30 over the test AOI: 0 exact-zero slope cells, so the package default + source is unaffected here. Bundled results confirmed unmoved (53,635 valley cells; attribution + unchanged). Recorded in `inst/notes/methodology.md` +- Phase 5 — NEWS + version bump 0.4.1, PR #45 opened +- Review round 1 (concurrent, landed post-PR) — 4 findings. Two real defects fixed on the branch + (weak over-correction guard, NEWS/methodology overclaiming the default DEM source); one disproved + by measuring the output rather than the intermediate (shipped `pars_valleys.tif` is bit-identical, + not stale); one minor test-fragility fix +- Next: archive PWF diff --git a/planning/archive/2026-08-issue-41-fl-cost-distance-zero-friction-seeds/task_plan.md b/planning/archive/2026-08-issue-41-fl-cost-distance-zero-friction-seeds/task_plan.md new file mode 100644 index 0000000..38179cf --- /dev/null +++ b/planning/archive/2026-08-issue-41-fl-cost-distance-zero-friction-seeds/task_plan.md @@ -0,0 +1,67 @@ +# Task: fl_cost_distance() seeds every zero-friction cell, not just stream cells (#41) + +## Problem + +`fl_cost_distance()` documents stream cells as the seed points, and does not deliver that. +`R/fl_cost_distance.R:51-52` sets stream cells to 0 and hands the raster to +`terra::costDist(target = 0)`, which finds **every** cell equal to 0 — including any cell whose +friction (percent slope) is already exactly zero. Those cells become free cost sources they were +never meant to be. + +Fix: floor non-positive friction to a negligible positive value *before* seeding, so zero is +constructed to mean "stream cell" and nothing else. + +## Phase 1: Branch + PWF baseline + +- [x] Branch `41-fl-cost-distance-seeds-every-zero-friction-c` off main +- [x] `planning/active/` — `task_plan.md`, `findings.md`, `progress.md` + +## Phase 2: Tests first — must be red before Phase 3 + +- [x] **Only stream cells are seeds** — `which(cost == 0)` is exactly the stream cell indices +- [x] **A flat patch is no longer a cost sink** — patch-centre cost exceeds near-stream cost +- [x] **The floor does not turn flat ground into a barrier** — crossing flat ground costs less than + an equal-length path over 10 %-friction ground +- [x] **Negative friction still errors** — revised during Phase 2. `costDist()` rejects a negative + cost surface, so flooring `<= 0` as the issue proposed would have disabled a real guard. + Floor `== 0` only; the test now asserts the guard survives +- [x] **Integer-typed friction survives the floor** (`INT2S`, no zeros after flooring) +- [x] **NA friction stays a barrier** +- [x] **Bundled data is unchanged** — premise (`sum(slope <= 0) == 0`) asserted in the same test +- [x] **The fix cannot move any bundled-data result** — revised: assert cost-raster equality plus + the no-exact-zeros premise, which implies the `fl_valley_confine()` claim without a second + full VCA run. Confirmed empirically too: 53,635 valley cells, unchanged +- [x] Restore the bug via namespace patch, confirm the core tests go red, record in `progress.md` + +## Phase 3: Fix + documentation + +- [x] `R/fl_cost_distance.R` — floor `friction == 0` before seeding, scale reasoning in a comment +- [x] Roxygen — `@details` on flooring, NA-as-barrier, and the negative-friction guard +- [x] `devtools::document()` + +## Phase 4: Verify + measure + +- [x] `devtools::test()` 246 pass / 0 fail; `check()` 0 errors 0 warnings 1 pre-existing NOTE + (`pkgdown/`); `lint_package()` leaves only two pre-existing lints outside this diff +- [x] MRDEM-30: **0** exact zeros on the small test-AOI clip but **80** on the shipped 20.9 Mcell + `pars_dem.tif`. First reading generalized from one clip and was corrected after review — the + default source *does* produce exact zeros at scale. Cost mask moves 2,289 cells (214 ha, + 0 added); the delineation does not move at all (521,028 cells both ways) +- [x] `inst/notes/methodology.md` — forward reference replaced with the finding and the DEM table +- [x] Confirmed unmoved: 53,635 valley cells, same 5 attribution rows, 0 fallback cells. Also + verified `costDist` target matching is exact equality (a 1e-14 cell reads 7.07e-14, not 0) + and that the zero set on bundled data is now identical to the 1,607 stream cells + +## Phase 5: Release + close + +- [x] `/code-check` — round 1 found 2 real defects (fixed on branch), 1 disproved by measurement, + 1 minor fragility (fixed) +- [x] `NEWS.md` + `DESCRIPTION` `0.4.0` → `0.4.1` as the final commit +- [x] PR #45 opened; `/planning-archive` on merge + +## Validation + +- [x] Tests pass +- [x] `/code-check` run; findings folded in +- [x] PWF checkboxes match landed work +- [ ] `/planning-archive` on completion diff --git a/tests/testthat/test-fl_cost_distance.R b/tests/testthat/test-fl_cost_distance.R index 62bff8a..cbec265 100644 --- a/tests/testthat/test-fl_cost_distance.R +++ b/tests/testthat/test-fl_cost_distance.R @@ -60,3 +60,157 @@ test_that("fl_cost_distance errors on mismatched grids", { expect_error(fl_cost_distance(r1, r2), "same extent") }) + +# --- Zero-friction cells must not seed (#41) ------------------------------- +# +# The bundled tile cannot reach this failure mode: slope.tif has no cell equal +# to zero (min 1.42e-14), and neither does slope derived from dem.tif. These +# fixtures are synthetic for that reason, not for convenience. + +# 50x50 at 10 m, friction 10% everywhere, one stream cell at [45, 45]. +# `flat` adds a 6x6 patch of exactly-zero friction at rows/cols 10-15, far from +# the stream — the shape a hydro-flattened lake or a void-filled plateau takes. +zero_friction_grid <- function(flat = TRUE) { + friction <- terra::rast(nrows = 50, ncols = 50, vals = 10, + xmin = 0, xmax = 500, ymin = 0, ymax = 500, + crs = "EPSG:3005") + if (flat) friction[10:15, 10:15] <- 0 + streams <- terra::rast(friction) + terra::values(streams) <- NA + streams[45, 45] <- 1 + list(friction = friction, streams = streams) +} + +test_that("only stream cells are cost-distance seeds", { + g <- zero_friction_grid() + + # Premise: the fixture must actually contain exact-zero friction, or the + # assertion below passes for nothing. + expect_gt(sum(terra::values(g$friction, mat = FALSE) == 0, na.rm = TRUE), 0) + + cd <- fl_cost_distance(g$friction, g$streams) + + zero_cells <- which(terra::values(cd, mat = FALSE) == 0) + seed_cells <- which(!is.na(terra::values(g$streams, mat = FALSE))) + expect_equal(zero_cells, seed_cells) +}) + +test_that("a flat patch is not a cost sink", { + g <- zero_friction_grid() + cd <- fl_cost_distance(g$friction, g$streams) + + # Under the bug the patch centre reads 0 — cheaper than ground adjacent to + # the stream itself, which is the tell. + expect_gt(cd[12, 12][[1]], cd[45, 44][[1]]) +}) + +test_that("the floor is negligible, not merely cheaper than sloped ground", { + # Guards the opposite over-correction: a floor high enough to make flat + # ground effectively impassable at scale. + # + # The obvious framing — "crossing the flat patch costs less than crossing + # sloped ground" — is far too weak to catch it, because any floor below the + # sloped friction satisfies it. Measured: a floor of 1 passes that version, + # while costing 1e5 over a 100 km path, or 40x a default cost_threshold of + # 2500. So assert negligibility against the sloped case instead. + # + # The ratio is the floor divided by the sloped friction, so 3e-5 passes any + # floor up to 1e-4 (300x headroom over the 1e-6 in use) and rejects 1e-3 and + # above, which is where a floor starts to consume the cost budget. + flat <- terra::rast(nrows = 50, ncols = 50, vals = 0, + xmin = 0, xmax = 500, ymin = 0, ymax = 500, + crs = "EPSG:3005") + sloped <- terra::setValues(flat, 10) + streams <- terra::rast(flat) + terra::values(streams) <- NA + streams[1, 1] <- 1 + + across_flat <- fl_cost_distance(flat, streams)[50, 50][[1]] + across_slope <- fl_cost_distance(sloped, streams)[50, 50][[1]] + + expect_gt(across_slope, 0) # premise: the reference is real + expect_lt(across_flat / across_slope, 3e-5) +}) + +test_that("negative friction still errors", { + # terra rejects a negative cost surface. Flooring `<= 0` rather than `== 0` + # would silently disable that guard and turn meaningless input into + # plausible-looking output. + friction <- terra::rast(nrows = 20, ncols = 20, vals = 10, + xmin = 0, xmax = 200, ymin = 0, ymax = 200, + crs = "EPSG:3005") + friction[5, 5] <- -3 + streams <- terra::rast(friction) + terra::values(streams) <- NA + streams[18, 18] <- 1 + + # The message is terra's, not this package's — matching loosely on "negative" + # so a reworded upstream error does not read as the guard having gone. + expect_error(fl_cost_distance(friction, streams), "negative") +}) + +test_that("integer-typed friction survives the floor", { + # An integer raster must promote to float, or the floor rounds back to zero + # and the fix silently does nothing. + friction <- terra::rast(nrows = 20, ncols = 20, + xmin = 0, xmax = 200, ymin = 0, ymax = 200, + crs = "EPSG:3005") + terra::values(friction) <- rep(10L, 400) + friction[5:8, 5:8] <- 0L + + tf <- tempfile(fileext = ".tif") + on.exit(unlink(tf), add = TRUE) + terra::writeRaster(friction, tf, datatype = "INT2S", overwrite = TRUE) + friction <- terra::rast(tf) + expect_equal(terra::datatype(friction), "INT2S") + + streams <- terra::rast(friction) + terra::values(streams) <- NA + streams[18, 18] <- 1 + + cd <- fl_cost_distance(friction, streams) + expect_equal(sum(terra::values(cd, mat = FALSE) == 0, na.rm = TRUE), 1L) +}) + +test_that("NA friction remains an impassable barrier", { + friction <- terra::rast(nrows = 20, ncols = 20, vals = 10, + xmin = 0, xmax = 200, ymin = 0, ymax = 200, + crs = "EPSG:3005") + friction[, 10] <- NA # a wall the full height of the grid + streams <- terra::rast(friction) + terra::values(streams) <- NA + streams[1, 1] <- 1 + + cd <- fl_cost_distance(friction, streams) + + expect_true(is.na(cd[10, 15][[1]])) # sealed off behind the wall + expect_true(is.na(cd[10, 10][[1]])) # the wall itself + expect_false(is.na(cd[10, 5][[1]])) # same side as the stream +}) + +test_that("the fix cannot move any bundled-data result", { + slope <- terra::rast(testdata_path("slope.tif")) + dem <- terra::rast(testdata_path("dem.tif")) + + # Premise, asserted here so a future test-data swap fails on this line — + # naming the real cause — rather than somewhere downstream. Every + # bundled-data baseline in this suite rests on it. + expect_equal(sum(terra::values(slope, mat = FALSE) == 0, na.rm = TRUE), 0L) + + # The slope = NULL path in fl_valley_confine() derives slope rather than + # reading it, so it needs the premise checked too. + slope_deg <- terra::terrain(dem, "slope", unit = "degrees") + derived <- tan(slope_deg * pi / 180) * 100 + expect_equal(sum(terra::values(derived, mat = FALSE) == 0, na.rm = TRUE), 0L) + + # Cost is the only route by which this change reaches fl_valley_confine() or + # fl_valley_attribute(), so equality here covers both. + streams_sf <- sf::st_read(testdata_path("streams.gpkg"), quiet = TRUE) + stream_r <- fl_stream_rasterize(streams_sf, dem, field = "channel_width") + + unfloored <- terra::costDist(terra::ifel(!is.na(stream_r), 0, slope), target = 0) + expect_equal( + as.vector(terra::values(fl_cost_distance(slope, stream_r))), + as.vector(terra::values(unfloored)) + ) +})