Report waterbody contribution to floodplain area for LULC analysis interpretation

[NewGraphEnvironment]

Problem

Waterbodies passed to fl_valley_confine() are stamped onto the output AFTER the VCA runs (line 218: ifel(wb_r == 1L, 1L, valleys)). They bypass all VCA criteria — flood_factor, slope_threshold, cost_threshold have no effect on waterbody pixels. This means:

1. Dilution of LULC change analysis: Lakes are always classified as "Water" in satellite imagery. When stamped as floodplain, they inflate the total floodplain area denominator without contributing to tree loss signal. Tree loss percentages are suppressed.

2. Waterbodies wider than VCA extent: A large lake on a small stream may be 500m wide while VCA at ff=4 only extends 50m from the channel. The entire lake gets stamped as floodplain — adding stable "Water" area that wouldn't be identified as valley floor by the algorithm.

3. Treed wetlands are different: FWA wetland polygons can contain treed land cover that shows real change over time. Excluding wetlands entirely would miss legitimate tree loss in riparian wetlands.

Proposed Solution

Add a waterbodies_report option or output that enables users to understand waterbody contribution:

• Report total floodplain area with and without waterbodies
• Report waterbody area as percentage of total floodplain
• Or: run VCA twice (with/without) and expose both as layers

This doesn't change the algorithm — it gives users the information to interpret their results correctly. When reporting LULC change within the floodplain, users can choose to exclude stable water pixels or adjust denominators.

Context

Discovered during restoration_wedzin_kwa_2024#123 when LULC tree loss percentages appeared diluted. The flooded docs already note: "Pre-filter to valley-bottom features before calling if headwater waterbodies are not wanted." The network filtering (limiting waterbodies to fish-accessible streams) addresses the headwater problem. This issue addresses the dilution problem for correctly-included waterbodies.

Relates to #28
Relates to NewGraphEnvironment/restoration_wedzin_kwa_2024#123
Relates to NewGraphEnvironment/restoration_wedzin_kwa_2024#138

Update 2026-08-27 — the number is cheap to get, and one plausible shortcut is wrong

Measured on the bundled Bulkley tile while building fl_valley_attribute() (#40):

valley without waterbodies : 53,635 cells (536.4 ha)
valley with waterbodies    : 55,345 cells (553.5 ha)
  added by the stamp       :  1,710 cells ( 17.1 ha) -> 3.1% of the total

waterbody cells inside the delineation : 1,842 cells (18.4 ha, 3.3% of delineated area)
  of which the VCA found independently :   132 cells ( 1.3 ha)

So on this tile the stamp adds ~3% of the reported floodplain area, and ~93% of the waterbody
area inside the delineation would not have been found by the VCA at all. Small here; the issue's
"large lake on a small stream" case is where it gets big, and nothing in the current output lets a
reader tell the two situations apart.

Both numbers this issue asks for are a two-line computation — no algorithm change, no second
VCA run needed for the second one:

# with / without
a_with <- fl_valley_confine(dem, streams, ..., waterbodies = wb)
a_without <- fl_valley_confine(dem, streams, ...)          # the only extra run

# waterbody share of the delineation, from a single run
wb_r <- terra::rasterize(terra::vect(wb), dem, field = 1L, background = 0L)
share <- sum(terra::values(wb_r == 1L & a_with == 1L), na.rm = TRUE) /
         sum(terra::values(a_with == 1L), na.rm = TRUE)

The second form is the one worth exposing, because it costs nothing on top of the run the user
already did.

Do not use attr(x, "fl_fallback_cells") from fl_valley_attribute() as this number. It is
close enough on this tile (1,643 cells) to look like the answer and it is not:

• it omits waterbody cells that satisfy the VCA criteria anyway (132 cells here), and
• it includes morphological-cleanup and channel-buffer cells that are not waterbodies at all.

It measures "valley cells no watercourse reaches within thresholds", which merely overlaps with
"waterbody cells". Recording that here because the resemblance is a trap.

What #40 does change: floodplain area is now reportable per watercourse
(fl_valley_attribute()), so the LULC denominator can be per river rather than per watershed group.
If waterbody dilution is being corrected, it should be corrected in the same denominator the change
statistics use — worth settling which granularity the reporting uses before adding a second knob.

Suggested shape, unchanged in spirit from the Proposed Solution above: return the share as an
attribute or an optional attribute column rather than adding a waterbodies_report mode, matching
how fl_valley_attribute() reports fl_fallback_cells.