Add physical-range attribute bounds check to invariant base layer

Adds a per-attribute physical-plausibility pass complementing the existing
topology/geometry invariants:

- hf_attr_bounds(): curated bounds table (inst/extdata/divide_attr_bounds.csv)
  giving the physically reasonable range for each divide/flowpath modeling
  parameter. A `caveat` column flags bounds that are too strict (strict-zero
  lowers that reject valid flat/zero values) or look like placeholders; those
  are excluded from checks by default. An `aliases` column (;-separated) lets a
  single bound match alternate column names, so one row covers schema naming
  differences (e.g. lengthkm matches length_km).
- hf_check_attr_bounds(): standalone validator; absent columns are skipped,
  canonical/alias names matched in order so the populated column wins.
- hf_check_invariants() gains an opt-in `attr_bounds=`/`domain=` hook. The pass
  is soft (warn-only) and off by default, so it never changes the returned `ok`
  for existing callers until explicitly enabled.

Bounds are aligned to our schema: flowpath length->lengthkm, incremental
area->divide_areasqkm, drainage->total_dasqkm, divide area->areasqkm, with the
source names retained as aliases.

Tests cover range detection, inclusive boundaries, caveat exclusion, alias
matching, and soft-under-strict behavior.

Author: Mike Johnson

DESCRIPTION

@@ -61,6 +61,7 @@ Collate:
     'duckdb.R'
     'geom.R'
     'gpkg_metadata.R'
+    'hf_attr_bounds.R'
     'hf_invariants.R'
     'hfio.R'
     'hfutils-package.R'

NAMESPACE

@@ -19,6 +19,8 @@ export(get_hydroseq)
 export(get_node)
 export(gpkg_get_version)
 export(gpkg_set_version)
+export(hf_attr_bounds)
+export(hf_check_attr_bounds)
 export(hf_check_invariants)
 export(layer_exists)
 export(lynker_spatial_auth)

R/hf_attr_bounds.R

