Add STAC DEM vignette with 1m lidar site-level analysis

.gitignore

@@ -6,3 +6,5 @@ inst/doc
 docs
 /doc/
 /Meta/
+*.html
+Rplots.pdf

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: flooded
 Title: Portable Floodplain Delineation from DEM and Stream Network
-Version: 0.1.0
+Version: 0.1.1
 Authors@R: c(
     person("Allan", "Irvine", , "al@newgraphenvironment.com", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0002-3495-2128")),
@@ -36,6 +36,8 @@ Suggests:
     whitebox,
     testthat (>= 3.0.0),
     knitr,
-    rmarkdown
+    rmarkdown,
+    rstac,
+    gdalcubes
 Config/testthat/edition: 3
 VignetteBuilder: knitr

NEWS.md

@@ -1,3 +1,12 @@
+# flooded 0.1.1
+
+* Add STAC DEM vignette comparing 25 m TRIM (resampled to 10 m) with native
+  1 m lidar — includes site-level zoom and pop-up analysis quantifying
+  anthropogenic barriers to floodplain connectivity.
+* Add `bcdata` reproducibility script (`data-raw/testdata_bcdata.R`).
+* Add resolution and restoration section to README.
+* Pre-build STAC vignette for fast pkgdown rendering.
+
 # flooded 0.1.0
 
 * Initial release. Valley Confinement Algorithm (VCA) pipeline for floodplain

README.md

@@ -44,6 +44,12 @@ flood_depth    = bankfull_depth * flood_factor
 
 Both `upstream_area` (hectares) and `precip` (mean annual precipitation, mm) are important — omitting precipitation underestimates flood depth by ~4x.
 
+## Resolution and restoration
+
+`flooded` is DEM-agnostic — any source works. But resolution changes what you can see. At 25 m (e.g. the provincial TRIM DEM), the floodplain appears as one continuous surface. At 1 m lidar, anthropogenic features emerge: roads, railway grades, dykes, and agricultural fill that sit above the flood surface and block lateral connectivity.
+
+The gap between coarse and fine results is a diagnostic: it shows **what is preventing floodplain from functioning** and **where to act** — removing fill, breaching dykes, or installing crossings to reconnect floodplain. See the [STAC DEM vignette](https://newgraphenvironment.github.io/flooded/articles/stac-dem.html) for a worked example comparing 25 m TRIM with 1 m lidar on the Neexdzii Kwah (Bulkley River).
+
 ## Performance
 
 Set `terra` threads before running the pipeline to enable multi-core processing:

data-raw/testdata_bcdata.R

@@ -0,0 +1,96 @@
+#!/usr/bin/env Rscript
+#
+# testdata_bcdata.R
+#
+# Reproduce the bundled 10m test DEM without the bcfishpass database.
+# Uses bcdata (Python CLI by Simon Norris) to fetch the 25m provincial DEM
+# from BC Data Catalogue WCS, then resamples to 10m with bilinear
+# interpolation — the same approach used by bcfishpass habitat_lateral.
+#
+# This is an internal reproducibility reference, not user-facing.
+# The bundled test data in inst/testdata/ is generated by testdata.R
+# (which reads from the bcfishpass BULK DEM). This script documents
+# an alternative path that doesn't require the database.
+#
+# Prerequisites:
+#   - Python with bcdata installed: pip install bcdata
+#   - R packages: terra, sf
+#
+# Pattern from:
+#   bcfishpass/model/03_habitat_lateral/valley_confinement.py
+#     bcdata.get_dem(bounds, path, as_rasterio=True, align=True)
+#     gdaldem slope -p dem.tif slope.tif
+#
+# Output:
+#   data-raw/dem_bcdata_25m.tif  — raw 25m DEM from BC Data Catalogue
+#   data-raw/dem_bcdata_10m.tif  — resampled to 10m (bilinear)
+#   data-raw/slope_bcdata_10m.tif — percent slope from 10m DEM
+
+library(terra)
+library(sf)
+
+# --- Test area extent (same as inst/testdata/) ---
+# Neexdzii Kwah near Topley, BC Albers EPSG:3005
+xmin <- 975728
+xmax <- 983728
+ymin <- 1054058
+ymax <- 1060538
+
+out_dir <- here::here("data-raw")
+
+# --- Download 25m DEM via bcdata CLI ---
+# bcdata get-dem fetches from BC Data Catalogue WCS
+# bounds format: xmin,ymin,xmax,ymax in EPSG:3005
+dem_25m_path <- file.path(out_dir, "dem_bcdata_25m.tif")
+
+if (!file.exists(dem_25m_path)) {
+  bounds <- paste(xmin, ymin, xmax, ymax, sep = ",")
+  cmd <- paste(
+    "bcdata", "get-dem",
+    "--bounds", bounds,
+    "--src_crs", "EPSG:3005",
+    "--dst_crs", "EPSG:3005",
+    dem_25m_path
+  )
+  message("Downloading 25m DEM via bcdata...")
+  message("  Command: ", cmd)
+  system(cmd)
+} else {
+  message("25m DEM already exists: ", dem_25m_path)
+}
+
+# --- Resample to 10m (bilinear, matching bcfishpass pattern) ---
+dem_10m_path <- file.path(out_dir, "dem_bcdata_10m.tif")
+
+message("Resampling 25m -> 10m (bilinear)...")
+dem_25m <- terra::rast(dem_25m_path)
+
+# Create 10m template matching the test data extent
+template_10m <- terra::rast(
+  xmin = xmin, xmax = xmax,
+  ymin = ymin, ymax = ymax,
+  resolution = 10,
+  crs = "EPSG:3005"
+)
+
+dem_10m <- terra::resample(dem_25m, template_10m, method = "bilinear")
+terra::writeRaster(dem_10m, dem_10m_path, overwrite = TRUE)
+
+# --- Derive slope (percent, matching gdaldem slope -p) ---
+slope_10m_path <- file.path(out_dir, "slope_bcdata_10m.tif")
+
+message("Deriving percent slope...")
+slope_deg <- terra::terrain(dem_10m, "slope", unit = "degrees")
+slope_pct <- tan(slope_deg * pi / 180) * 100
+terra::writeRaster(slope_pct, slope_10m_path, overwrite = TRUE)
+
+# --- Summary ---
+message("\nOutputs:")
+message("  25m DEM: ", dem_25m_path, " (",
+        round(file.size(dem_25m_path) / 1e6, 1), " MB)")
+message("  10m DEM: ", dem_10m_path, " (",
+        round(file.size(dem_10m_path) / 1e6, 1), " MB)")
+message("  Slope:   ", slope_10m_path, " (",
+        round(file.size(slope_10m_path) / 1e6, 1), " MB)")
+message("  Extent:  ", paste(xmin, xmax, ymin, ymax, sep = ", "), " (EPSG:3005)")
+message("  Res:     ", paste(res(dem_10m), collapse = " x "), " m")