diff --git a/.Rbuildignore b/.Rbuildignore index 664c8d5..79b3921 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -11,3 +11,5 @@ ^\.lintr$ ^doc$ ^Meta$ +^planning$ +^\.claude$ diff --git a/.gitignore b/.gitignore index ab70af3..d48b14c 100644 --- a/.gitignore +++ b/.gitignore @@ -8,3 +8,5 @@ docs /Meta/ *.html Rplots.pdf +*.Rcheck/ +*.tar.gz diff --git a/DESCRIPTION b/DESCRIPTION index 64f45d5..fd875ba 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: flooded Title: Portable Floodplain Delineation from DEM and Stream Network -Version: 0.2.1 +Version: 0.3.0 Authors@R: c( person("Allan", "Irvine", , "al@newgraphenvironment.com", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-3495-2128")), @@ -36,11 +36,17 @@ Imports: tibble Suggests: bookdown, - whitebox, - testthat (>= 3.0.0), + DBI, + fresh, + gdalcubes, knitr, rmarkdown, rstac, - gdalcubes + testthat (>= 3.0.0), + whitebox, + xciter +Remotes: + NewGraphEnvironment/fresh, + NewGraphEnvironment/xciter Config/testthat/edition: 3 VignetteBuilder: knitr diff --git a/NAMESPACE b/NAMESPACE index 7cd9c3f..3250178 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,6 +1,7 @@ # Generated by roxygen2: do not edit by hand export(fl_cost_distance) +export(fl_dem_aoi) export(fl_flood_assemble) export(fl_flood_depth) export(fl_flood_model) @@ -16,6 +17,8 @@ export(fl_stream_rasterize) export(fl_valley_confine) export(fl_valley_poly) import(terra) +importFrom(sf,st_buffer) importFrom(sf,st_crs) importFrom(sf,st_read) importFrom(sf,st_transform) +importFrom(sf,st_zm) diff --git a/NEWS.md b/NEWS.md index 0593606..6fc95b0 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,9 @@ +# flooded 0.3.0 + +- New `fl_dem_aoi()` — AOI-driven DEM fetch helper. Defaults to MRDEM-30 via `/vsicurl/` (the rtj-doc-settled choice for watershed-scale BC work) but accepts any local path, `/vsicurl/` URL, or `/vsis3/` S3 URL via `source =`. Buffered crop happens in the source raster's CRS, reprojection after crop. Replaces hand-rolled per-project DEM plumbing (#34). +- New `vignettes/pars-floodplain.Rmd` — watershed-scale showcase running `fl_dem_aoi()` + `fl_valley_confine()` end-to-end on the Parsnip River WSG (5,597 km²). Designed to port to a bookdown report appendix (#34). +- New `data-raw/wsg_vignette_data.R` — generic, parameterised by `wsg <- "PARS"`. Re-runs the full pipeline for any 4-letter BC watershed group, namespaces outputs by WSG code (#34). + # flooded 0.2.1 - Startup quote ritual: `library(flooded)` prints a random fact-checked quote on attach. Italic quote, grey attribution, clickable blue `source` hyperlink (OSC 8). Suppress via `options(flooded.quote_show_source = FALSE)`. diff --git a/R/fl_dem_aoi.R b/R/fl_dem_aoi.R new file mode 100644 index 0000000..462d9c4 --- /dev/null +++ b/R/fl_dem_aoi.R @@ -0,0 +1,106 @@ +#' Fetch a DEM cropped to an AOI +#' +#' Returns a `SpatRaster` of elevation cropped to a buffered AOI. Source is a +#' local file, an HTTP COG (`/vsicurl/`), or an S3 COG (`/vsis3/`). Defaults to +#' MRDEM-30 — NRCan's Medium-Resolution Digital Elevation Model (30 m, all of +#' Canada, public S3, no auth) — the source-of-truth choice for watershed-scale +#' BC floodplain work as decided in `rtj/docs/dem-sources.md`. +#' +#' @param aoi An `sf` or `sfc` object — polygon, lines, or points. Buffered by +#' `buffer` (via [sf::st_buffer()]) before crop. To crop tightly along a +#' stream corridor (memory-efficient when the watershed-scale AOI has sparse +#' stream coverage), pass the streams `sf` here rather than the full WSG +#' polygon. +#' @param source Character. Path or URL for the source raster. `NULL` (default) +#' uses MRDEM-30 DTM via `/vsicurl/`. Local file paths, `/vsicurl/https://...` +#' URLs, and `/vsis3/...` S3 paths all work — `terra::rast()` handles them +#' identically. +#' @param buffer Numeric. Distance in metres to buffer `aoi` before crop. +#' Default `2000`. +#' @param target_crs CRS spec recognised by [sf::st_crs()] (EPSG code, WKT, +#' PROJ string, or `crs` object). `NULL` (default) returns the raster in the +#' input AOI's CRS. **Reprojection happens after crop**, never before — the +#' 84 GB MRDEM COG must not be reprojected as a whole. +#' +#' @return A `SpatRaster` cropped to the buffered AOI extent, projected to +#' `target_crs` if it differs from the source raster's CRS. +#' +#' @details +#' MRDEM-30 (`s3://canelevation-dem/mrdem-30/mrdem-30-dtm.tif`) is a single +#' 84 GB Cloud-Optimized GeoTIFF in EPSG:3979 covering all of Canada. The +#' `/vsicurl/` access pattern range-reads only the bytes intersecting the AOI, +#' so total bandwidth scales with AOI size, not the COG size. +#' +#' For sub-10 m riparian-scale work where lidar coverage exists, query the +#' `stac-dem-bc` STAC catalog and pass an item's COG URL as `source`. See the +#' example below. +#' +#' @seealso [fl_valley_confine()] +#' +#' @examples +#' aoi <- sf::st_read( +#' system.file("testdata/streams.gpkg", package = "flooded"), +#' quiet = TRUE +#' ) +#' +#' # Local file (any path or URL works the same way) +#' dem <- fl_dem_aoi( +#' aoi, +#' source = system.file("testdata/dem.tif", package = "flooded"), +#' buffer = 200 +#' ) +#' terra::plot(dem) +#' +#' \dontrun{ +#' # Default: MRDEM-30 DTM via /vsicurl/ — fetched lazily from NRCan S3 +#' dem <- fl_dem_aoi(aoi, buffer = 2000) +#' terra::plot(dem, main = "MRDEM-30 over AOI") +#' +#' # LidarBC via stac-dem-bc — sub-10 m where lidar coverage exists +#' bbox_4326 <- sf::st_bbox(sf::st_transform(aoi, 4326)) +#' items <- rstac::stac("https://images.a11s.one/") |> +#' rstac::stac_search( +#' collections = "stac-dem-bc", +#' bbox = unname(bbox_4326) +#' ) |> +#' rstac::post_request() |> +#' rstac::items_fetch() +#' if (length(items$features) > 0L) { +#' cog <- paste0("/vsicurl/", items$features[[1]]$assets$image$href) +#' dem_lidar <- fl_dem_aoi(aoi, source = cog, buffer = 100) +#' } +#' } +#' +#' @export +fl_dem_aoi <- function(aoi, source = NULL, buffer = 2000, target_crs = NULL) { + if (!(inherits(aoi, "sf") || inherits(aoi, "sfc"))) { + stop("`aoi` must be an sf or sfc object.", call. = FALSE) + } + stopifnot(is.numeric(buffer), length(buffer) == 1L, buffer >= 0) + + if (is.null(source)) { + source <- paste0( + "/vsicurl/https://canelevation-dem.s3.ca-central-1.amazonaws.com/", + "mrdem-30/mrdem-30-dtm.tif" + ) + } + stopifnot(is.character(source), length(source) == 1L) + + aoi_crs <- sf::st_crs(aoi) + target_crs <- if (is.null(target_crs)) aoi_crs else sf::st_crs(target_crs) + + # GEOS can't buffer XYZM/XYZ — fwapg streams carry route-measure M values + aoi_buf <- sf::st_buffer(sf::st_zm(aoi), dist = buffer) + + r <- terra::rast(source) + r_crs <- sf::st_crs(terra::crs(r)) + + aoi_buf_in_r <- if (r_crs == aoi_crs) aoi_buf else sf::st_transform(aoi_buf, r_crs) + r_clip <- terra::crop(r, terra::vect(aoi_buf_in_r), snap = "out") + + if (r_crs != target_crs) { + r_clip <- terra::project(r_clip, target_crs$wkt) + } + + r_clip +} diff --git a/R/flooded-package.R b/R/flooded-package.R index 5a23719..fb1273f 100644 --- a/R/flooded-package.R +++ b/R/flooded-package.R @@ -1,6 +1,6 @@ #' @keywords internal #' @import terra -#' @importFrom sf st_crs st_transform st_read +#' @importFrom sf st_crs st_transform st_read st_buffer st_zm "_PACKAGE" ## usethis namespace: start diff --git a/data-raw/bib_regenerate.R b/data-raw/bib_regenerate.R new file mode 100644 index 0000000..6752df6 --- /dev/null +++ b/data-raw/bib_regenerate.R @@ -0,0 +1,50 @@ +#!/usr/bin/env Rscript +# +# data-raw/regenerate_bib.R +# +# Regenerate vignettes/references.bib from the union of [@key] markers in +# vignettes/pars-floodplain.Rmd and the citation_keys columns in +# `flooded::fl_params()` and `flooded::fl_scenarios()`. Pulls source +# records from Zotero via Better BibTeX (rbbt -> BBT). +# +# Run after adding/removing [@key] markers in the vignette or after +# changing flood_params.csv / flood_scenarios.csv: +# Rscript data-raw/regenerate_bib.R +# +# Prerequisites: +# - Zotero desktop running with BBT plugin enabled +# - All keys must resolve to items in the Zotero library +# +# CI does not run this script — vignettes/references.bib is committed +# and pkgdown reads the static file. Re-run + commit whenever cites +# change. + +stopifnot(requireNamespace("rbbt", quietly = TRUE)) +stopifnot(requireNamespace("flooded", quietly = TRUE)) + +vignette_path <- here::here("vignettes", "pars-floodplain.Rmd") +out_path <- here::here("vignettes", "references.bib") + +# Keys cited inline in the vignette +keys_vignette <- rbbt::bbt_detect_citations( + paste(readLines(vignette_path), collapse = "\n") +) +message(" vignette: ", length(keys_vignette), " keys") + +# Keys referenced from the package's parameter / scenario tables +keys_tables <- unique(unlist(strsplit( + c(flooded::fl_params()$citation_keys, + flooded::fl_scenarios()$citation_keys), + ";" +))) +keys_tables <- keys_tables[!is.na(keys_tables) & nzchar(keys_tables)] +message(" fl_params + fl_scenarios: ", length(keys_tables), " keys") + +all_keys <- sort(unique(c(keys_vignette, keys_tables))) +message("\nUnion: ", length(all_keys), " unique keys") + +bib_text <- rbbt::bbt_bib(all_keys, .action = rbbt::bbt_return) +writeLines(bib_text, out_path) + +n_entries <- length(grep("^@", readLines(out_path))) +message("Wrote ", n_entries, " entries to ", out_path) diff --git a/data-raw/wsg_vignette_data.R b/data-raw/wsg_vignette_data.R new file mode 100644 index 0000000..1f8932d --- /dev/null +++ b/data-raw/wsg_vignette_data.R @@ -0,0 +1,274 @@ +# Pre-compute heavy data for the WSG floodplain showcase vignette. +# +# Generic — runs for any BC watershed group. Set `wsg` and `species_view` +# below. +# +# Generates (under inst/vignette-data/, namespaced by WSG code): +# .gpkg multi-layer GeoPackage with the sf layers: +# aoi, streams, waterbodies, floodplain, +# railways, roads, reserves, parks, +# named_streams +# _dem.tif MRDEM-30 clip cropped to streams + 2 km buffer +# _valleys.tif fl_valley_confine() binary output (1=valley) +# +# Bundling the sf layers into a single multi-layer gpkg matches how QGIS +# expects a per-WSG project bundle and keeps the cache to three files per +# WSG. Rasters stay separate (terra writes COG/GeoTIFF, not into gpkg). +# +# Streams come from `bcfishpass.streams_bt_vw` (bull trout accessible +# network — every row has `access IN (1, 2)` by construction). Filtering +# to `access IN (1, 2) AND stream_order >= min_order` plus the +# `watershed_group_code` filter gives "best accessible habitat order 3+" +# without needing a working-table classification pipeline. +# +# Waterbodies (lakes + wetlands) are picked up via `waterbody_key` +# linkage from the filtered stream list — only those physically anchored +# to the network are included. They feed `fl_valley_confine(waterbodies =)` +# so the final valley raster fills lake / wetland cells the gradient and +# cost-distance masks would otherwise carve donut holes around. +# +# To swap species: change `species_view` to e.g. `"streams_co_vw"` for +# coho or `"streams_st_vw"` for steelhead — same column shape (`access`, +# `spawning`, `rearing`). +# +# Prerequisites: +# - SSH tunnel to fwapg PostgreSQL up; PG_*_SHARE env vars set +# (defaults of fresh::frs_db_conn()) +# - Outbound HTTPS to s3.ca-central-1.amazonaws.com for MRDEM /vsicurl/ +# - Run from package root so devtools::load_all() finds dev fl_dem_aoi() + +# ---- params ------------------------------------------------------------- + +wsg <- "PARS" # any 4-letter BC watershed group code +species_view <- "streams_bt_vw" # bcfishpass species-accessible view +min_order <- 3 # minimum stream order to keep + +# ---- env ---------------------------------------------------------------- + +Sys.setenv( + GDAL_HTTP_MAX_RETRY = "3", + GDAL_HTTP_RETRY_DELAY = "2", + VSI_CACHE = "TRUE" +) + +devtools::load_all(quiet = TRUE) +library(fresh) +library(sf) +library(terra) + +out_dir <- file.path("inst", "vignette-data") +dir.create(out_dir, recursive = TRUE, showWarnings = FALSE) + +stub <- tolower(wsg) +out <- function(name) file.path(out_dir, paste0(stub, "_", name)) # rasters + +# Single multi-layer GeoPackage for all sf layers. Wipe at the start of +# each run so we never carry stale layers across re-runs. +gpkg <- file.path(out_dir, paste0(stub, ".gpkg")) +if (file.exists(gpkg)) file.remove(gpkg) + +write_layer <- function(x, layer) { + sf::st_write(x, gpkg, layer = layer, delete_layer = TRUE, + append = TRUE, quiet = TRUE) +} + +# ---- 1. WSG boundary ---------------------------------------------------- + +message("Fetching ", wsg, " watershed group boundary ...") +conn <- fresh::frs_db_conn() +on.exit(try(DBI::dbDisconnect(conn), silent = TRUE), add = TRUE) + +aoi <- sf::st_read( + conn, + query = sprintf( + "SELECT watershed_group_code, watershed_group_name, area_ha, geom + FROM whse_basemapping.fwa_watershed_groups_poly + WHERE watershed_group_code = '%s'", + wsg + ), + quiet = TRUE +) +write_layer(aoi, "aoi") +message(sprintf(" area: %.0f km²", aoi$area_ha[1] / 100)) + +# ---- 2. Streams: accessible, order 3+ ----------------------------------- + +message(sprintf("Fetching %s accessible order %d+ streams from bcfishpass.%s ...", + wsg, min_order, species_view)) +streams <- sf::st_read( + conn, + query = sprintf(" + SELECT segmented_stream_id, linear_feature_id, blue_line_key, + waterbody_key, downstream_route_measure, upstream_area_ha, + map_upstream, channel_width, gnis_name, stream_order, + gradient, mapping_code, access, spawning, rearing, + dam_dnstr_ind, watershed_group_code, geom + FROM bcfishpass.%s + WHERE watershed_group_code = '%s' + AND access IN (1, 2) + AND stream_order >= %d", species_view, wsg, min_order), + quiet = TRUE +) +streams <- sf::st_zm(streams) +write_layer(streams, "streams") +message(sprintf(" %d segments (assessed: %d, modelled: %d)", + nrow(streams), + sum(streams$access == 1L, na.rm = TRUE), + sum(streams$access == 2L, na.rm = TRUE))) + +# ---- 3. Waterbodies on the accessible network --------------------------- + +wb_keys <- unique(streams$waterbody_key) +wb_keys <- wb_keys[!is.na(wb_keys) & wb_keys != 0] +message(sprintf("Fetching waterbodies on %d distinct waterbody_keys ...", + length(wb_keys))) + +if (length(wb_keys) > 0L) { + keys_sql <- paste(wb_keys, collapse = ", ") + lakes <- sf::st_read( + conn, + query = sprintf( + "SELECT waterbody_key, geom FROM whse_basemapping.fwa_lakes_poly + WHERE waterbody_key IN (%s)", keys_sql), + quiet = TRUE + ) + wetlands <- sf::st_read( + conn, + query = sprintf( + "SELECT waterbody_key, geom FROM whse_basemapping.fwa_wetlands_poly + WHERE waterbody_key IN (%s)", keys_sql), + quiet = TRUE + ) + lakes_min <- sf::st_sf(waterbody_type = rep("L", nrow(lakes)), + geom = sf::st_geometry(lakes)) + wetlands_min <- sf::st_sf(waterbody_type = rep("W", nrow(wetlands)), + geom = sf::st_geometry(wetlands)) + waterbodies <- rbind(lakes_min, wetlands_min) + # Drop any non-polygon geoms (degenerate features in source data + # occasionally come through as MULTILINESTRING and break terra::rasterize). + keep <- sf::st_geometry_type(waterbodies) %in% c("POLYGON", "MULTIPOLYGON") + if (any(!keep)) { + message(sprintf(" dropping %d non-polygon waterbody geom(s)", sum(!keep))) + waterbodies <- waterbodies[keep, ] + } + message(sprintf(" lakes: %d, wetlands: %d, total: %d", + nrow(lakes), nrow(wetlands), nrow(waterbodies))) +} else { + waterbodies <- sf::st_sf(waterbody_type = character(0), + geom = sf::st_sfc(crs = sf::st_crs(streams))) + message(" no waterbodies on streams") +} +waterbodies <- sf::st_zm(waterbodies) +write_layer(waterbodies, "waterbodies") + +# ---- 4. Context layers (railways, roads, reserves, parks) -------------- + +# Spatial intersect with the WSG boundary. Each query becomes a layer in +# the multi-layer gpkg. Geometry-only attributes plus a name field where +# the source carries one. + +fetch_layer <- function(query_sql, layer_name, label) { + layer <- try(sf::st_read(conn, query = query_sql, quiet = TRUE), + silent = TRUE) + if (inherits(layer, "try-error") || nrow(layer) == 0L) { + message(sprintf(" %s: 0 features (skipping)", label)) + return(invisible(NULL)) + } + layer <- sf::st_zm(layer) + write_layer(layer, layer_name) + message(sprintf(" %s: %d features", label, nrow(layer))) +} + +aoi_wkt <- sf::st_as_text(sf::st_geometry(aoi)) +intersect_clause <- function(geom_col = "geom") { + sprintf("ST_Intersects(%s, ST_GeomFromText('%s', 3005))", geom_col, aoi_wkt) +} + +message("Fetching context layers (railways, roads, reserves, parks) ...") + +fetch_layer( + sprintf("SELECT track_name, geom FROM whse_basemapping.gba_railway_tracks_sp + WHERE %s", intersect_clause()), + "railways", "railways") + +fetch_layer( + sprintf("SELECT transport_line_id, structured_name_1, transport_line_type_code, + highway_route_1, geom + FROM whse_basemapping.transport_line + WHERE %s", intersect_clause()), + "roads", "roads") + +fetch_layer( + sprintf("SELECT english_name, band_name, geom + FROM whse_admin_boundaries.adm_indian_reserves_bands_sp + WHERE %s", intersect_clause()), + "reserves", "First Nations reserves") + +fetch_layer( + sprintf("SELECT protected_lands_name, protected_lands_designation, geom + FROM whse_tantalis.ta_park_ecores_pa_svw + WHERE %s", intersect_clause()), + "parks", "parks / protected") + +# Named streams come pre-filtered by watershed_group_code — no spatial join. +fetch_layer( + sprintf("SELECT gnis_name, blue_line_key, stream_order, geom + FROM whse_basemapping.fwa_named_streams + WHERE watershed_group_code = '%s'", wsg), + "named_streams", "named streams") + +DBI::dbDisconnect(conn) + +# ---- 5. MRDEM-30 clip via fl_dem_aoi() ---------------------------------- + +message("Fetching MRDEM-30 clip (heaviest network step) ...") +dem <- fl_dem_aoi(streams, buffer = 2000) +terra::writeRaster( + terra::as.int(dem), out("dem.tif"), + overwrite = TRUE, + datatype = "INT2S", + gdal = c("COMPRESS=DEFLATE", "PREDICTOR=2", "TILED=YES") +) +message(sprintf(" %d x %d cells (%.1f Mcells)", + ncol(dem), nrow(dem), + ncol(dem) * nrow(dem) / 1e6)) + +# ---- 6. Run flooded VCA pipeline ---------------------------------------- + +message("Running fl_valley_confine() with waterbodies ...") +terra::terraOptions(threads = max(1L, parallel::detectCores() - 2L)) + +valleys <- fl_valley_confine( + dem = dem, + streams = streams, + field = "upstream_area_ha", + precip = fl_stream_rasterize(streams, dem, field = "map_upstream"), + waterbodies = waterbodies, + flood_factor = 4 # ff04 — functional floodplain (recurrent inundation) +) +terra::writeRaster( + valleys, out("valleys.tif"), + overwrite = TRUE, + datatype = "INT1U", + gdal = c("COMPRESS=DEFLATE", "TILED=YES") +) + +message("Polygonizing valleys ...") +floodplain <- fl_valley_poly(valleys) +write_layer(floodplain, "floodplain") + +# ---- 7. Report cache size ----------------------------------------------- + +cache_files <- c(gpkg, list.files(out_dir, + pattern = paste0("^", stub, "_.*\\.tif$"), + full.names = TRUE)) +cache_sizes <- file.info(cache_files)$size +total_mb <- sum(cache_sizes) / 1024^2 +message(sprintf("Cache for %s total: %.1f MB across %d files", + wsg, total_mb, length(cache_files))) +for (i in seq_along(cache_files)) { + message(sprintf(" %-40s %.2f MB", basename(cache_files[i]), + cache_sizes[i] / 1024^2)) +} +message(sprintf(" gpkg layers: %s", + paste(sf::st_layers(gpkg)$name, collapse = ", "))) diff --git a/inst/vignette-data/pars.gpkg b/inst/vignette-data/pars.gpkg new file mode 100644 index 0000000..43458d9 Binary files /dev/null and b/inst/vignette-data/pars.gpkg differ diff --git a/inst/vignette-data/pars_dem.tif b/inst/vignette-data/pars_dem.tif new file mode 100644 index 0000000..fb1bb08 Binary files /dev/null and b/inst/vignette-data/pars_dem.tif differ diff --git a/inst/vignette-data/pars_valleys.tif b/inst/vignette-data/pars_valleys.tif new file mode 100644 index 0000000..ae6e284 Binary files /dev/null and b/inst/vignette-data/pars_valleys.tif differ diff --git a/man/fl_dem_aoi.Rd b/man/fl_dem_aoi.Rd new file mode 100644 index 0000000..55f6908 --- /dev/null +++ b/man/fl_dem_aoi.Rd @@ -0,0 +1,87 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/fl_dem_aoi.R +\name{fl_dem_aoi} +\alias{fl_dem_aoi} +\title{Fetch a DEM cropped to an AOI} +\usage{ +fl_dem_aoi(aoi, source = NULL, buffer = 2000, target_crs = NULL) +} +\arguments{ +\item{aoi}{An \code{sf} or \code{sfc} object — polygon, lines, or points. Buffered by +\code{buffer} (via \code{\link[sf:geos_unary]{sf::st_buffer()}}) before crop. To crop tightly along a +stream corridor (memory-efficient when the watershed-scale AOI has sparse +stream coverage), pass the streams \code{sf} here rather than the full WSG +polygon.} + +\item{source}{Character. Path or URL for the source raster. \code{NULL} (default) +uses MRDEM-30 DTM via \verb{/vsicurl/}. Local file paths, \verb{/vsicurl/https://...} +URLs, and \verb{/vsis3/...} S3 paths all work — \code{terra::rast()} handles them +identically.} + +\item{buffer}{Numeric. Distance in metres to buffer \code{aoi} before crop. +Default \code{2000}.} + +\item{target_crs}{CRS spec recognised by \code{\link[sf:st_crs]{sf::st_crs()}} (EPSG code, WKT, +PROJ string, or \code{crs} object). \code{NULL} (default) returns the raster in the +input AOI's CRS. \strong{Reprojection happens after crop}, never before — the +84 GB MRDEM COG must not be reprojected as a whole.} +} +\value{ +A \code{SpatRaster} cropped to the buffered AOI extent, projected to +\code{target_crs} if it differs from the source raster's CRS. +} +\description{ +Returns a \code{SpatRaster} of elevation cropped to a buffered AOI. Source is a +local file, an HTTP COG (\verb{/vsicurl/}), or an S3 COG (\verb{/vsis3/}). Defaults to +MRDEM-30 — NRCan's Medium-Resolution Digital Elevation Model (30 m, all of +Canada, public S3, no auth) — the source-of-truth choice for watershed-scale +BC floodplain work as decided in \code{rtj/docs/dem-sources.md}. +} +\details{ +MRDEM-30 (\verb{s3://canelevation-dem/mrdem-30/mrdem-30-dtm.tif}) is a single +84 GB Cloud-Optimized GeoTIFF in EPSG:3979 covering all of Canada. The +\verb{/vsicurl/} access pattern range-reads only the bytes intersecting the AOI, +so total bandwidth scales with AOI size, not the COG size. + +For sub-10 m riparian-scale work where lidar coverage exists, query the +\code{stac-dem-bc} STAC catalog and pass an item's COG URL as \code{source}. See the +example below. +} +\examples{ +aoi <- sf::st_read( + system.file("testdata/streams.gpkg", package = "flooded"), + quiet = TRUE +) + +# Local file (any path or URL works the same way) +dem <- fl_dem_aoi( + aoi, + source = system.file("testdata/dem.tif", package = "flooded"), + buffer = 200 +) +terra::plot(dem) + +\dontrun{ +# Default: MRDEM-30 DTM via /vsicurl/ — fetched lazily from NRCan S3 +dem <- fl_dem_aoi(aoi, buffer = 2000) +terra::plot(dem, main = "MRDEM-30 over AOI") + +# LidarBC via stac-dem-bc — sub-10 m where lidar coverage exists +bbox_4326 <- sf::st_bbox(sf::st_transform(aoi, 4326)) +items <- rstac::stac("https://images.a11s.one/") |> + rstac::stac_search( + collections = "stac-dem-bc", + bbox = unname(bbox_4326) + ) |> + rstac::post_request() |> + rstac::items_fetch() +if (length(items$features) > 0L) { + cog <- paste0("/vsicurl/", items$features[[1]]$assets$image$href) + dem_lidar <- fl_dem_aoi(aoi, source = cog, buffer = 100) +} +} + +} +\seealso{ +\code{\link[=fl_valley_confine]{fl_valley_confine()}} +} diff --git a/planning/active/.gitkeep b/planning/active/.gitkeep new file mode 100644 index 0000000..e69de29 diff --git a/planning/archive/2026-05-issue-34-dem-source-helpers/README.md b/planning/archive/2026-05-issue-34-dem-source-helpers/README.md new file mode 100644 index 0000000..0516224 --- /dev/null +++ b/planning/archive/2026-05-issue-34-dem-source-helpers/README.md @@ -0,0 +1,29 @@ +## Outcome + +Shipped `fl_dem_aoi()` — an AOI-driven DEM fetch helper that defaults to MRDEM-30 via `/vsicurl/` but accepts any local path, HTTP COG, or `/vsis3/` URL. Replaces hand-rolled per-project DEM plumbing across the ecosystem. Reproject-after-crop is the load-bearing detail (84 GB MRDEM COG must not be reprojected as a whole). + +Built the **Parsnip River Watershed Group** showcase vignette to exercise the full pipeline end-to-end on a real WSG (5,597 km², 8,436 accessible order-3+ streams from `bcfishpass.streams_bt_vw`, 664 waterbodies, 481.9 km² floodplain at `flood_factor = 4`). The vignette: + +- Loads from a single multi-layer GeoPackage (`pars.gpkg`) plus two GeoTIFFs — generated once locally by the generic `data-raw/wsg_vignette_data.R` (parameterized by `wsg <- "PARS"`) +- Methodology section uses live `fl_params()` / `fl_scenarios()` tables with `knitr::kable` + bookdown captions +- Inline citations rendered via `xciter::xct_keys_to_inline_table_col` against `vignettes/references.bib` (regenerated by `data-raw/bib_regenerate.R` via rbbt → BBT) +- Two maps: full-WSG hero + south-east detail inset, with hydro / floodplain / parks / FN reserves (black diamond + halo'd labels) / roads (FSRs only on hero) / railways / named-stream labels (gq registry colours) +- Direct GitHub raw download links for non-R users + +R CMD check: 0/0/0. 172/172 tests pass. + +**Key learnings worth carrying forward:** + +- Multi-layer GeoPackage is the right shape for a per-WSG bundle (matches QGIS conventions); rasters stay as GeoTIFFs (gpkg raster spec is for tile pyramids, wrong format for analytical DTM) +- `bcfishpass.streams_bt_vw` already classifies access; no need to rebuild the wedzin `frs_break_find` + `frs_classify` pipeline for species-level showcases +- `xciter::xct_keys_to_inline_table_col` + a tiny `keys_to_pandoc()` pre-formatter is the cleanest way to render citation keys from registry CSVs as inline citations within kable tables +- Reproject-after-crop on `/vsicurl/` COGs is the load-bearing pattern for AOI-driven raster fetch — full-COG reprojection is catastrophic at 84 GB + +**Follow-ups filed:** + +- flooded#36 — Harden `fl_valley_confine()` against XYZM streams from fwapg (worked around in data-raw via `sf::st_zm()`) +- soul (TBD) — `/planning-init` slug regex uses GNU `sed`-only `\+`; fix to `sed -E` for macOS BSD +- fresh (TBD) — `frs_break` forwards `points_table`/`points_where`/`aoi` args into `frs_break_find` which doesn't accept them; either re-add them to find or branch the dispatch +- xciter (optional) — promote the `keys_to_pandoc()` helper into the package so `xct_keys_to_inline_table_col` can accept bare semicolon-separated keys without callers pre-formatting + +Closed by: commits c84f301 (PWF baseline), 1304fd4 (Phase 1 — `fl_dem_aoi`), 0b62424 (Phase 2a — data-raw + initial cache), fd8bc98 (Phase 2b — vignette), 4c36261 (Phase 3 — release polish), c20f680 (Phase 4 — vignette refinement). PR pending. diff --git a/planning/archive/2026-05-issue-34-dem-source-helpers/findings.md b/planning/archive/2026-05-issue-34-dem-source-helpers/findings.md new file mode 100644 index 0000000..f439065 --- /dev/null +++ b/planning/archive/2026-05-issue-34-dem-source-helpers/findings.md @@ -0,0 +1,90 @@ +# Findings — Add DEM source helpers (file, STAC) for AOI-driven workflows (#34) + +## Issue context + +### Problem + +`fl_valley_confine()` and related flooded functions accept `SpatRaster` — the user is responsible for sourcing the DEM. Two patterns surface in real use: + +- **Manual file copy** — `restoration_wedzin_kwa_2024/scripts/floodplain_lcc/02_floodplain_model.R:87-98` copies a pre-clipped DEM from a bcfishpass habitat-lateral output (`/Users/airvine/Projects/repo/bcfishpass/model/habitat_lateral/data/temp/BULK/dem.tif`) to the project's GIS folder. One-time, hardcoded to the upper Bulkley project. +- **STAC fetch** — `vignettes/stac-dem.Rmd` shows the rstac + gdalcubes pattern. Illustrative, but production drivers re-roll it. + +Scaling to FWCP peace (Pars WSG to start, then more) means each project shouldn't have to hand-roll AOI-driven DEM fetch. Need helper(s) that take an AOI and return a `SpatRaster`, with the source abstracted. + +### Acceptance (from issue) + +- `fl_dem_aoi()` exported, tested with both file and STAC sources +- Vignette section showing AOI-driven DEM fetch (replaces ad-hoc per-project code) +- `buffer` arg with sensible default + docstring guidance +- Memory note in docstring: very large AOIs may need caller-side chunking (deferred to #B issue) + +### Out of scope (from issue) + +- DEM source registry / config module — defer until N=2+ projects need it +- Slope derivation — wedzin pre-derives it from the same source; separate concern +- Multi-AOI orchestration — covered in companion issue (WSG-scale pipeline) + +## Architecture decisions resolved during planning + +### Source-of-truth: `Projects/repo/rtj/docs/dem-sources.md` + +After exploring the surface (potential new package, fresh-as-DEM-grabber, dedicated dem-bc package, S3 versioning), the rtj doc settled the architecture: + +- **MRDEM-30** is the primary path for watershed-scale flooded work. Single 84 GB COG on `s3://canelevation-dem/mrdem-30/mrdem-30-dtm.tif`, public, no auth, `/vsicurl/` access. `terra::rast(URL)` opens lazily; `terra::crop()` range-reads only the AOI bytes. **No `rstac` needed for the primary path.** +- **LidarBC via `stac-dem-bc`** is the secondary/optional path for sub-10m riparian-scale work. `rstac` already in flooded's Suggests handles this. +- **Don't mirror MRDEM** — NRCan-hosted S3 is more durable than anything we'd build (rtj doc:67). +- **Reproject AFTER crop** (rtj doc:156) — never reproject the full 84 GB COG. Load-bearing gotcha. +- **TRIM archive** (29 WSGs only) is a one-time tarball at `s3://fresh-bc/archive/`, not STAC. Already-extracted local files work via `fl_dem_aoi(source = "/path/to/extracted/dem.tif")`. + +### Buffer semantics resolved by user + +`buffer` applies to whatever sf geometry is passed as `aoi`. Pass a polygon AOI for region-bounded crops; pass streams (sf LINESTRING) for tight, memory-efficient crops along the stream corridor. Wedzin's production optimization (`02_floodplain_model.R:43, 105` — `buf <- 2000; stream_ext <- ext(vect(streams)) + buf`) is exactly the streams-as-aoi case, now first-class via this API. Docstring must call out the memory/time tradeoff. + +### Why no new package, why not fresh + +- fresh is a **network engine** (FWA streams, snapping, traversal) — not a DEM grabber. +- MRDEM is `terra::rast(URL)` with no STAC needed — a single function in flooded is the right home. +- New package (`stacdembc`, `dempot`, etc.) is overhead without payoff when the primary path is one URL constant + 15 LOC. +- LidarBC STAC fetch stays as an example pattern in vignette + roxygen, not a wrapper function. If a second STAC-driven workflow appears later, promote to function. + +## flooded codebase patterns confirmed + +From the Explore agent survey: + +- `fl_*` exports, one function per file, mirrored test in `tests/testthat/` +- SpatRaster-first arg convention; `stopifnot()` validation; return SpatRaster +- `@import terra` and `@importFrom sf ...` centralized in `R/flooded-package.R:2-3` +- Examples use `system.file()` via `testdata_path()` helper in `tests/testthat/setup.R:2-4` +- `inst/testdata/` has `dem.tif` (10m BC Albers), `slope.tif`, `streams.gpkg`, `waterbodies.gpkg` +- `vignettes/stac-dem.Rmd` already exists with the LidarBC STAC pattern at lines 95-184 — leave untouched + +## Vignette caching pattern from `cd` + +`~/Projects/repo/cd/data-raw/peace_fwcp_vignette_data.R` and `~/Projects/repo/cd/vignettes/peace-fwcp.Rmd` are the reference: + +- data-raw script runs heavy network/computation work locally (regional climate analysis with ~144 COG range requests) +- Caches outputs as small .rds (271 KB for Peace) + small .tif (6 KB for Peace) to `inst/vignette-data/` +- **Total bundled output target <1 MB** — cd's full Peace + Kootenay analysis fits in 692 KB +- Vignette frontmatter: `bookdown::html_document2`, `bibliography: references.bib`, `link-citations: true` — directly portable to a Peace report appendix + +For flooded's PARS vignette, the cached outputs are mostly sf objects (highly compressible) plus a small DEM tile for the hero map. Realistic cache budget 5-15 MB. + +## Critical reference paths + +- `Projects/repo/rtj/docs/dem-sources.md:50-65` — exact MRDEM access pattern to lift verbatim +- `Projects/repo/rtj/docs/dem-sources.md:67` — "do not mirror MRDEM" rationale +- `Projects/repo/rtj/docs/dem-sources.md:156` — reproject-after-crop gotcha +- `Projects/repo/cd/data-raw/peace_fwcp_vignette_data.R` — data-raw cache pattern +- `Projects/repo/cd/vignettes/peace-fwcp.Rmd` — vignette structure for report-appendix portability +- `Projects/repo/restoration_wedzin_kwa_2024/scripts/floodplain_lcc/02_floodplain_model.R:43, 87-107` — wedzin's manual pattern this replaces +- `Projects/repo/flooded/vignettes/stac-dem.Rmd:95-184` — existing LidarBC STAC pattern (left untouched, referenced in roxygen example) + +## Bugs to file separately + +### Planning-init slug regex (soul) + +Planning-init skill's slug regex uses GNU `sed`-only `\+` syntax — fails silently on macOS BSD sed. Produced `34-add dem source helpers (file, stac) for` with literal spaces and parens until I switched to `sed -E`. Out of scope for #34; will file against soul. + +### `fl_valley_confine()` XYZM robustness (flooded follow-up) + +Streams fetched from fwapg (e.g., via `fresh::frs_stream_fetch()`) carry XYZM geometries — route measures live in M. `fl_valley_confine()` calls `sf::st_buffer()` internally (channel_buffer feature, R/fl_valley_confine.R:202) which fails on XYZM with `GEOS does not support XYM or XYZM geometries`. Every caller pulling from fwapg has to remember to `sf::st_zm()` first. One-line fix: drop Z/M at the top of `fl_valley_confine()` once streams type is sf. Worth filing as its own issue; out of scope for #34, fixed by `st_zm()` call in the data-raw script for now. (`fl_dem_aoi()` was already hardened during this work.) diff --git a/planning/archive/2026-05-issue-34-dem-source-helpers/progress.md b/planning/archive/2026-05-issue-34-dem-source-helpers/progress.md new file mode 100644 index 0000000..84e9f0b --- /dev/null +++ b/planning/archive/2026-05-issue-34-dem-source-helpers/progress.md @@ -0,0 +1,36 @@ +# Progress — Add DEM source helpers (file, STAC) for AOI-driven workflows (#34) + +## Session 2026-05-05 + +- Plan-mode exploration: + - Surveyed flooded `fl_*` conventions, STAC vignette pattern, test scaffolding via Explore agent + - Surveyed wedzin DEM handling pattern (`02_floodplain_model.R`) via Explore agent + - Read `Projects/repo/rtj/docs/dem-sources.md` to ground architecture in settled rtj decisions + - Discussed scope/architecture with user: rejected fresh-as-DEM-grabber, rejected new package, settled on single `fl_dem_aoi()` in flooded with MRDEM-30 inline default + - Confirmed PARS vignette deliverable feasible (cd cache pattern proves region-scale at <1 MB) +- Phases approved by user +- Created branch `34-add-dem-source-helpers-file-stac-for-aoi` off main +- Scaffolded PWF baseline (`task_plan.md`, `findings.md`, `progress.md`) with approved phases +- Phase 1 complete: `fl_dem_aoi()` implemented with MRDEM-30 default, file/`/vsicurl/`/`/vsis3/` source dispatch, CRS-safe crop-before-reproject. 8 tests pass including live MRDEM range read; lintr clean; full suite 172 PASS. +- Phase 2a complete: `data-raw/wsg_vignette_data.R` is generic (any WSG via `wsg <- "PARS"`), pulls boundary from `whse_basemapping.fwa_watershed_groups_poly` and streams from `bcfishpass.streams_vw` (carries channel_width / upstream_area_ha / map_upstream), runs `fl_dem_aoi()` against MRDEM-30 via /vsicurl/, runs `fl_valley_confine()` end-to-end, polygonizes. Two XYZM bugs surfaced: fixed in `fl_dem_aoi()` itself, worked around for `fl_valley_confine()` via `sf::st_zm()` in the script (logged as follow-up). +- Cache for PARS: 11 MB total — `pars_aoi.gpkg` 0.31 MB, `pars_streams.gpkg` 2.7 MB, `pars_dem.tif` 6.2 MB (INT2S, was 29 MB FLT4S), `pars_valleys.tif` 88 KB, `pars_floodplain.gpkg` 1.2 MB. +- Phase 2b complete: `vignettes/pars-floodplain.Rmd` written using `bookdown::html_vignette2` (matches existing flooded vignettes; YAML swappable for Peace report appendix). Loads cached PARS data, walks through `fl_dem_aoi()` + `fl_valley_confine()`, renders a hillshade + floodplain hero map, summary stats (449 km² floodplain on order 4+ streams, 8% of PARS). Renders cleanly to 604 KB HTML; full test suite still 172 PASS. +- Phase 3 in progress: document + lintr + tests clean. devtools::check originally reported 2 NOTEs (`top-level files`, `vignettes/figure leftover`). Fixed: added `^planning$` + `^\.claude$` to `.Rbuildignore` (cleared top-level NOTE), deleted `vignettes/figure/` (knitr leftover from earlier vignette renders, cleared vignettes NOTE). Filed `fl_valley_confine` XYZM hardening as flooded#36 — out of scope for #34, worked around in data-raw script. NEWS entry under 0.3.0; DESCRIPTION bumped 0.2.1 → 0.3.0. + +## Session 2026-05-07 (Phase 4 — vignette refinement) + +Iterated heavily on the PARS vignette and the data-raw script in response to user feedback: + +- **Stream source switched** from generic `bcfishpass.streams_vw` to **`bcfishpass.streams_bt_vw`** (bull trout accessible network — every row already has `access IN (1, 2)` by construction). Replaces an attempted wedzin-style `frs_network` + `frs_break_find` + `frs_classify` pipeline that hit a bug in `fresh::frs_break` (forwards `points_table`/`points_where` args into `frs_break_find` which doesn't accept them — filed for upstream fresh fix). +- **min_order lowered** 4 → 3 to capture more habitat-relevant tributaries. +- **Waterbodies added** (lakes + wetlands via `waterbody_key` linkage) and passed to `fl_valley_confine(waterbodies = ...)` so the binary valley raster fills lake / wetland cells. +- **Flood factor set to 4** (ff04 — functional floodplain). Floodplain area: 481.9 km² (8.6% of PARS). +- **Context layers added**: railways (`gba_railway_tracks_sp`), roads (`transport_line` — filtered to `RRS`/`RRD`/`RRN` resource roads on hero map, all on inset), First Nations reserves (`adm_indian_reserves_bands_sp`, with `english_name` for labels and `band_name` retained), parks / protected (`ta_park_ecores_pa_svw`), named streams (`fwa_named_streams`, labels on inset only). +- **Multi-layer GeoPackage consolidation**: 9 sf layers now live in a single `pars.gpkg` instead of nine separate files. Rasters stay as standalone GeoTIFFs (the gpkg raster spec uses tile pyramids — wrong format for analytical DEM / binary valleys). Cache: 14 MB across 3 files. +- **Methodology section bumped to top** with ecological context + algorithm provenance; `fl_params()` and `fl_scenarios()` rendered as `knitr::kable` tables with `bookdown::html_vignette2` captions. Citation keys converted to inline citations via `xciter::xct_keys_to_inline_table_col()` with a `keys_to_pandoc()` helper that pre-formats bare semicolon-separated keys to pandoc syntax. +- **References** wired via `bibliography: references.bib` + `data-raw/bib_regenerate.R` (rbbt → BBT, mirrors the `cd` repo pattern). Generates `vignettes/references.bib` from the union of vignette `[@key]` markers and the `citation_keys` columns in `fl_params()` / `fl_scenarios()`. +- **Cartography polish**: gq registry colours for parks (`#639b5f`) and reserves (`#b2b2b2`); reserves drawn with black diamond markers + thin-halo'd `english_name` labels; layers drawn ABOVE hydro/transport so they read on overlap; all layers `fresh::frs_clip()`'d to the WSG and DEM `terra::mask()`'d to the WSG so nothing renders outside the basin; `xpd = TRUE` on labels so edge-of-AOI text isn't cropped. +- **Direct download links** for `pars.gpkg`, `pars_dem.tif`, `pars_valleys.tif` from the GitHub raw URL — readers can grab the bundle without installing R. +- **DESCRIPTION** Suggests now includes `DBI`, `fresh`, `xciter` (vignette uses) plus `Remotes:` for the two NewGraphEnvironment-only deps. Cleared the `unstated dependencies in vignettes` R CMD check NOTE. + +Final R CMD check: 0 errors, 0 warnings, 0 notes. 172/172 tests pass. Vignette renders to ~1.7 MB HTML. diff --git a/planning/archive/2026-05-issue-34-dem-source-helpers/task_plan.md b/planning/archive/2026-05-issue-34-dem-source-helpers/task_plan.md new file mode 100644 index 0000000..d992c8e --- /dev/null +++ b/planning/archive/2026-05-issue-34-dem-source-helpers/task_plan.md @@ -0,0 +1,94 @@ +# Task: Add DEM source helpers (file, STAC) for AOI-driven workflows (#34) + +## Problem + +`fl_valley_confine()` and related flooded functions accept `SpatRaster` — the +user is responsible for sourcing the DEM. Two patterns surface in real use +(manual file copy hardcoded to project paths, and an illustrative STAC vignette +that production drivers re-roll). Scaling to FWCP Peace (Pars WSG to start, +then more) means each project shouldn't have to hand-roll AOI-driven DEM fetch. + +## Resolved approach + +Single function — `fl_dem_aoi(aoi, source = NULL, buffer = 2000, target_crs = NULL)` +— that defaults to MRDEM-30 via `/vsicurl/` (the architecture settled in +`Projects/repo/rtj/docs/dem-sources.md`). LidarBC via STAC stays as an example +pattern in roxygen/vignettes, not a wrapper function. PARS WSG floodplain +showcase vignette is the big deliverable, built so it can port to the Peace +report appendix. + +## Phase 1 — Implement `fl_dem_aoi()` + +- [x] Create `R/fl_dem_aoi.R` with function body lifting `rtj/docs/dem-sources.md:50-65` +- [x] Roxygen with two examples: default MRDEM-30 fetch, and `\dontrun{}` LidarBC-via-STAC override +- [x] Update `R/flooded-package.R` if any new `@importFrom` needed +- [x] Create `tests/testthat/test-fl_dem_aoi.R`: + - [x] File mode happy path with bundled `dem.tif` + - [x] CRS mismatch (AOI in 4326, raster in 3005) — function transforms correctly + - [x] `target_crs` override reprojects after crop + - [x] MRDEM `/vsicurl/` mode gated by `skip_on_cran()` + `skip_if_offline()` + - [x] Buffer semantics — larger buffer yields larger extent +- [x] `lintr::lint_package()` clean +- [x] `devtools::document()` regenerates NAMESPACE / man cleanly + +## Phase 2 — Floodplain showcase vignette + data-raw cache + +Generic data-raw script that runs for any WSG by changing one variable +(`wsg <- "PARS"`); PARS is the worked example for the vignette and the +Peace report appendix. + +- [x] Create `data-raw/wsg_vignette_data.R` (generic — `wsg` parameter at top, output filenames namespaced): + - [x] Fetch WSG boundary + streams from fwapg via `fresh::frs_*` + - [x] Run `fl_dem_aoi()` against MRDEM-30 (streams-as-aoi pattern for tight crop) + - [x] Drop XYZM in `fl_dem_aoi()` so fwapg streams flow through GEOS — bug fix + - [x] Drop XYZM on streams in script before `fl_valley_confine()` (workaround pending follow-up issue) + - [x] Run `fl_valley_confine()` end-to-end + - [x] Cache outputs to `inst/vignette-data/` as INT2S DEM (target <15 MB total) +- [x] Run data-raw locally for PARS; cache 14 MB across 3 files (multi-layer gpkg + 2 GeoTIFFs) +- [x] Refactor data-raw to bcfishpass-direct accessible-habitat extraction: + - [x] Use `bcfishpass.streams_bt_vw` directly — every row already has `access IN (1, 2)` so no working-table classification needed + - [x] Lower `min_order` from 4 to 3 + - [x] Pull waterbodies via `waterbody_key` linkage on the accessible network (lakes + wetlands) + - [x] Pass `waterbodies` to `fl_valley_confine()` (fills lake/wetland donut holes) + - [x] Set `flood_factor = 4` (ff04 — functional floodplain scenario) +- [x] Update vignette to: + - [x] Add ecological context paragraph + methodology section at top with `fl_params()` + `fl_scenarios()` tables (kable + bookdown captions) + - [x] Inline citations + `vignettes/references.bib` regenerated by `data-raw/bib_regenerate.R` (rbbt → BBT) + - [x] Show waterbodies, parks, FN reserves, roads, railways, named-stream labels on hero + zoom maps + - [x] FN reserves with black diamond markers + halo'd `english_name` labels (gq registry colours) + - [x] Add zoomed inset map of the south-east corner with full detail + - [x] Direct GitHub raw download links for `pars.gpkg` + 2 GeoTIFFs +- [x] Consolidate sf cache into single multi-layer gpkg (`pars.gpkg` with 9 layers); rasters stay as GeoTIFF +- [x] Create `vignettes/pars-floodplain.Rmd` (using `bookdown::html_vignette2` to match existing flooded style; YAML easily swappable for Peace report appendix bookdown config) + - [x] Why MRDEM-30 (link to rtj doc rationale) + - [x] AOI definition (PARS WSG, framed as "any WSG via the generic data-raw script") + - [x] `fl_dem_aoi(aoi = streams, buffer = 2000)` — streams-as-aoi pattern + - [x] `fl_valley_confine()` walkthrough + - [x] Result map: hillshade + valleys overlay + streams + WSG boundary via terra::plot + - [x] Summary stats: floodplain area (449 km², 8% of PARS), stream + polygon counts +- [x] Confirm vignette renders cleanly (`rmarkdown::render` produces 604 KB HTML) +- [x] Existing `vignettes/stac-dem.Rmd` left untouched (LidarBC reference) + +## Phase 3 — Polish + release prep + +- [x] `devtools::document()` — clean +- [x] `lintr::lint_package()` — 0 lints in new files (3 pre-existing in fl_flood_surface, fl_valley_poly, zzz; not touching adjacent code) +- [x] `devtools::test()` — 172/172 PASS including live MRDEM `/vsicurl/` fetch +- [x] `devtools::check()` — Status: OK (0 errors, 0 warnings, 0 notes) after `.Rbuildignore` updates and `vignettes/figure/` cleanup +- [x] `du -sh inst/vignette-data/` — 11 MB, under 15 MB budget +- [x] `NEWS.md` entry: `fl_dem_aoi()`, PARS vignette, generic data-raw script — all linked to #34 +- [x] `DESCRIPTION` version bump `0.2.1` → `0.3.0` +- [x] `.Rbuildignore` — added `^planning$` and `^\.claude$` + +## Validation + +- [x] All tests pass — 172/172 PASS +- [x] `/code-check` clean on each commit (c84f301, 1304fd4, 0b62424, fd8bc98) +- [x] PWF checkboxes match landed work +- [ ] `/planning-archive` on completion +- [ ] `/gh-pr-push` opens PR with `Fixes #34` + +## Follow-ups filed + +- flooded#36 — Harden `fl_valley_confine()` against XYZM streams from fwapg (worked around for #34 via `sf::st_zm()` in data-raw script) +- soul (TBD) — `/planning-init` slug regex uses GNU `sed`-only `\+`; fix to `sed -E` for macOS BSD compatibility diff --git a/tests/testthat/test-fl_dem_aoi.R b/tests/testthat/test-fl_dem_aoi.R new file mode 100644 index 0000000..8f047c6 --- /dev/null +++ b/tests/testthat/test-fl_dem_aoi.R @@ -0,0 +1,86 @@ +test_that("fl_dem_aoi crops local file to buffered AOI", { + aoi <- sf::st_read(testdata_path("streams.gpkg"), quiet = TRUE) + + result <- fl_dem_aoi(aoi, source = testdata_path("dem.tif"), buffer = 100) + + expect_s4_class(result, "SpatRaster") + expect_true(terra::ncell(result) > 0) + + # Result extent must cover the buffered AOI + aoi_buf_ext <- terra::ext(terra::vect(sf::st_buffer(aoi, 100))) + res_ext <- terra::ext(result) + expect_true(res_ext$xmin <= aoi_buf_ext$xmin) + expect_true(res_ext$xmax >= aoi_buf_ext$xmax) + expect_true(res_ext$ymin <= aoi_buf_ext$ymin) + expect_true(res_ext$ymax >= aoi_buf_ext$ymax) +}) + +test_that("fl_dem_aoi handles AOI/source CRS mismatch", { + aoi_3005 <- sf::st_read(testdata_path("streams.gpkg"), quiet = TRUE) + aoi_4326 <- sf::st_transform(aoi_3005, 4326) + + result <- fl_dem_aoi(aoi_4326, source = testdata_path("dem.tif"), buffer = 100) + + expect_s4_class(result, "SpatRaster") + expect_true(terra::ncell(result) > 0) + # target_crs defaulted to aoi_4326's CRS (EPSG:4326) + expect_equal(sf::st_crs(terra::crs(result))$epsg, 4326L) +}) + +test_that("fl_dem_aoi target_crs override reprojects after crop", { + aoi <- sf::st_read(testdata_path("streams.gpkg"), quiet = TRUE) + + result <- fl_dem_aoi( + aoi, + source = testdata_path("dem.tif"), + buffer = 100, + target_crs = 4326 + ) + + expect_s4_class(result, "SpatRaster") + expect_equal(sf::st_crs(terra::crs(result))$epsg, 4326L) +}) + +test_that("fl_dem_aoi larger buffer yields equal-or-larger extent", { + aoi <- sf::st_read(testdata_path("streams.gpkg"), quiet = TRUE) + + small <- fl_dem_aoi(aoi, source = testdata_path("dem.tif"), buffer = 100) + big <- fl_dem_aoi(aoi, source = testdata_path("dem.tif"), buffer = 500) + + expect_true(terra::ncell(big) >= terra::ncell(small)) +}) + +test_that("fl_dem_aoi accepts polygon AOI", { + waterbodies <- sf::st_read(testdata_path("waterbodies.gpkg"), quiet = TRUE) + + result <- fl_dem_aoi( + waterbodies, + source = testdata_path("dem.tif"), + buffer = 100 + ) + + expect_s4_class(result, "SpatRaster") + expect_true(terra::ncell(result) > 0) +}) + +test_that("fl_dem_aoi rejects non-sf aoi", { + expect_error( + fl_dem_aoi("not_sf", source = testdata_path("dem.tif")), + "sf" + ) +}) + +test_that("fl_dem_aoi fetches MRDEM-30 via /vsicurl/", { + testthat::skip_on_cran() + testthat::skip_if_offline() + + # Single-feature AOI keeps the range read tiny (~100s of KB) + aoi <- sf::st_read(testdata_path("streams.gpkg"), quiet = TRUE)[1, ] + + result <- fl_dem_aoi(aoi, buffer = 100) + + expect_s4_class(result, "SpatRaster") + expect_true(terra::ncell(result) > 0) + # Default target_crs is the input AOI's CRS + expect_equal(sf::st_crs(terra::crs(result)), sf::st_crs(aoi)) +}) diff --git a/vignettes/figure/plot-compare-1.png b/vignettes/figure/plot-compare-1.png deleted file mode 100644 index ccbf43d..0000000 Binary files a/vignettes/figure/plot-compare-1.png and /dev/null differ diff --git a/vignettes/figure/site-compare-1.png b/vignettes/figure/site-compare-1.png deleted file mode 100644 index 95c1fce..0000000 Binary files a/vignettes/figure/site-compare-1.png and /dev/null differ diff --git a/vignettes/pars-floodplain.Rmd b/vignettes/pars-floodplain.Rmd new file mode 100644 index 0000000..a9a0163 --- /dev/null +++ b/vignettes/pars-floodplain.Rmd @@ -0,0 +1,463 @@ +--- +title: "Watershed-scale floodplain delineation — Parsnip River Watershed Group" +output: + bookdown::html_vignette2 +bibliography: references.bib +link-citations: false +vignette: > + %\VignetteIndexEntry{Watershed-scale floodplain delineation — Parsnip River Watershed Group} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r setup, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + fig.width = 7, + fig.height = 6, + dpi = 150, + message = FALSE, + warning = FALSE +) +library(flooded) +``` + +This vignette runs the `flooded` pipeline on the Parsnip River Watershed +Group (`PARS`, 5,597 km², north-eastern BC) and demonstrates the +AOI-driven helper [`flooded::fl_dem_aoi()`][fl_dem_aoi]. + +PARS sits between Prince George and Mackenzie, BC. The Parsnip River +flows north and enters the southern arm of Williston Reservoir at +Mackenzie, joining the Peace River system. From there the drainage runs +Peace → Slave → Mackenzie River, ultimately discharging to the Arctic +Ocean via the Mackenzie Delta. + +[fl_dem_aoi]: ../reference/fl_dem_aoi.html + +## Why floodplains + +Floodplains are the low-lying, periodically inundated lands flanking a +stream — the parts of the valley that water reclaims when discharge +exceeds bankfull. They store flood water and dissipate flood energy, +recharge shallow groundwater, sort sediment, and host the off-channel +ponds, side channels, sloughs, lakes, and wetlands that drive +disproportionate ecological productivity per unit area. For salmonids +in particular, floodplain habitat — slow-water margins, beaver +complexes, off-channel rearing — is often the limiting factor for +freshwater survival, especially through winter low-flow and high-flow +refugia. Mapping the *extent* of the floodplain at watershed scale is +therefore the first step in scoping where restoration, protection, or +flood-risk planning has the most leverage. + +## Modelling parameters + +`flooded::fl_valley_confine()` combines four masks — slope, distance +from stream, cost-distance through the terrain, and a bankfull flood +model — and applies morphological cleanup (closing, hole-fill, patch +removal) before returning a binary valley raster. The algorithm is +adapted from the +[USDA Valley Confinement Algorithm](https://github.com/Ecotrust/BlueGeo) +[@nagel_etal2014LandscapeScale] (BlueGeo R port by Devin Cairns, MIT) +and the [bcfishpass](https://github.com/smnorris/bcfishpass) lateral +habitat assembly (Simon Norris, Apache 2.0). The bankfull flood model +follows the regional regression of @hall_etal2007Predictingriver. The +parameter legend (units, defaults, citations, ecological effect) lives +in `flooded::fl_params()`. + +This run uses the **`ff04`** scenario from `flooded::fl_scenarios()` — +*functional floodplain*, the recurrent-inundation footprint. The active +parameter values are shown in Table \@ref(tab:params-stamp). + +```{r citation-helper, include = FALSE} +# Convert "key1;key2" bare-key strings (the format used in +# `fl_params()$citation_keys` and `fl_scenarios()$citation_keys`) into +# pandoc citation markers `[@key1; @key2]`. The output is then handed to +# `xciter::xct_keys_to_inline_table_col()` which calls pandoc to render +# them as proper inline citations, against the same references.bib that +# pandoc uses for the rest of the vignette. +keys_to_pandoc <- function(s) { + vapply(s, function(x) { + if (is.na(x) || !nzchar(x)) return("") + paste0("[", paste0("@", trimws(strsplit(x, ";")[[1]]), collapse = "; "), "]") + }, character(1), USE.NAMES = FALSE) +} +``` + +```{r params-stamp, echo = FALSE} +ff_used <- 4 +pars_used <- c("flood_factor", "slope_threshold", "max_width", "cost_threshold") +p_tab <- flooded::fl_params()[match(pars_used, flooded::fl_params()$parameter), + c("parameter", "default", "unit", "effect", + "citation_keys")] +p_tab$value <- p_tab$default +p_tab$value[p_tab$parameter == "flood_factor"] <- ff_used +p_tab$citation_keys <- keys_to_pandoc(p_tab$citation_keys) +p_tab <- xciter::xct_keys_to_inline_table_col( + p_tab, col_format = "citation_keys", + path_bib = "references.bib" +) +p_tab <- p_tab[, c("parameter", "value", "unit", "effect", "citation_keys")] +names(p_tab)[names(p_tab) == "citation_keys"] <- "source" +knitr::kable( + p_tab, row.names = FALSE, + caption = "Active parameter values for the PARS run. `flood_factor = 4` (`ff04`); other values are package defaults from `flooded::fl_params()`." +) +``` + +The package ships three pre-baked scenarios (Table +\@ref(tab:scenarios-stamp)). They differ only in `flood_factor`; all +other parameters are held constant so output differences isolate the +ecological signal. + +```{r scenarios-stamp, echo = FALSE} +s_tab <- flooded::fl_scenarios()[, c("scenario_id", "flood_factor", + "description", "ecological_process", + "citation_keys")] +s_tab$citation_keys <- keys_to_pandoc(s_tab$citation_keys) +s_tab <- xciter::xct_keys_to_inline_table_col( + s_tab, col_format = "citation_keys", + path_bib = "references.bib" +) +names(s_tab)[names(s_tab) == "citation_keys"] <- "source" +knitr::kable( + s_tab, row.names = FALSE, + caption = "Pre-baked flood-factor scenarios shipped with the package. Switch by passing `flood_factor =` to `flooded::fl_valley_confine()`." +) +``` + +For a side-by-side comparison of all three scenarios on a smaller reach, +see the +[valley-confinement vignette](https://newgraphenvironment.github.io/flooded/articles/valley-confinement.html). + +## Cached inputs + +The vignette renders against pre-computed outputs from +`data-raw/wsg_vignette_data.R` so it builds fast and doesn't touch the +network or database. The hydro / context layers ship as a single +multi-layer GeoPackage and the rasters as separate GeoTIFFs (raster +tiles in GPKG would lose the continuous DTM precision and the binary +valleys semantics — wrong format for analytical layers). + +Direct downloads of the cached PARS bundle from the repo (open in QGIS +or any GDAL-aware tool): + +- [`pars.gpkg`](https://github.com/NewGraphEnvironment/flooded/raw/main/inst/vignette-data/pars.gpkg) + — vectors: `aoi`, `streams`, `waterbodies`, `floodplain`, `railways`, + `roads`, `reserves`, `parks`, `named_streams` +- [`pars_dem.tif`](https://github.com/NewGraphEnvironment/flooded/raw/main/inst/vignette-data/pars_dem.tif) + — MRDEM-30 clip (Int16) +- [`pars_valleys.tif`](https://github.com/NewGraphEnvironment/flooded/raw/main/inst/vignette-data/pars_valleys.tif) + — `fl_valley_confine()` binary output + +```{r load, include = FALSE} +library(terra) +library(sf) + +gpkg <- system.file("vignette-data/pars.gpkg", package = "flooded", mustWork = TRUE) +dem_tif <- system.file("vignette-data/pars_dem.tif", package = "flooded", mustWork = TRUE) +val_tif <- system.file("vignette-data/pars_valleys.tif", package = "flooded", mustWork = TRUE) + +present <- sf::st_layers(gpkg)$name +load_layer <- function(name) { + if (name %in% present) sf::st_read(gpkg, layer = name, quiet = TRUE) else NULL +} + +aoi <- load_layer("aoi") +streams <- load_layer("streams") +waterbodies <- load_layer("waterbodies") +floodplain <- load_layer("floodplain") +railways <- load_layer("railways") +roads <- load_layer("roads") +reserves <- load_layer("reserves") +parks <- load_layer("parks") +named_streams <- load_layer("named_streams") +dem <- terra::rast(dem_tif) +valleys <- terra::rast(val_tif) + +# Trim the DEM to the WSG and clip every layer to the same boundary so +# nothing renders outside the basin. +dem <- terra::mask(terra::crop(dem, terra::vect(aoi)), terra::vect(aoi)) +valleys <- terra::mask(terra::crop(valleys, terra::vect(aoi)), terra::vect(aoi)) +streams <- fresh::frs_clip(streams, aoi) +waterbodies <- fresh::frs_clip(waterbodies, aoi) +floodplain <- fresh::frs_clip(floodplain, aoi) +if (!is.null(roads)) roads <- fresh::frs_clip(roads, aoi) +if (!is.null(railways)) railways <- fresh::frs_clip(railways, aoi) +if (!is.null(reserves)) reserves <- fresh::frs_clip(reserves, aoi) +if (!is.null(parks)) parks <- fresh::frs_clip(parks, aoi) +if (!is.null(named_streams)) named_streams <- fresh::frs_clip(named_streams, aoi) +``` + +## DEM + +The underlying elevation data is +[MRDEM-30](https://open.canada.ca/data/en/dataset/18752265-bda3-498c-a4ba-9dfe68cb98da), +NRCan's 30 m Medium-Resolution Digital Elevation Model — a single +Cloud-Optimized GeoTIFF covering all of Canada, hosted on public S3 +with no authentication required, integrated from LidarBC where lidar +coverage exists with Copernicus TanDEM-X / CDEM-derived fallback +elsewhere. `flooded::fl_dem_aoi()` reads it via `/vsicurl/`, so the only +bytes transferred are those intersecting the AOI — bandwidth scales +with the AOI, not with the COG size. + +The DEM was fetched in `data-raw/wsg_vignette_data.R` with a single +call, passing the streams as the AOI so the crop hugs the stream +corridor (a memory-efficient pattern for large watersheds with sparse +stream networks): + +```{r dem-fetch, eval = FALSE} +dem <- flooded::fl_dem_aoi(streams, buffer = 2000) +``` + +The default `source = NULL` resolves to the canonical MRDEM-30 DTM +`/vsicurl/` URL inside `flooded::fl_dem_aoi()`. `buffer = 2000` extends +the AOI by 2 km in metres before crop. To override the source — e.g., +to fetch a LidarBC COG — pass `source = "/vsicurl/https://.../tile.tif"`. + +## Streams and waterbodies + +Streams are habitat segments modelled as accessible to bull trout from +[bcfishpass](https://github.com/smnorris/bcfishpass) outputs. bcfishpass +runs weekly on a hosted virtual machine and republishes the +`streams_bt_vw` view, so the data is current to within a week of any +upstream FWA, observation, or barrier change. + +```{r bcfishpass-version, echo = FALSE} +conn <- fresh::frs_db_conn() +ver <- fresh::frs_db_query(conn, + "SELECT model_version, date_completed + FROM bcfishpass.log + WHERE model_type = 'LINEAR' + ORDER BY date_completed DESC LIMIT 1") +DBI::dbDisconnect(conn) +``` + +This run is built against bcfishpass `r ver$model_version` (LINEAR model, +completed `r format(ver$date_completed, "%Y-%m-%d")`). + +The streams query is a single +[`fresh::frs_db_query()`](https://newgraphenvironment.github.io/fresh/reference/frs_db_query.html) +call — every row in `bcfishpass.streams_bt_vw` is already classified +accessible (`access IN (1, 2)` = assessed or modelled), so filtering to +"best accessible habitat order 3+" needs no working-table classification: + +```r +streams <- fresh::frs_db_query(conn, " + SELECT segmented_stream_id, blue_line_key, waterbody_key, + upstream_area_ha, map_upstream, channel_width, stream_order, + gradient, mapping_code, access, spawning, rearing, geom + FROM bcfishpass.streams_bt_vw + WHERE watershed_group_code = 'PARS' + AND access IN (1, 2) + AND stream_order >= 3") +``` + +Waterbodies (lakes + wetlands) come from `whse_basemapping.fwa_lakes_poly` +and `whse_basemapping.fwa_wetlands_poly` joined on `waterbody_key` — only +those physically anchored to the streams above are pulled in. They feed +`flooded::fl_valley_confine(waterbodies = ...)` so the final valley +raster fills lake / wetland cells the gradient and cost-distance masks +would otherwise carve donut holes around. + +## Valley confinement + +`flooded::fl_valley_confine()` runs the full VCA pipeline on the cached +DEM, streams, and waterbodies. On PARS at 30 m (~20 Mcells) it takes a +couple of minutes single-threaded; `terra::terraOptions(threads = N)` +parallelises the heavy raster ops. + +```{r vca, eval = FALSE} +valleys <- flooded::fl_valley_confine( + dem = dem, + streams = streams, + field = "upstream_area_ha", + precip = flooded::fl_stream_rasterize(streams, dem, field = "map_upstream"), + waterbodies = waterbodies, + flood_factor = 4 # ff04 — functional floodplain +) +floodplain <- flooded::fl_valley_poly(valleys) +``` + +## Floodplain map — full WSG + +```{r map-floodplain, echo = FALSE, fig.cap = "PARS unconfined valleys (green) over MRDEM-30 hillshade. Parks (light green polygon), First Nations reserves (light grey polygon, black diamond marker at centroid + label), accessible order 3+ streams in blue, lakes and wetlands in light blue, forest service / resource roads grey, railways black-dashed, watershed boundary heavy black."} +hill <- terra::shade( + terra::terrain(dem, "slope", unit = "radians"), + terra::terrain(dem, "aspect", unit = "radians") +) + +terra::plot(hill, col = grey(0:100 / 100), legend = FALSE, axes = FALSE, + mar = c(2, 2, 2, 2), + main = "PARS — functional floodplain (ff04), accessible order 3+") + +# Hydro / floodplain +terra::plot(valleys, col = c(NA, "#2c7a2c"), alpha = 0.55, + add = TRUE, legend = FALSE) +plot(sf::st_geometry(waterbodies), col = "#a3cdb966", border = "#2171B5", + lwd = 0.3, add = TRUE) +plot(sf::st_geometry(streams), col = "#1f77b4", lwd = 0.4, add = TRUE) + +# Transport — FSRs only on the wide view; full network too cluttered +if (!is.null(roads)) { + fsr <- roads[roads$transport_line_type_code %in% c("RRS", "RRD", "RRN"), ] + if (nrow(fsr) > 0L) { + plot(sf::st_geometry(fsr), col = "#88888888", lwd = 0.25, add = TRUE) + } +} +if (!is.null(railways)) { + plot(sf::st_geometry(railways), col = "black", lwd = 0.8, lty = "dashed", + add = TRUE) +} + +# Land-use polygons drawn ABOVE hydro/transport so they're visible where +# they overlap. Colours from the gq registry (`reg_qgis_restoration.json`). +if (!is.null(parks)) { + plot(sf::st_geometry(parks), col = "#639b5f55", border = "#33a02c", + lwd = 0.7, add = TRUE) +} +if (!is.null(reserves)) { + plot(sf::st_geometry(reserves), col = "#b2b2b288", border = "#232323", + lwd = 0.6, add = TRUE) + # Reserves are tiny at WSG scale — overlay a point marker at each + # centroid so they read on the map. + pts <- suppressWarnings(sf::st_centroid(sf::st_geometry(reserves), + of_largest_polygon = TRUE)) + plot(pts, pch = 23, bg = "#222222", col = "white", cex = 0.9, add = TRUE) +} + +plot(sf::st_geometry(aoi), border = "black", lwd = 1.5, add = TRUE) + +# Reserve labels — drawn last and with xpd = TRUE so a label whose +# centroid sits near the AOI edge after clipping isn't cropped at the +# plot region. Tiny font with a thin white halo for legibility on +# hillshade. +if (!is.null(reserves) && "english_name" %in% names(reserves) && + nrow(reserves) > 0L) { + pts <- suppressWarnings(sf::st_centroid(sf::st_geometry(reserves), + of_largest_polygon = TRUE)) + coords <- sf::st_coordinates(pts) + ofs_y <- par("cxy")[2] * 0.45 * 0.9 + halo_x <- par("cxy")[1] * 0.45 * 0.05 + halo_y <- par("cxy")[2] * 0.45 * 0.05 + ly <- coords[, 2] - ofs_y + for (dx in c(-1, 1)) { + text(coords[, 1] + dx * halo_x, ly, labels = reserves$english_name, + cex = 0.45, col = "white", font = 2, xpd = TRUE) + } + for (dy in c(-1, 1)) { + text(coords[, 1], ly + dy * halo_y, labels = reserves$english_name, + cex = 0.45, col = "white", font = 2, xpd = TRUE) + } + text(coords[, 1], ly, labels = reserves$english_name, + cex = 0.45, col = "#111111", font = 2, xpd = TRUE) +} +``` + +## Detail map — south-east corner + +The full-WSG view compresses a lot of detail. Cropping to the bottom-right +quadrant — the lower confluences where the trunk approaches Williston — +shows the per-reach floodplain pattern, individual channels, lakes / +wetlands, and the named First Nations reserves at full resolution. + +```{r map-detail, echo = FALSE, fig.cap = "South-east corner of PARS at full resolution. Parks (light green), First Nations reserves (light grey polygon with black diamond marker + formal english_name label at centroid), waterbodies, valleys, named streams (italic blue labels), roads (grey), railways (black dashed)."} +e <- sf::st_bbox(aoi) +inset_bbox <- sf::st_bbox(c( + xmin = unname(e["xmin"] + (e["xmax"] - e["xmin"]) * 0.55), + ymin = unname(e["ymin"]), + xmax = unname(e["xmax"]), + ymax = unname(e["ymin"] + (e["ymax"] - e["ymin"]) * 0.45) +), crs = sf::st_crs(streams)) +inset_ext <- terra::ext(inset_bbox[c("xmin", "xmax", "ymin", "ymax")]) + +crop_sf <- function(x) { + if (is.null(x) || nrow(x) == 0L) return(NULL) + out <- suppressWarnings(sf::st_crop(x, inset_bbox)) + if (nrow(out) == 0L) NULL else out +} + +dem_inset <- terra::crop(dem, inset_ext) +valleys_inset <- terra::crop(valleys, inset_ext) +streams_inset <- crop_sf(streams) +wb_inset <- crop_sf(waterbodies) +aoi_inset <- crop_sf(aoi) +roads_inset <- crop_sf(roads) +railways_inset <- crop_sf(railways) +reserves_inset <- crop_sf(reserves) +parks_inset <- crop_sf(parks) +named_streams_inset <- crop_sf(named_streams) + +hill_inset <- terra::shade( + terra::terrain(dem_inset, "slope", unit = "radians"), + terra::terrain(dem_inset, "aspect", unit = "radians") +) + +terra::plot(hill_inset, col = grey(0:100 / 100), legend = FALSE, axes = FALSE, + mar = c(2, 2, 2, 2), main = "PARS — south-east detail") + +# Hydro / floodplain +terra::plot(valleys_inset, col = c(NA, "#2c7a2c"), alpha = 0.55, + add = TRUE, legend = FALSE) +if (!is.null(wb_inset)) { + plot(sf::st_geometry(wb_inset), col = "#a3cdb988", border = "#2171B5", + lwd = 0.4, add = TRUE) +} +plot(sf::st_geometry(streams_inset), col = "#1f77b4", lwd = 0.7, add = TRUE) + +# Transport +if (!is.null(roads_inset)) { + plot(sf::st_geometry(roads_inset), col = "#666666aa", lwd = 0.4, add = TRUE) +} +if (!is.null(railways_inset)) { + plot(sf::st_geometry(railways_inset), col = "black", lwd = 1.0, lty = "dashed", + add = TRUE) +} + +# Land-use polygons on top so they read on overlap. +if (!is.null(parks_inset)) { + plot(sf::st_geometry(parks_inset), col = "#639b5f55", border = "#33a02c", + lwd = 0.8, add = TRUE) +} +if (!is.null(reserves_inset)) { + plot(sf::st_geometry(reserves_inset), col = "#b2b2b2aa", border = "#232323", + lwd = 0.7, add = TRUE) +} + +if (!is.null(aoi_inset)) { + plot(sf::st_geometry(aoi_inset), border = "black", lwd = 1.5, add = TRUE) +} + +# Named-stream labels — one per unique gnis_name at the segment midpoint +if (!is.null(named_streams_inset)) { + dedup <- named_streams_inset[!duplicated(named_streams_inset$gnis_name), ] + pts <- suppressWarnings(sf::st_centroid(sf::st_geometry(dedup))) + text(sf::st_coordinates(pts), labels = dedup$gnis_name, + cex = 0.5, font = 3, col = "#0d3a6c") +} + +# First Nations reserve markers + labels (inset only) +if (!is.null(reserves_inset) && "english_name" %in% names(reserves_inset) && + nrow(reserves_inset) > 0L) { + pts <- suppressWarnings(sf::st_centroid(sf::st_geometry(reserves_inset), + of_largest_polygon = TRUE)) + coords <- sf::st_coordinates(pts) + plot(pts, pch = 23, bg = "#222222", col = "white", cex = 1.0, add = TRUE) + ofs_y <- par("cxy")[2] * 0.45 * 0.9 + halo_x <- par("cxy")[1] * 0.45 * 0.05 + halo_y <- par("cxy")[2] * 0.45 * 0.05 + ly <- coords[, 2] - ofs_y + for (dx in c(-1, 1)) { + text(coords[, 1] + dx * halo_x, ly, labels = reserves_inset$english_name, + cex = 0.45, col = "white", font = 2, xpd = TRUE) + } + for (dy in c(-1, 1)) { + text(coords[, 1], ly + dy * halo_y, labels = reserves_inset$english_name, + cex = 0.45, col = "white", font = 2, xpd = TRUE) + } + text(coords[, 1], ly, labels = reserves_inset$english_name, + cex = 0.45, col = "#111111", font = 2, xpd = TRUE) +} +``` + +## References diff --git a/vignettes/references.bib b/vignettes/references.bib new file mode 100644 index 0000000..37591e1 --- /dev/null +++ b/vignettes/references.bib @@ -0,0 +1,26 @@ +@article{hall_etal2007Predictingriver, + title = {Predicting River Floodplain and Lateral Channel Migration for Salmon Habitat Conservation}, + author = {Hall, Jason E. and Holzer, Diane M. and Beechie, Timothy J.}, + date = {2007}, + journaltitle = {Journal of the American Water Resources Association}, + volume = {43}, + number = {3}, + pages = {786--797}, + doi = {10.1111/j.1752-1688.2007.00063.x}, + url = {https://doi.org/10.1111/j.1752-1688.2007.00063.x}, + abstract = {In this article, we describe a method for predicting floodplain locations and potential lateral channel migration across 82,900 km (491 km2 by bankfull area) of streams in the Columbia River basin. Predictions are based on channel confinement, channel slope, bankfull width, and bankfull depth derived from digital elevation and precipitation data. Half of the 367 km2 (47,900 km by length) of low-gradient channels (≤ 4\% channel slope) were classified as floodplain channels with a high likelihood of lateral channel migration (182 km2, 50\%). Classification agreement between modeled and field-measured floodplain confinement was 85\% (κ = 0.46, p {$<$} 0.001) with the largest source of error being the misclassification of unconfined channels as confined (55\% omission error). Classification agreement between predicted channel migration and lateral migration determined from aerial photographs was 76\% (κ = 0.53, p {$<$} 0.001) with the largest source of error being the misclassification of laterally migrating channels as non-migrating (35\% omission error). On average, more salmon populations were associated with laterally migrating channels and floodplains than with confined or nonmigrating channels. These data are useful for many river basin planning applications, including identification of land use impacts to floodplain habitats and locations with restoration potential for listed salmonids or other species of concern.}, + file = {/Users/airvine/Zotero/storage/W6LD4RRG/hall_et_al._2007-predicting_river_flo.pdf} +} + +@report{nagel_etal2014LandscapeScale, + title = {A {{Landscape Scale Valley Confinement Algorithm}}: {{Delineating Unconfined Valley Bottoms}} for {{Geomorphic}}, {{Aquatic}}, and {{Riparian Applications}}}, + author = {Nagel, David E. and Buffington, John M. and Parkes, Sharon L. and Wenger, Seth and Goode, John R.}, + date = {2014}, + number = {Gen. Tech. Rep. RMRS-GTR-321}, + pages = {12}, + institution = {U.S. Department of Agriculture, Forest Service, Rocky Mountain Research Station}, + location = {Fort Collins, CO}, + url = {https://www.fs.usda.gov/rmrs/publications/landscape-scale-valley-confinement-algorithm-delineating-unconfined-valley-bottoms}, + file = {/Users/airvine/Zotero/storage/TFBBPKGI/Nagel_et_al_2014_RMRS-GTR-321.pdf} +} +