@@ -0,0 +1,111 @@
+#' Physical-range bounds for divide and flowpath attributes
+#'
+#' Returns the attribute plausibility-bounds table shipped with the package
+#' (`inst/extdata/divide_attr_bounds.csv`). Each row gives the physically
+#' reasonable range for a divide or flowpath modeling parameter. The `caveat`
+#' column flags bounds that are demonstrably too strict or look like
+#' placeholders; those are excluded from checks by default. The `aliases`
+#' column (`;`-separated) lists alternate column names a bound also matches, so
+#' a single bound covers schema naming differences (e.g. `lengthkm` matches
+#' `length_km`).
+#'
+#' @return A data.frame of bounds (one row per attribute, plus per-domain rows
+#'   for `lat`/`lon`).
+#' @export
+hf_attr_bounds <- function() {
+  path <- system.file("extdata", "divide_attr_bounds.csv", package = "hfutils")
+  if (!nzchar(path)) stop("divide_attr_bounds.csv not found in installed package", call. = FALSE)
+  utils::read.csv(path, stringsAsFactors = FALSE, na.strings = c("", "NA"))
+}
+
+# Build a named list of bound checks for a single layer. Returns results shaped
+# like .hf_ok() but marked soft = TRUE so they warn rather than abort even under
+# strict reporting. `prefix` disambiguates names when divides and flowpaths
+# checks are concatenated (both carry e.g. `area_sqkm`).
+.hf_attr_bounds_checks <- function(layer, which = c("divides", "flowpaths"),
+                                   domain = NULL, trust_caveated = FALSE,
+                                   prefix = NULL, max_report = 5L) {
+  which <- match.arg(which)
+  if (is.null(layer)) return(list())
+  df <- as.data.frame(layer)
+  b  <- hf_attr_bounds()
+  b  <- b[b$layer == which, , drop = FALSE]
+
+  is_domain_row <- b$domain != "all"
+  if (is.null(domain)) {
+    b <- b[!is_domain_row, , drop = FALSE]
+  } else {
+    b <- b[!is_domain_row | b$domain == domain, , drop = FALSE]
+  }
+  if (!trust_caveated) b <- b[is.na(b$caveat), , drop = FALSE]
+
+  has_aliases <- "aliases" %in% names(b)
+  checks <- list()
+  for (i in seq_len(nrow(b))) {
+    a <- b$attribute[i]
+    # Match our canonical name first, then any alias (handles schema naming
+    # differences, e.g. length_km vs lengthkm, and flowline_* variants).
+    cand <- a
+    if (has_aliases && !is.na(b$aliases[i]) && nzchar(b$aliases[i]))
+      cand <- c(cand, trimws(strsplit(b$aliases[i], ";", fixed = TRUE)[[1]]))
+    col <- cand[cand %in% names(df)]
+    if (!length(col)) next
+    col <- col[1]
+    v <- suppressWarnings(as.numeric(df[[col]]))
+    v <- v[!is.na(v)]
+    if (!length(v)) next
+
+    lo <- b$lower[i]; hi <- b$upper[i]
+    bad <- logical(length(v))
+    if (!is.na(lo)) bad <- bad | (if (isTRUE(b$lower_inclusive[i])) v < lo else v <= lo)
+    if (!is.na(hi)) bad <- bad | (if (isTRUE(b$upper_inclusive[i])) v > hi else v >= hi)
+
+    nm  <- if (is.null(prefix)) a else paste0(prefix, ".", a)
+    rng <- sprintf("[%s, %s]",
+                   if (is.na(lo)) "-Inf" else lo,
+                   if (is.na(hi)) "Inf" else hi)
+    if (sum(bad) == 0L) {
+      checks[[nm]] <- list(ok = TRUE, msg = sprintf("in range %s", rng),
+                           kind = "check", soft = TRUE)
+    } else {
+      ex <- utils::head(sort(unique(v[bad])), max_report)
+      checks[[nm]] <- list(
+        ok = FALSE, kind = "check", soft = TRUE,
+        msg = sprintf("%d/%d out of range %s (e.g. %s)",
+                      sum(bad), length(v), rng, paste(signif(ex, 4), collapse = ", ")))
+    }
+  }
+  checks
+}
+
+#' Check divide/flowpath attribute values against physical-range bounds
+#'
+#' Complements the topology/geometry checks in [hf_check_invariants()] with a
+#' per-attribute physical-plausibility pass: every modeling parameter present
+#' in the layer is checked against its expected range. Attributes not present
+#' in the layer are skipped; bounds carrying a `caveat` are skipped unless
+#' `trust_caveated = TRUE`. Failures are reported but never abort, even when
+#' `strict = TRUE`.
+#'
+#' @param layer A data.frame or sf object (the Divides or Flowpaths layer).
+#' @param which One of `"divides"` or `"flowpaths"`.
+#' @param domain Domain code for `lat`/`lon` bounds: one of `"CONUS"`, `"AK"`,
+#'   `"HI"`, `"PRVI"`, or `NULL` to skip lat/lon.
+#' @param strict Logical. Passed through to the reporter; attribute-bound
+#'   failures are soft and warn regardless.
+#' @param trust_caveated Logical. Include bounds flagged in the `caveat` column
+#'   (default `FALSE`).
+#' @param max_report Integer. Max example offending values to show per attribute.
+#'
+#' @return Invisibly a list with `ok` and `checks` (same shape as
+#'   [hf_check_invariants()]).
+#' @export
+hf_check_attr_bounds <- function(layer, which = c("divides", "flowpaths"),
+                                 domain = NULL, strict = FALSE,
+                                 trust_caveated = FALSE, max_report = 5L) {
+  which  <- match.arg(which)
+  checks <- .hf_attr_bounds_checks(layer, which, domain = domain,
+                                   trust_caveated = trust_caveated,
+                                   max_report = max_report)
+  .hf_report_checks(checks, stage = paste0("attr_bounds:", which), strict = strict)
+}

R/hf_invariants.R

@@ -15,6 +15,15 @@
 #'   coverage check (aggregated and ngen stages only). Default `FALSE`.
 #' @param coverage_min Minimum fraction of a flowpath that must lie inside
 #'   its assigned catchment. Default `0.90`.
+#' @param attr_bounds Logical. If `TRUE`, append a per-attribute physical-range
+#'   plausibility pass over the `divides` and `flowpaths` layers (see
+#'   [hf_check_attr_bounds()]). Soft/warn-only and off by default, so it never
+#'   changes the returned `ok` for existing callers until enabled. Absent
+#'   attribute columns are skipped.
+#' @param domain Domain code (`"CONUS"`, `"AK"`, `"HI"`, `"PRVI"`) selecting the
+#'   `lat`/`lon` bounds for the attribute pass; `NULL` skips lat/lon.
+#' @param attr_trust_caveated Logical. Include attribute bounds flagged in the
+#'   `caveat` column of the bounds table (default `FALSE`).
 #'
 #' @details
 #' Expected arguments per stage:
