Continuous index-trajectory change detection (Sentinel-2 + BFAST) (#30)

.Rbuildignore

@@ -11,3 +11,5 @@
 ^doc$
 ^Meta$
 ^CITATION\.cff$
+^planning$
+^\.claude$

.gitignore

@@ -7,3 +7,6 @@ docs
 /doc/
 /Meta/
 vignettes/*.html
+
+# local STAC cube cache (generated by data-raw/vignette_data_break.R)
+data-raw/.break_cache/

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: drift
 Title: Detecting Riparian and Inland Floodplain Transitions
-Version: 0.2.4
-Date: 2026-07-07
+Version: 0.3.0
+Date: 2026-07-08
 Authors@R: c(
     person("Allan", "Irvine", , "[email protected]", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0002-3495-2128")),
@@ -31,6 +31,7 @@ Imports:
     terra (>= 1.8-10),
     tibble
 Suggests:
+    bfast,
     bookdown,
     DT,
     flooded,

NAMESPACE

@@ -4,12 +4,17 @@ export(dft_cache_clear)
 export(dft_cache_info)
 export(dft_cache_path)
 export(dft_class_table)
+export(dft_index_expr)
+export(dft_index_table)
 export(dft_map_interactive)
+export(dft_rast_break)
 export(dft_rast_classify)
 export(dft_rast_consensus)
 export(dft_rast_summarize)
 export(dft_rast_transition)
+export(dft_rast_trend)
 export(dft_stac_classes)
 export(dft_stac_config)
+export(dft_stac_cube)
 export(dft_stac_fetch)
 export(dft_transition_vectors)

NEWS.md

@@ -1,3 +1,11 @@
+# drift 0.3.0
+
+- Continuous index-trajectory change detection for floodplain reaches (#30). A new fetch-and-reduce pipeline complements the categorical `dft_stac_fetch()` path. `dft_stac_cube()` builds a cloud-masked monthly spectral-index stack from Sentinel-2 (via a new `"sentinel-2-l2a"` source); `dft_rast_break()` reduces it per pixel with `bfast::bfastmonitor()` into a two-band raster of *abrupt* break date and magnitude; and `dft_rast_trend()` reduces it to a per-pixel *gradual* trend — a robust Theil-Sen slope (index change per year) with Mann-Kendall significance — for degradation/recovery monitoring the annual labels cannot show. Together they let a continuous trajectory validate categorical land-cover transitions (confirming which mapped losses carry a real spectral decline) and detect gradual change. See the "Trajectories as a Check on Land-Cover Change" vignette.
+- `dft_index_expr()` and `dft_index_table()` add a table-driven spectral-index registry (NDVI, kNDVI, NDMI) whose formulas are written over band *roles*, so one index resolves against any reflectance source; the reflectance scale/offset is folded into each expression.
+- Sentinel-2 handling is correctness-focused: `dft_stac_cube()` masks cloud/shadow/cirrus/snow, restricts to caller-chosen calendar `months` (e.g. the growing season) to sharpen the signal and cut scenes streamed, and — because the +1000 DN reflectance offset only applies from processing baseline 04.00 (2022-01-25) — splits items at that boundary and corrects each side, so a multi-year series carries no artificial index step at 2022.
+- `dft_stac_config()` gains a role-based schema for reflectance cube sources (band roles, mask classes, scale/offset, offset boundary), leaving the categorical `io-lulc`/`esa-worldcover` sources unchanged. `bfast` added to Suggests.
+- Known limitation tracked as a follow-up: the cube spans the AOI bounding box rather than the polygon (a gdalcubes `filter_geom` limitation, #32); labelling breaks with from/to land-cover classes is #31.
+
 # drift 0.2.4
 
 - `dft_transition_vectors()` no longer exhausts memory on large-extent rasters (#27). The per-class loop allocated full-grid vectors per class and per patch — ncell × n_patches churn that OOM-killed a 102.6M-cell, 56-class floodplain. Replaced by a single `terra::patches(values = TRUE)` pass plus a sparse patch-to-label map. Output is identical (verified patch-by-patch against the old implementation); only `patch_id` numbering / row order changes, to raster scan order. Benchmark at 24M cells: 1.9 s for a 4,799-patch raster; the old code took 122 s on a milder 1,232-patch raster of the same size.

R/dft_index_expr.R

@@ -0,0 +1,143 @@
+#' Apply a spectral index to a data cube
+#'
+#' Resolve a named spectral index (e.g. `"kndvi"`) into a per-pixel arithmetic
+#' expression over a source's band roles and apply it to a `gdalcubes` data
+#' cube, returning a single-band cube named after the index.
+#'
+#' Index formulas are stored in a shipped registry ([dft_index_table()]) written
+#' over band **roles** (`red`, `nir`, `swir16`), not literal asset names. The
+#' roles are resolved to per-source asset names via [dft_stac_config()], so the
+#' same `"kndvi"` works against Sentinel-2 (`B04`/`B08`) or any future
+#' reflectance source without changing the formula.
+#'
+#' Reflectance `scale`/`offset` are folded **into** the expression as a per-band
+#' affine transform `(asset * scale + offset)`. This matters for ratio indices:
+#' a non-zero offset does not cancel in `(nir - red)/(nir + red)`, so computing
+#' NDVI on raw digital numbers is wrong for sources with an offset (Landsat C2
+#' L2, or Sentinel-2 processing baseline 04.00).
+#'
+#' @param cube A `gdalcubes` data cube (e.g. the lazy cube built inside
+#'   [dft_stac_cube()]) whose bands are the source's assets.
+#' @param index Character. An index name present in [dft_index_table()]
+#'   (default `"kndvi"`).
+#' @param source Character. Source name passed to [dft_stac_config()] to resolve
+#'   the role→asset map and reflectance scale/offset (default
+#'   `"sentinel-2-l2a"`).
+#' @param roles Named list mapping roles to asset names. When `NULL`, taken from
+#'   `dft_stac_config(source)$roles`.
+#' @param scale,offset Numeric reflectance affine transform. When `NULL`, taken
+#'   from the source config (falling back to `1` / `0`).
+#'
+#' @return A single-band `gdalcubes` cube with the band named `index`.
+#'
+#' @seealso [dft_index_table()] for the registry, [dft_stac_cube()] for the
+#'   caller that builds the input cube.
+#'
+#' @examples
+#' # The registry the resolver reads:
+#' dft_index_table()
+#'
+#' \dontrun{
+#' # Applied to a lazy Sentinel-2 cube (requires network + gdalcubes):
+#' aoi <- sf::st_read(system.file("extdata", "example_aoi.gpkg", package = "drift"))
+#' cube <- dft_stac_cube(aoi, index = "kndvi")  # dft_stac_cube calls this internally
+#' }
+#'
+#' @export
+dft_index_expr <- function(cube,
+                           index = "kndvi",
+                           source = "sentinel-2-l2a",
+                           roles = NULL,
+                           scale = NULL,
+                           offset = NULL) {
+  rlang::check_installed("gdalcubes", reason = "to apply an index to a cube")
+  cfg <- dft_stac_config(source)
+  roles <- roles %||% cfg$roles
+  scale <- scale %||% cfg$scale %||% 1
+  offset <- offset %||% cfg$offset %||% 0
+  expr <- index_resolve_expr(index, roles, scale, offset)
+  gdalcubes::apply_pixel(cube, expr, names = index)
+}
+
+#' Load the shipped spectral-index registry
+#'
+#' Reads the CSV index registry bundled with the package. Each row defines one
+#' index as a `gdalcubes`/tinyexpr formula written over band roles.
+#'
+#' @return A tibble with columns `index`, `formula`, `roles` (comma-separated
+#'   role names), and `description`.
+#'
+#' @examples
+#' dft_index_table()
+#'
+#' @export
+dft_index_table <- function() {
+  path <- system.file("indices", "indices.csv", package = "drift", mustWork = TRUE)
+  tibble::as_tibble(utils::read.csv(path, stringsAsFactors = FALSE))
+}
+
+#' Look up one index registry row, erroring on unknown index
+#' @noRd
+index_row <- function(index) {
+  tbl <- dft_index_table()
+  row <- tbl[tbl$index == index, ]
+  if (nrow(row) == 0) {
+    cli::cli_abort(c(
+      "Unknown index {.val {index}}.",
+      "i" = "Available indices: {.val {tbl$index}}."
+    ))
+  }
+  row
+}
+
+#' Band roles required by an index (internal)
+#'
+#' Used by [dft_stac_cube()] to decide which assets to pull.
+#' @noRd
+index_roles <- function(index) {
+  trimws(strsplit(index_row(index)$roles[[1]], ",")[[1]])
+}
+
+#' Build a per-band affine reflectance token for an expression (internal)
+#'
+#' Returns the bare asset name when the transform is identity (`scale == 1`,
+#' `offset == 0`); otherwise `(asset * scale +/- |offset|)`. Numbers are
+#' formatted without scientific notation so the tinyexpr C parser accepts them.
+#' @noRd
+scale_token <- function(asset, scale, offset) {
+  if (scale == 1 && offset == 0) return(asset)
+  core <- if (scale == 1) {
+    asset
+  } else {
+    sprintf("%s * %s", asset, format(scale, scientific = FALSE, trim = TRUE))
+  }
+  if (offset == 0) return(sprintf("(%s)", core))
+  sign <- if (offset < 0) "-" else "+"
+  sprintf("(%s %s %s)", core, sign,
+          format(abs(offset), scientific = FALSE, trim = TRUE))
+}
+
+#' Resolve an index name to a per-pixel expression string (internal)
+#'
+#' Substitutes each role token in the registry formula with its scaled asset
+#' token. Roles are substituted longest-name-first with word boundaries so a
+#' shorter role name cannot clobber part of a longer one.
+#' @noRd
+index_resolve_expr <- function(index, roles, scale = 1, offset = 0) {
+  row <- index_row(index)
+  formula <- row$formula[[1]]
+  needed <- trimws(strsplit(row$roles[[1]], ",")[[1]])
+  needed <- needed[order(nchar(needed), decreasing = TRUE)]
+  for (role in needed) {
+    asset <- roles[[role]]
+    if (is.null(asset)) {
+      cli::cli_abort(c(
+        "Index {.val {index}} needs role {.val {role}}, absent from the role map.",
+        "i" = "Available roles: {.val {names(roles)}}."
+      ))
+    }
+    formula <- gsub(paste0("\\b", role, "\\b"),
+                    scale_token(asset, scale, offset), formula)
+  }
+  formula
+}

R/dft_rast_break.R

@@ -0,0 +1,153 @@
+#' Detect per-pixel index-trajectory breakpoints
+#'
+#' Reduce a monthly index stack (from [dft_stac_cube()]) over time with
+#' [bfast::bfastmonitor()], returning a two-band `SpatRaster` of the break date
+#' and magnitude for every pixel. Where categorical differencing compares two
+#' land-cover labels, this asks a stronger question of a continuous index
+#' trajectory: *when* did the pixel's spectral history break, and by how much?
+#'
+#' `bfastmonitor` fits a season-trend model to a stable *history* period, then
+#' watches the *monitoring* period (from `start` onward) for a structural break.
+#' The returned `break_mag` is the median monitoring-period residual: **negative
+#' means the index dropped** (e.g. vegetation loss / channel scour), positive
+#' means it rose (establishment). `break_date` is a decimal year (e.g.
+#' `2022.42`) or `NA` where no break was detected.
+#'
+#' Pixels are reduced in parallel with `parallel::mclapply()` (forked workers, so
+#' the per-pixel logic and its parameters are inherited directly). Pixels with
+#' fewer than `min_obs` valid observations short-circuit to `NA`.
+#'
+#' @param cube A monthly index `SpatRaster` (the return value of
+#'   [dft_stac_cube()]): one layer per time step, with a time value per layer.
+#' @param history Character. `bfastmonitor` history-selection method: `"all"`
+#'   (default), `"ROC"`, or `"BP"`.
+#' @param start Numeric `c(year, period)`. Start of the monitoring period, in the
+#'   stack's temporal frequency (e.g. `c(2022, 1)` = Jan 2022 for a monthly
+#'   stack). Everything before it is the stable history.
+#' @param frequency Numeric or `NULL`. Seasonal frequency of the time series
+#'   (12 for monthly, 1 for annual). When `NULL`, derived from the layer time
+#'   spacing; when supplied, it must agree with that spacing or the call errors.
+#' @param order Integer. Harmonic order of the season-trend model passed to
+#'   [bfast::bfastmonitor()] (default 3). Lower it (1-2) when the series samples
+#'   only part of the year (e.g. a growing-season-only cube from
+#'   [dft_stac_cube()] `months`), where a high order overfits sparse seasonal
+#'   coverage.
+#' @param level Numeric. Significance level passed to [bfast::bfastmonitor()]
+#'   (default 0.01).
+#' @param min_obs Integer. Minimum non-`NA` observations required to attempt a
+#'   fit; pixels with fewer return `NA` (default 6).
+#' @param cores Integer or `NULL`. Forked workers for the per-pixel reduction.
+#'   When `NULL`, uses one fewer than the detected cores.
+#'
+#' @return A two-band [terra::SpatRaster] with layers `break_date` (decimal year
+#'   or `NA`) and `break_mag` (signed index change; negative = index drop).
+#'
+#' @seealso [dft_stac_cube()] (builds the input stack), [dft_index_expr()].
+#'
+#' @examples
+#' \dontrun{
+#' # Requires network + gdalcubes + bfast
+#' aoi <- sf::st_read(system.file("extdata", "example_aoi.gpkg", package = "drift"))
+#' cube <- dft_stac_cube(aoi, index = "kndvi", datetime = "2019-01-01/2023-12-31")
+#' breaks <- dft_rast_break(cube, start = c(2022, 1))
+#' terra::plot(breaks[["break_mag"]])  # negative (blue) = scour / veg loss
+#' }
+#'
+#' @export
+dft_rast_break <- function(cube,
+                           history = "all",
+                           start = c(2022, 1),
+                           frequency = NULL,
+                           order = 3,
+                           level = 0.01,
+                           min_obs = 6,
+                           cores = NULL) {
+  rlang::check_installed("bfast", reason = "for trajectory breakpoint detection")
+  if (!inherits(cube, "SpatRaster")) {
+    cli::cli_abort("`cube` must be a SpatRaster time stack from {.fn dft_stac_cube}.")
+  }
+  tm <- terra::time(cube)
+  if (length(tm) != terra::nlyr(cube) || anyNA(tm)) {
+    cli::cli_abort(c(
+      "`cube` must carry a time value for every layer.",
+      "i" = "Pass the stack returned by {.fn dft_stac_cube}."
+    ))
+  }
+
+  # derive the ts() start and seasonal frequency from the layer times
+  cadence_freq <- cadence_frequency(tm)
+  if (is.na(cadence_freq)) {
+    cli::cli_abort("Unsupported layer cadence; use a monthly or annual stack.")
+  }
+  if (is.null(frequency)) {
+    frequency <- cadence_freq
+  } else if (!isTRUE(all.equal(as.numeric(frequency), as.numeric(cadence_freq)))) {
+    cli::cli_abort(c(
+      "`frequency` ({frequency}) disagrees with the layer cadence (= {cadence_freq}).",
+      "i" = "Leave `frequency = NULL` to derive it from the stack."
+    ))
+  }
+  t0 <- as.Date(tm[1])
+  yr <- as.integer(format(t0, "%Y"))
+  mo <- as.integer(format(t0, "%m"))
+  ts_start <- c(yr, floor((mo - 1) / (12 / frequency)) + 1)
+
+  if (is.null(cores)) {
+    dc <- parallel::detectCores()
+    cores <- if (is.na(dc)) 2L else max(1L, dc - 1L)
+  }
+
+  # reduce only pixels with a usable series; the rest stay NA
+  vals <- terra::values(cube)
+  usable <- which(rowSums(!is.na(vals)) >= min_obs)
+  res <- matrix(NA_real_, nrow(vals), 2)
+  if (length(usable)) {
+    chunks <- split(usable, (seq_along(usable) - 1) %% cores)
+    parts <- parallel::mclapply(chunks, function(ii) {
+      t(vapply(ii, function(i) {
+        .dft_break_pixel(vals[i, ], ts_start, frequency, start, history, order,
+                         level, min_obs)
+      }, numeric(2)))
+    }, mc.cores = cores)
+    for (k in seq_along(chunks)) res[chunks[[k]], ] <- parts[[k]]
+  }
+
+  out <- cube[[1:2]]
+  terra::values(out) <- res
+  names(out) <- c("break_date", "break_mag")
+  out
+}
+
+
+#' Per-pixel breakpoint reducer logic (internal, unit-testable)
+#'
+#' The degenerate branches (all-`NA` or fewer than `min_obs` observations) return
+#' `c(NA, NA)` before any `bfast` symbol is touched, so they are testable without
+#' bfast installed.
+#' @noRd
+.dft_break_pixel <- function(v, ts_start, frequency, start, history, order,
+                             level, min_obs) {
+  if (all(is.na(v)) || sum(!is.na(v)) < min_obs) return(c(NA_real_, NA_real_))
+  ts_v <- stats::ts(v, start = ts_start, frequency = frequency)
+  tryCatch({
+    m <- bfast::bfastmonitor(ts_v, start = start, history = history,
+                             order = order, level = level)
+    c(m$breakpoint, m$magnitude)
+  }, error = function(e) c(NA_real_, NA_real_))
+}
+
+
+#' Seasonal frequency implied by a stack's layer times (internal)
+#'
+#' Monthly spacing -> 12, quarterly -> 4, annual -> 1. Returns `NA` for
+#' unsupported cadences.
+#' @noRd
+cadence_frequency <- function(tm) {
+  if (length(tm) < 2) return(NA_real_)
+  d <- stats::median(as.numeric(diff(as.Date(tm))))
+  if (is.na(d)) return(NA_real_)
+  if (d >= 26 && d <= 32) return(12)
+  if (d >= 85 && d <= 95) return(4)
+  if (d >= 360 && d <= 370) return(1)
+  NA_real_
+}