@@ -39,7 +48,9 @@
 #' @importFrom stats setNames
 #' @export
 hf_check_invariants <- function(stage, ..., strict = TRUE,
-                                coverage = FALSE, coverage_min = 0.90) {
+                                coverage = FALSE, coverage_min = 0.90,
+                                attr_bounds = FALSE, domain = NULL,
+                                attr_trust_caveated = FALSE) {
   stage <- match.arg(stage, c("refactored", "reconciled", "aggregated", "ngen"))
   args <- list(...)
   checks <- switch(stage,
@@ -50,6 +61,18 @@ hf_check_invariants <- function(stage, ..., strict = TRUE,
                                       coverage_min = coverage_min),
     ngen       = .hf_check_ngen(args, coverage = coverage,
                                 coverage_min = coverage_min))
+
+  # Optional per-attribute physical-plausibility pass. Soft (warn-only) and
+  # off by default, so it never disturbs callers gating on the returned `ok`
+  # until explicitly enabled. Absent attribute columns are skipped.
+  if (isTRUE(attr_bounds) && stage %in% c("reconciled", "aggregated", "ngen")) {
+    checks <- c(checks,
+      .hf_attr_bounds_checks(args$divides, "divides", domain = domain,
+                             trust_caveated = attr_trust_caveated, prefix = "div"),
+      .hf_attr_bounds_checks(args$flowpaths, "flowpaths", domain = domain,
+                             trust_caveated = attr_trust_caveated, prefix = "fp"))
+  }
+
   .hf_report_checks(checks, stage, strict)
 }
 
@@ -355,7 +378,7 @@ hf_check_invariants <- function(stage, ..., strict = TRUE,
     tag <- sprintf("[invariants:%s]", stage)
     line <- sprintf("%s %s %-32s %s", tag, icon, nm, res$msg)
     if (isTRUE(res$ok) || identical(res$kind, "info")) message(line)
-    else if (strict) stop(line, call. = FALSE)
+    else if (strict && !isTRUE(res$soft)) stop(line, call. = FALSE)
     else message(line)
   }
   invisible(list(ok = ok_all, stage = stage, checks = checks))

_pkgdown.yml

@@ -126,6 +126,8 @@ reference:
   - gpkg_set_version
   - gpkg_get_version
   - hf_check_invariants
+  - hf_check_attr_bounds
+  - hf_attr_bounds
 
 - title: DBI / OGR-SQL backend
   desc: >

inst/extdata/divide_attr_bounds.csv

@@ -0,0 +1,67 @@
+layer,domain,attribute,lower,upper,lower_inclusive,upper_inclusive,description,caveat,aliases
+divides,all,areasqkm,0,,FALSE,FALSE,catchment area (sqkm),,area_sqkm;divide_areasqkm
+divides,all,bexp_mode,2,15,FALSE,FALSE,Clapp-Hornberger beta exponent,,
+divides,all,isltyp_mode,1,16,TRUE,TRUE,dominant soil type category,,
+divides,all,ivgtyp_mode,1,16,TRUE,TRUE,dominant vegetation type category,,
+divides,all,dksat_geomean,1.95e-07,1.41e-03,FALSE,FALSE,saturated hydraulic conductivity,,
+divides,all,psisat_geomean,0.036,0.955,FALSE,FALSE,saturated capillary head,,
+divides,all,cwpvt_mean,0.09,0.36,FALSE,FALSE,empirical wind canopy parameter,,
+divides,all,mp_mean,3.6,12.6,FALSE,FALSE,slope of conductance-to-photosynthesis,,
+divides,all,mfsno_mean,0.5,4.0,FALSE,FALSE,melt factor for snow depletion curve,,
+divides,all,quartz_mean,0,1.0,FALSE,FALSE,mean soil quartz content,LOWER_TOO_STRICT: a legitimately zero-quartz divide fails gt 0; use ge 0,
+divides,all,refkdt_mean,0.1,4.0,FALSE,FALSE,reference soil infiltration parameter,,
+divides,all,slope1km_mean,0,1.0,FALSE,FALSE,hydraulic-head gradient at soil bottom,LOWER_TOO_STRICT: flat divide fails gt 0; use ge 0,
+divides,all,smcmax_mean,0.16,0.9,FALSE,FALSE,saturated soil moisture content,,
+divides,all,smcwlt_mean,0.05,0.3,FALSE,FALSE,wilting-point soil moisture content,,
+divides,all,vcmx_mean,24.0,112.0,FALSE,FALSE,max carboxylation rate,,
+divides,all,imperv_mean,0,1.0,FALSE,FALSE,impervious fraction,LOWER_TOO_STRICT: 0% impervious is valid and fails gt 0; use ge 0,
+divides,all,elevation_mean,-86.0,4422.0,FALSE,FALSE,mean terrain elevation (m),,
+divides,all,slope250m_mean,0,90.0,FALSE,FALSE,terrain slope (deg),LOWER_TOO_STRICT: flat terrain fails gt 0; use ge 0,
+divides,all,aspect_circmean,0,360.0,FALSE,FALSE,circular-mean terrain aspect (deg),LOWER_TOO_STRICT: due-north aspect=0 fails gt 0; use ge 0 / wrap,
+divides,all,lzfpm_mean,40.0,600.0,FALSE,FALSE,SAC max lower-zone free water primary,,
+divides,all,lzpk_mean,0.001,0.015,FALSE,FALSE,SAC lower-zone recession primary,,
+divides,all,lztwm_mean,75.0,300.0,FALSE,FALSE,SAC max lower-zone tension water,,
+divides,all,rexp_mean,1.4,3.5,FALSE,FALSE,SAC percolation exponent,,
+divides,all,uzk_mean,0.2,0.5,FALSE,FALSE,SAC upper-zone recession,,
+divides,all,zperc_mean,0,360.0,FALSE,FALSE,SAC min percolation rate coefficient,UPPER_SUSPECT: 360 ceiling looks like a placeholder, not a real limit,
+divides,all,lzfsm_mean,0,360.0,FALSE,FALSE,SAC max lower-zone free water secondary,UPPER_SUSPECT: 360 ceiling looks like a placeholder, not a real limit,
+divides,all,lzsk_mean,0.03,0.2,FALSE,FALSE,SAC lower-zone recession secondary,,
+divides,all,pfree_mean,0,0.5,FALSE,FALSE,SAC fraction percolating direct to LZ free water,,
+divides,all,uzfwm_mean,10.0,100.0,FALSE,FALSE,SAC max upper-zone free water,,
+divides,all,uztwm_mean,25.0,125.0,FALSE,FALSE,SAC max upper-zone tension water,,
+divides,all,mfmin_mean,0.01,0.6,FALSE,FALSE,Snow17 min non-rain melt factor,,
+divides,all,mfmax_mean,0,360.0,FALSE,FALSE,Snow17 max non-rain melt factor,UPPER_SUSPECT: 360 ceiling implausible; realistic mfmax is single digits,
+divides,all,uadj_mean,0.01,0.2,FALSE,FALSE,Snow17 avg wind function rain-on-snow,,
+divides,all,a_xinanjiang_inflection_point_parameter,-0.5,0.5,FALSE,FALSE,Xinanjiang inflection point,,
+divides,all,b_xinanjiang_shape_parameter,0.01,10.0,FALSE,FALSE,Xinanjiang shape b,,
+divides,all,x_xinanjiang_shape_parameter,0.01,10.0,FALSE,FALSE,Xinanjiang shape x,,
+divides,all,temp_delta_jan_mean,0,,FALSE,FALSE,mean temp difference January,,
+divides,all,temp_delta_feb_mean,0,,FALSE,FALSE,mean temp difference February,,
+divides,all,temp_delta_mar_mean,0,,FALSE,FALSE,mean temp difference March,,
+divides,all,temp_delta_apr_mean,0,,FALSE,FALSE,mean temp difference April,,
+divides,all,temp_delta_may_mean,0,,FALSE,FALSE,mean temp difference May,,
+divides,all,temp_delta_jun_mean,0,,FALSE,FALSE,mean temp difference June,,
+divides,all,temp_delta_jul_mean,0,,FALSE,FALSE,mean temp difference July,,
+divides,all,temp_delta_aug_mean,0,,FALSE,FALSE,mean temp difference August,,
+divides,all,temp_delta_sep_mean,0,,FALSE,FALSE,mean temp difference September,,
+divides,all,temp_delta_oct_mean,0,,FALSE,FALSE,mean temp difference October,,
+divides,all,temp_delta_nov_mean,0,,FALSE,FALSE,mean temp difference November,,
+divides,all,temp_delta_dec_mean,0,,FALSE,FALSE,mean temp difference December,,
+divides,all,glacier_percent,0,1,TRUE,TRUE,glacier cover fraction,,
+divides,all,cgw,1.80e-06,0.0018,FALSE,FALSE,groundwater coefficient,,
+divides,all,expon,1.0,8.0,FALSE,FALSE,groundwater exponent,,
+divides,all,max_gw_storage,0.01,0.25,FALSE,FALSE,baseflow bucket height,,
+divides,CONUS,lat,24,55,TRUE,TRUE,divide centroid latitude (CONUS),,
+divides,CONUS,lon,-125,-66,TRUE,TRUE,divide centroid longitude (CONUS),,
+divides,AK,lat,51,72,TRUE,TRUE,divide centroid latitude (AK),,
+divides,AK,lon,-172,-125,TRUE,TRUE,divide centroid longitude (AK),,
+divides,HI,lat,18,23,TRUE,TRUE,divide centroid latitude (HI),,
+divides,HI,lon,-161,-155,TRUE,TRUE,divide centroid longitude (HI),,
+divides,PRVI,lat,17,19,TRUE,TRUE,divide centroid latitude (PRVI),,
+divides,PRVI,lon,-68,-64,TRUE,TRUE,divide centroid longitude (PRVI),,
+flowpaths,all,lengthkm,0,,FALSE,FALSE,flowpath length (km),,length_km;flowline_lengthkm
+flowpaths,all,divide_areasqkm,0,,FALSE,FALSE,incremental divide area (sqkm),,areasqkm;area_sqkm;flowline_areasqkm
+flowpaths,all,total_dasqkm,0,,FALSE,FALSE,total upstream drainage area (sqkm),,total_da_sqkm;flowline_total_dasqkm
+flowpaths,all,path_length,0,,FALSE,FALSE,downstream path length,,
+flowpaths,all,mean_elevation,-86.0,4422.0,FALSE,FALSE,mean terrain elevation (m),,
+flowpaths,all,slope,0,90.0,FALSE,FALSE,terrain slope (deg),LOWER_TOO_STRICT: flat reach fails gt 0; use ge 0,

man/hf_attr_bounds.Rd

@@ -169,3 +169,35 @@ test_that("ngen stage detects a cycle in the fp->nexus->fp graph", {
     "cycle"
   )
 })
+
+test_that("hf_check_attr_bounds flags out-of-range and skips absent/zero-valid", {
+  df <- data.frame(
+    smcmax_mean     = c(0.4, 0.95, 0.2),   # 0.95 > 0.9 -> 1 bad
+    bexp_mode       = c(5, 5, 20),         # 20 > 15 -> 1 bad
+    glacier_percent = c(0, 0.5, 1),        # inclusive 0..1 -> all ok
+    area_sqkm       = c(1, 2, 3),          # all ok
+    not_a_bound     = c(9, 9, 9)           # no bound -> skipped
+  )
+  res <- suppressMessages(hf_check_attr_bounds(df, "divides", strict = FALSE))
+  expect_false(res$ok)
+  expect_false(res$checks[["smcmax_mean"]]$ok)
+  expect_false(res$checks[["bexp_mode"]]$ok)
+  expect_true(res$checks[["glacier_percent"]]$ok)
+  expect_true(res$checks[["areasqkm"]]$ok)  # data col `area_sqkm` matched via alias, reported canonically
+  expect_null(res$checks[["not_a_bound"]])
+})
+
+test_that("caveated bounds excluded by default; attr failures are soft under strict", {
+  df <- data.frame(imperv_mean = c(0, 0.2, 0.5), smcmax_mean = c(0.5, 0.95, 0.5))
+  res <- suppressMessages(hf_check_attr_bounds(df, "divides", strict = TRUE))
+  expect_null(res$checks[["imperv_mean"]])      # caveated -> skipped
+  expect_false(res$checks[["smcmax_mean"]]$ok)  # soft fail did not abort
+})
+
+test_that("hf_check_attr_bounds matches via aliases (schema name differences)", {
+  # bound canonical name is `lengthkm`; data uses the source alias `length_km`
+  df <- data.frame(length_km = c(0.5, -1, 2))   # -1 <= 0 -> 1 bad
+  res <- suppressMessages(hf_check_attr_bounds(df, "flowpaths", strict = FALSE))
+  expect_false(res$ok)
+  expect_false(res$checks[["lengthkm"]]$ok)     # matched via alias, reported canonically
+})