diff --git a/.Rbuildignore b/.Rbuildignore index d4909df..159d74f 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -11,3 +11,5 @@ ^doc$ ^Meta$ ^CITATION\.cff$ +^planning$ +^\.claude$ diff --git a/.gitignore b/.gitignore index fbc5e6d..e3e6e2c 100644 --- a/.gitignore +++ b/.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/ diff --git a/DESCRIPTION b/DESCRIPTION index 449a7da..b2af226 100644 --- a/DESCRIPTION +++ b/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", , "al@newgraphenvironment.com", 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, diff --git a/NAMESPACE b/NAMESPACE index 7a39c4a..a9fa599 100644 --- a/NAMESPACE +++ b/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) diff --git a/NEWS.md b/NEWS.md index a07e71d..a834e27 100644 --- a/NEWS.md +++ b/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. diff --git a/R/dft_index_expr.R b/R/dft_index_expr.R new file mode 100644 index 0000000..f117fca --- /dev/null +++ b/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 +} diff --git a/R/dft_rast_break.R b/R/dft_rast_break.R new file mode 100644 index 0000000..530cd02 --- /dev/null +++ b/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_ +} diff --git a/R/dft_rast_trend.R b/R/dft_rast_trend.R new file mode 100644 index 0000000..64c6574 --- /dev/null +++ b/R/dft_rast_trend.R @@ -0,0 +1,102 @@ +#' Per-pixel index trend over time +#' +#' Reduce an index stack (from [dft_stac_cube()]) to a per-pixel **trend** — the +#' rate at which the index rises or falls across the whole record. Where +#' [dft_rast_break()] finds an abrupt structural break, this measures the slow, +#' monotonic direction: gradual degradation (a declining stand) or recovery (a +#' restoration planting greening up) that annual categorical labels cannot show. +#' +#' The slope is a **Theil-Sen** estimate — the median of all pairwise slopes — +#' which is resistant to a single anomalous season (e.g. a smoke/drought year), +#' unlike an ordinary least-squares fit. Significance is the non-parametric +#' **Mann-Kendall** test for a monotonic trend. +#' +#' @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 min_obs Integer. Minimum non-`NA` observations required; 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 `trend` (index change per +#' year — negative = declining greenness, positive = recovering) and `trend_p` +#' (Mann-Kendall two-sided p-value; small = a significant monotonic trend). +#' +#' @seealso [dft_rast_break()] (abrupt breaks), [dft_stac_cube()]. +#' +#' @examples +#' \dontrun{ +#' aoi <- sf::st_read(system.file("extdata", "example_aoi.gpkg", package = "drift")) +#' cube <- dft_stac_cube(aoi, index = "kndvi", datetime = "2017-01-01/2023-12-31") +#' trend <- dft_rast_trend(cube) +#' terra::plot(trend[["trend"]]) # negative (red) = declining, positive = recovering +#' } +#' +#' @export +dft_rast_trend <- function(cube, min_obs = 6, cores = NULL) { + if (!inherits(cube, "SpatRaster")) { + cli::cli_abort("`cube` must be a SpatRaster time stack from {.fn dft_stac_cube}.") + } + if (terra::nlyr(cube) < 2) { + cli::cli_abort("`cube` needs at least 2 time layers to estimate a trend.") + } + 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}." + )) + } + min_obs <- max(2L, as.integer(min_obs)) # a slope needs at least two points + t_yr <- as.numeric(as.Date(tm)) / 365.25 # decimal years -> slope is per year + + if (is.null(cores)) { + dc <- parallel::detectCores() + cores <- if (is.na(dc)) 2L else max(1L, dc - 1L) + } + + 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_trend_pixel(vals[i, ], t_yr, 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("trend", "trend_p") + out +} + + +#' Theil-Sen slope + Mann-Kendall p for one pixel time series (internal) +#' +#' Robust to a minority of anomalous points. Returns `c(NA, NA)` on fewer than +#' `min_obs` valid observations, before any statistics — testable without data. +#' @noRd +.dft_trend_pixel <- function(y, t, min_obs) { + ok <- !is.na(y) + if (sum(ok) < min_obs) return(c(NA_real_, NA_real_)) + y <- y[ok] + t <- t[ok] + n <- length(y) + + ij <- utils::combn(n, 2) + i <- ij[1, ] + j <- ij[2, ] + # Theil-Sen: median pairwise slope + slope <- stats::median((y[j] - y[i]) / (t[j] - t[i])) + + # Mann-Kendall S = sum of signs of pairwise differences (ordered by time) + s <- sum(sign(y[j] - y[i])) + var_s <- n * (n - 1) * (2 * n + 5) / 18 + z <- if (s > 0) (s - 1) / sqrt(var_s) else if (s < 0) (s + 1) / sqrt(var_s) else 0 + p <- 2 * stats::pnorm(-abs(z)) + + c(slope, p) +} diff --git a/R/dft_stac_config.R b/R/dft_stac_config.R index 4fdc261..dbc110d 100644 --- a/R/dft_stac_config.R +++ b/R/dft_stac_config.R @@ -1,26 +1,39 @@ -#' Get STAC configuration for a known land cover source +#' Get STAC configuration for a known source #' #' Returns connection details for pre-configured STAC collections. -#' Used as a convenience wrapper around [dft_stac_fetch()] so users -#' don't need to remember STAC URLs and collection IDs. +#' Used as a convenience wrapper around [dft_stac_fetch()] (categorical +#' sources) and [dft_stac_cube()] (continuous index-trajectory sources) so +#' users don't need to remember STAC URLs, collection IDs, and band names. #' -#' @param source Character. One of `"io-lulc"` (Esri IO LULC annual v02) -#' or `"esa-worldcover"` (ESA WorldCover). +#' Sources are of two kinds. **Categorical** sources (`"io-lulc"`, +#' `"esa-worldcover"`) host single-band classified rasters and carry a flat +#' `asset` name for [dft_stac_fetch()]. **Cube** sources (`"sentinel-2-l2a"`) +#' host multi-band reflectance imagery and instead carry a role-based band map +#' (`red`/`nir`/`swir16`/`mask`), mask values, and reflectance scale/offset for +#' [dft_stac_cube()]; they are marked with `cube = TRUE`. The role-based schema +#' means a new reflectance source (e.g. Landsat C2 L2) drops in with no API +#' change — only the role→asset map and scale/offset differ. #' -#' @return A list with elements: -#' \describe{ -#' \item{stac_url}{STAC API endpoint} -#' \item{collection}{Collection ID} -#' \item{asset}{Asset name to download} -#' \item{available_years}{Integer vector of available years} -#' } +#' @param source Character. One of `"io-lulc"` (Esri IO LULC annual v02), +#' `"esa-worldcover"` (ESA WorldCover), or `"sentinel-2-l2a"` (Sentinel-2 +#' L2A surface reflectance, a cube source). +#' +#' @return A list. Categorical sources have elements `stac_url`, `collection`, +#' `asset`, `available_years`. Cube sources have `stac_url`, `collection`, +#' `cube = TRUE`, `roles` (a named list mapping `red`/`nir`/`swir16`/`mask` +#' to asset names), `mask_values` (integer mask classes to exclude), +#' `scale`/`offset` (DN → reflectance affine transform), and +#' `available_datetime` (an ISO 8601 interval string). The `cube` field is +#' absent (not `FALSE`) for categorical sources; test with `isTRUE(cfg$cube)`. #' #' @examples #' dft_stac_config("io-lulc") #' dft_stac_config("esa-worldcover") +#' dft_stac_config("sentinel-2-l2a") #' #' @export -dft_stac_config <- function(source = c("io-lulc", "esa-worldcover")) { +dft_stac_config <- function(source = c("io-lulc", "esa-worldcover", + "sentinel-2-l2a")) { source <- match.arg(source) switch(source, "io-lulc" = list( @@ -34,6 +47,46 @@ dft_stac_config <- function(source = c("io-lulc", "esa-worldcover")) { collection = "esa-worldcover", asset = "map", available_years = c(2020L, 2021L) + ), + "sentinel-2-l2a" = list( + # Cost: free — Microsoft Planetary Computer (Azure Open Data), no egress + # charge. Check egress terms before adding a requester-pays source. + stac_url = "https://planetarycomputer.microsoft.com/api/stac/v1", + collection = "sentinel-2-l2a", + cube = TRUE, + # Planetary Computer asset names (NOT the Element84/AWS red/nir/scl names) + roles = list( + red = "B04", + nir = "B08", + swir16 = "B11", + mask = "SCL" + ), + # Scene Classification Layer classes to mask out: + # 3 cloud shadow, 8 cloud medium, 9 cloud high, 10 thin cirrus, 11 snow/ice + # (snow is masked by default: it is never a vegetation signal and otherwise + # dominates the winter trajectory at high latitudes) + mask_values = c(3L, 8L, 9L, 10L, 11L), + # DN -> reflectance. The +1000 DN offset only applies from processing + # baseline 04.00 (2022-01-25 on); earlier scenes have no offset. dft_stac_cube + # splits items at the boundary and applies `offset` after it, `offset_before` + # before it, so a multi-year series has no artificial step at 2022. + scale = 1e-4, + offset = -0.1, + offset_boundary = "2022-01-25", + offset_before = 0, + available_datetime = "2017-01-01/2024-12-31" ) + # "landsat-c2-l2" drops in with the same cube shape, no API change: + # list( + # # Cost: free — Planetary Computer (Azure Open Data), no egress charge. + # stac_url = "https://planetarycomputer.microsoft.com/api/stac/v1", + # collection = "landsat-c2-l2", + # cube = TRUE, + # roles = list(red = "red", nir = "nir08", + # swir16 = "swir16", mask = "qa_pixel"), + # mask_values = , + # scale = 0.0000275, offset = -0.2, # non-zero offset breaks ratio + # available_datetime = "1984-01-01/2024-12-31" # indices on raw DNs + # ) ) } diff --git a/R/dft_stac_cube.R b/R/dft_stac_cube.R new file mode 100644 index 0000000..711ebc0 --- /dev/null +++ b/R/dft_stac_cube.R @@ -0,0 +1,302 @@ +#' Fetch a masked spectral-index cube from a STAC catalog +#' +#' Sibling of [dft_stac_fetch()] for continuous change detection. Where +#' [dft_stac_fetch()] materializes one categorical raster per year, this builds +#' a sub-annual reflectance cube, masks clouds, computes a spectral index over +#' band roles, and returns the index time series as a `SpatRaster` (one layer +#' per time step) — the input to [dft_rast_break()] for per-pixel trajectory +#' breakpoint detection. +#' +#' The index stack is materialized once to a GeoTIFF under [dft_cache_path()] +#' as `/cube_.tif`, keyed by a hash of the AOI geometry and every +#' cube-affecting parameter. Because it is invariant to [dft_rast_break()]'s +#' parameters, caching it here makes bfast parameter sweeps cheap — they re-read +#' the local raster instead of re-streaming COGs. +#' +#' Three STAC-query specifics distinguish cube mode from [dft_stac_fetch()]: +#' pagination via [rstac::items_fetch()] is mandatory (a monthly multi-year query +#' returns hundreds of items; a single page silently truncates); the query uses +#' `intersects` with the AOI geometry, not a bounding box (floodplain polygons +#' are highly non-rectangular); and a scene-level `eo:cloud_cover` pre-filter +#' shrinks the collection before any pixel is read, complementing per-pixel mask +#' filtering. +#' +#' @param aoi An `sf` polygon defining the area of interest. +#' @param source Character. A cube source name for [dft_stac_config()] (default +#' `"sentinel-2-l2a"`). Must be a source with `cube = TRUE`. +#' @param index Character. Spectral index from [dft_index_table()] (default +#' `"kndvi"`). Determines which band roles (and thus assets) are fetched. +#' @param datetime Character. ISO 8601 interval `"start/end"`. When `NULL`, uses +#' `available_datetime` from [dft_stac_config()]. +#' @param res Numeric. Output pixel size in CRS units (default 10). +#' @param crs Character. Target CRS as an EPSG string. When `NULL`, auto-detected +#' from the AOI centroid's UTM zone. +#' @param dt Character. ISO 8601 duration for the temporal aggregation window +#' (default `"P1M"`, monthly). The cadence [dft_rast_break()]'s `frequency` +#' must agree with. +#' @param aggregation Character. Temporal aggregation for multiple scenes in one +#' `dt` window (default `"median"`). +#' @param resampling Character. Spatial resampling (default `"bilinear"`). +#' @param cloud_cover_max Numeric. Scene-level `eo:cloud_cover` maximum percent +#' for the STAC pre-filter (default 60). +#' @param months Integer vector of calendar months (1-12) to keep, or `NULL` +#' (default) for all. Restricting to the growing season (e.g. `6:9`) both +#' sharpens the vegetation signal — snow and low-sun winter scenes carry no +#' vegetation information — and cuts the number of scenes streamed. Months with +#' no retained scenes become `NA` in the monthly cube, so the per-pixel series +#' stays regular at `frequency = 12` for [dft_rast_break()]. Prefer a longer +#' `datetime` window when using this, so enough growing-season history remains +#' to fit a stable BFAST baseline. +#' @param mask_values Integer vector of mask-band classes to exclude. When +#' `NULL`, uses `mask_values` from [dft_stac_config()] (e.g. Sentinel-2 SCL +#' cloud / shadow / cirrus classes). +#' @param cache_dir Character. Cache directory. When `NULL`, uses +#' [dft_cache_path()]. +#' @param force Logical. Re-fetch even if cached, overwriting the cached raster +#' (default `FALSE`). +#' @param sign_fn A signing function for STAC assets. Default is +#' [rstac::sign_planetary_computer()]. +#' +#' @return A [terra::SpatRaster] index stack — one layer per time step, with a +#' time value per layer — cached as a GeoTIFF. The stack spans the AOI +#' **bounding box** (cloud-masked but not clipped to the AOI polygon); clip the +#' reduced raster from [dft_rast_break()] with `terra::mask()` if a tight AOI is +#' needed. For sources with a reflectance-offset baseline boundary (Sentinel-2), +#' items are split at the boundary and offset-corrected per side, so a series +#' crossing it carries no artificial index step. +#' +#' @seealso [dft_rast_break()] (the reducer that consumes this cube), +#' [dft_index_expr()] (the index applied), [dft_stac_fetch()] (categorical +#' sibling). +#' +#' @examples +#' \dontrun{ +#' # Monthly kNDVI cube for a floodplain reach (requires network + gdalcubes) +#' aoi <- sf::st_read(system.file("extdata", "example_aoi.gpkg", package = "drift")) +#' cube <- dft_stac_cube( +#' aoi, +#' source = "sentinel-2-l2a", +#' index = "kndvi", +#' datetime = "2019-01-01/2023-12-31", +#' dt = "P1M" +#' ) +#' breaks <- dft_rast_break(cube, start = c(2022, 1)) +#' } +#' +#' @export +dft_stac_cube <- function(aoi, + source = "sentinel-2-l2a", + index = "kndvi", + datetime = NULL, + res = 10, + crs = NULL, + dt = "P1M", + aggregation = "median", + resampling = "bilinear", + cloud_cover_max = 60, + months = NULL, + mask_values = NULL, + cache_dir = NULL, + force = FALSE, + sign_fn = rstac::sign_planetary_computer()) { + rlang::check_installed("gdalcubes", reason = "to fetch STAC cubes") + + # GDAL cloud-read tuning for /vsicurl COG streaming (biggest win: + # DISABLE_READDIR_ON_OPEN avoids a remote directory listing on every open). + # Restored on exit so we don't mutate the caller's session. + gdal_cfg <- c( + GDAL_DISABLE_READDIR_ON_OPEN = "EMPTY_DIR", + GDAL_HTTP_MULTIPLEX = "YES", + GDAL_HTTP_VERSION = "2", + VSI_CACHE = "TRUE", + CPL_VSIL_CURL_ALLOWED_EXTENSIONS = ".tif" + ) + old_cfg <- Sys.getenv(names(gdal_cfg), unset = NA) + do.call(Sys.setenv, as.list(gdal_cfg)) + on.exit({ + set_again <- old_cfg[!is.na(old_cfg)] + if (length(set_again)) do.call(Sys.setenv, as.list(set_again)) + unset <- names(old_cfg)[is.na(old_cfg)] + if (length(unset)) Sys.unsetenv(unset) + }, add = TRUE) + + cfg <- dft_stac_config(source) + if (!isTRUE(cfg$cube)) { + cli::cli_abort(c( + "Source {.val {source}} is not a cube source.", + "i" = "Use {.fn dft_stac_fetch} for categorical rasters." + )) + } + datetime <- datetime %||% cfg$available_datetime + mask_values <- mask_values %||% cfg$mask_values + scale <- cfg$scale %||% 1 + offset <- cfg$offset %||% 0 + # sources whose reflectance offset changes at a processing-baseline boundary + # (e.g. Sentinel-2 +1000 DN from 2022-01-25) carry the boundary date and the + # pre-boundary offset; the fetch splits items at the boundary and corrects + # each side, so a series crossing it has no artificial index step. + offset_boundary <- cfg$offset_boundary + offset_before <- cfg$offset_before %||% 0 + + # Ensure aoi is sf + if (inherits(aoi, "SpatVector")) aoi <- sf::st_as_sf(aoi) + stopifnot(inherits(aoi, c("sf", "sfc"))) + + target_crs <- if (is.null(crs)) auto_utm_epsg(aoi) else crs + aoi_wgs84 <- sf::st_transform(aoi, 4326) + aoi_target <- sf::st_transform(aoi, as.integer(gsub("EPSG:", "", target_crs))) + bbox_target <- sf::st_bbox(aoi_target) + + # Assets: the index's required roles + the mask role (index_roles errors on + # an unknown index before any network call) + roles_needed <- index_roles(index) + band_assets <- unlist(cfg$roles[roles_needed], use.names = FALSE) + mask_asset <- cfg$roles$mask + + # datetime interval -> cube_view time bounds + dr <- strsplit(datetime, "/", fixed = TRUE)[[1]] + if (length(dr) != 2) { + cli::cli_abort("`datetime` must be an ISO 8601 interval {.val start/end}.") + } + t0 <- dr[1] + t1 <- dr[2] + + # Cache + cache_base <- dft_cache_path(cache_dir) + cache_source_dir <- file.path(cache_base, source) + dir.create(cache_source_dir, recursive = TRUE, showWarnings = FALSE) + cache_key <- stac_cube_cache_key( + aoi_target, res, target_crs, dt, aggregation, resampling, + cfg$stac_url, cfg$collection, band_assets, datetime, index, + cloud_cover_max, mask_values, scale, offset, months, offset_before + ) + cache_file <- file.path(cache_source_dir, paste0("cube_", cache_key, ".tif")) + + # monthly layer times, derived from the datetime window start + month_times <- function(n) { + seq(as.Date(paste0(substr(dr[1], 1, 7), "-01")), by = "month", length.out = n) + } + + if (!force && file.exists(cache_file)) { + message(" cube: cached") + r <- terra::rast(cache_file) + terra::time(r) <- month_times(terra::nlyr(r)) + return(r) + } + + # STAC query: intersects (not bbox) + scene cloud pre-filter + pagination + message("Querying STAC: ", cfg$collection, " (", datetime, ")...") + items <- rstac::stac(cfg$stac_url) |> + rstac::stac_search( + collections = cfg$collection, + # union so a multi-feature AOI queries its whole footprint, matching the + # cube extent and filter_geom clip (a single first-feature geometry would + # leave silent NoData holes over the other features) + intersects = sf::st_geometry(sf::st_union(aoi_wgs84))[[1]], + datetime = datetime, + limit = 500 + ) |> + rstac::ext_filter(`eo:cloud_cover` <= {{cloud_cover_max}}) |> + rstac::post_request() |> + rstac::items_fetch() |> + rstac::items_sign(sign_fn = sign_fn) + + # Restrict to growing-season (or any) calendar months. Fetching fewer, better + # months both sharpens the vegetation signal (drops snow/low-sun winter noise) + # and cuts the number of scenes streamed. Months with no scenes become NA in + # the monthly cube, so the ts() stays regular at frequency 12. + if (!is.null(months)) { + item_dt <- vapply(items$features, function(f) f$properties$datetime %||% NA_character_, "") + item_mo <- as.integer(format(as.Date(substr(item_dt, 1, 10)), "%m")) + items$features <- items$features[!is.na(item_mo) & item_mo %in% months] + } + + n_items <- length(items$features) + message(" ", n_items, " items returned") + if (n_items == 0) stop("No STAC items found for ", cfg$collection) + + v <- gdalcubes::cube_view( + srs = target_crs, + extent = list( + left = bbox_target[["xmin"]], right = bbox_target[["xmax"]], + bottom = bbox_target[["ymin"]], top = bbox_target[["ymax"]], + t0 = t0, t1 = t1 + ), + dx = res, dy = res, dt = dt, + aggregation = aggregation, resampling = resampling + ) + + # Build the index cube for one item subset with one offset, materialize it, and + # read it back as a terra stack. The cube spans the AOI bounding box: + # gdalcubes::filter_geom() to clip to the polygon yields an all-NA cube (and can + # crash the compute worker) on the pinned build, so we mask clouds here and + # leave polygon clipping to the caller, as the sibling dft_stac_fetch() does. + build_index_stack <- function(features, offset_use) { + img_col <- gdalcubes::stac_image_collection( + features, asset_names = c(band_assets, mask_asset) + ) + cube <- gdalcubes::raster_cube( + img_col, v, mask = gdalcubes::image_mask(mask_asset, values = mask_values) + ) + idx <- dft_index_expr(cube, index = index, source = source, + roles = cfg$roles, scale = scale, offset = offset_use) + tmp <- tempfile(fileext = ".nc") + gdalcubes::write_ncdf(idx, tmp, overwrite = TRUE) + terra::rast(tmp) + } + + # Baseline-conditional offset: split items at the boundary and correct each + # side with its own offset, then coalesce onto the shared monthly grid. Both + # subcubes are built over the full view so their layers align for terra::cover. + is_pre <- rep(FALSE, length(items$features)) + if (!is.null(offset_boundary)) { + item_date <- as.Date(substr( + vapply(items$features, function(f) f$properties$datetime %||% NA_character_, ""), + 1, 10 + )) + is_pre <- !is.na(item_date) & item_date < as.Date(offset_boundary) + } + + if (any(is_pre) && !all(is_pre)) { + message(" offset split at ", offset_boundary, ": ", + sum(is_pre), " pre / ", sum(!is_pre), " post") + stk <- terra::cover( + build_index_stack(items$features[is_pre], offset_before), + build_index_stack(items$features[!is_pre], offset) + ) + } else { + stk <- build_index_stack(items$features, if (all(is_pre)) offset_before else offset) + } + + terra::time(stk) <- month_times(terra::nlyr(stk)) + names(stk) <- rep(index, terra::nlyr(stk)) + terra::writeRaster(stk, cache_file, overwrite = TRUE) + stk +} + + +#' Cache key for one STAC index-cube parameter set +#' +#' Cube-mode analogue of `stac_cache_key()` (kept separate so the fetch key +#' stays byte-for-byte stable). Hashes the AOI geometry as WKB plus every +#' parameter that changes the written index cube. `res` is coerced to double so +#' `10L` and `10` key alike; `mask_values` is sorted so order does not matter. +#' `scale`/`offset` are included because they change pixel values. +#' @noRd +stac_cube_cache_key <- function(aoi_target, res, target_crs, dt, aggregation, + resampling, stac_url, collection, band_assets, + datetime, index, cloud_cover_max, mask_values, + scale, offset, months = NULL, + offset_before = 0) { + geom_wkb <- sf::st_as_binary(sf::st_geometry(aoi_target), endian = "little") + substr( + rlang::hash(list( + geom_wkb, as.numeric(res), target_crs, dt, aggregation, resampling, + stac_url, collection, band_assets, datetime, index, + as.numeric(cloud_cover_max), sort(as.numeric(mask_values)), + as.numeric(scale), as.numeric(offset), sort(as.numeric(months)), + as.numeric(offset_before) + )), + 1, 12 + ) +} diff --git a/R/dft_stac_fetch.R b/R/dft_stac_fetch.R index e72e202..392f116 100644 --- a/R/dft_stac_fetch.R +++ b/R/dft_stac_fetch.R @@ -61,6 +61,12 @@ dft_stac_fetch <- function(aoi, # Resolve config if (is.null(stac_url) || is.null(collection) || is.null(asset)) { cfg <- dft_stac_config(source) + if (isTRUE(cfg$cube)) { + cli::cli_abort(c( + "Source {.val {source}} is a cube source, not a categorical raster.", + "i" = "Use {.fn dft_stac_cube} for continuous index-trajectory sources." + )) + } stac_url <- stac_url %||% cfg$stac_url collection <- collection %||% cfg$collection asset <- asset %||% cfg$asset diff --git a/R/drift-package.R b/R/drift-package.R index 0cc9deb..f7e8a6a 100644 --- a/R/drift-package.R +++ b/R/drift-package.R @@ -1,2 +1,6 @@ #' @keywords internal "_PACKAGE" + +# NSE symbols: `eo:cloud_cover` in rstac::ext_filter() CQL2 (dft_stac_cube); +# `.data` is the rlang pronoun used in dplyr pipelines (dft_rast_classify) +utils::globalVariables(c("eo:cloud_cover", ".data")) diff --git a/data-raw/vignette_data_break.R b/data-raw/vignette_data_break.R new file mode 100644 index 0000000..db1065f --- /dev/null +++ b/data-raw/vignette_data_break.R @@ -0,0 +1,53 @@ +# Generate the committed artifact for vignettes/trajectory-break-detection.Rmd. +# +# NOT part of the package build. Requires network (Microsoft Planetary +# Computer), `gdalcubes`, and `bfast`. Follows drift's data-raw pattern: do the +# network + bfast work here, save a small artifact to inst/testdata/, and let +# the vignette load that artifact so every chunk runs live under pkgdown with no +# network or bfast dependency. +# +# Window (2017-2023) and reach match the land-cover vignette. The artifact +# carries the two continuous reductions (abrupt break, gradual trend), the LULC +# Trees->* transition polygons (for the interactive map + overlap check), and one +# real harvested patch's trajectory. The Sentinel-2 fetch is ~10-15 min (once; +# cached); the reductions are seconds. +# Rscript data-raw/vignette_data_break.R + +library(drift) + +aoi <- sf::st_read( + system.file("extdata", "example_aoi.gpkg", package = "drift"), + quiet = TRUE +) +cache <- file.path("data-raw", ".break_cache") + +cube <- dft_stac_cube( + aoi, source = "sentinel-2-l2a", index = "kndvi", + datetime = "2017-01-01/2023-12-31", dt = "P1M", months = 6:9, + cloud_cover_max = 60, cache_dir = cache +) + +breaks <- dft_rast_break(cube, start = c(2022, 1), order = 1) # abrupt: the cuts +trend <- dft_rast_trend(cube) # gradual: which way + +# ---- LULC Trees -> * transition polygons (the categorical "harvest" signal) +l17 <- terra::rast(system.file("extdata", "example_2017.tif", package = "drift")) +l23 <- terra::rast(system.file("extdata", "example_2023.tif", package = "drift")) +classified <- dft_rast_classify(list("2017" = l17, "2023" = l23), source = "io-lulc") +tt <- dft_rast_transition(classified, from = "2017", to = "2023") +tv <- dft_transition_vectors(tt$raster) +# just the actual loss transition (exclude stable "Trees -> Trees") +tree_trans <- tv[tv$transition == "Trees -> Rangeland", ] + +aoi_v <- terra::vect(sf::st_transform(aoi, terra::crs(breaks))) + +dir.create("inst/testdata", recursive = TRUE, showWarnings = FALSE) +saveRDS( + list(trend_r = terra::wrap(terra::mask(trend, aoi_v)), + break_r = terra::wrap(terra::mask(breaks[["break_mag"]], aoi_v)), + tree_trans = tree_trans, aoi = aoi), + "inst/testdata/neexdzii_break.rds" +) + +message("Wrote inst/testdata/neexdzii_break.rds") +message(" Trees->* transition polygons: ", nrow(tree_trans)) diff --git a/inst/indices/indices.csv b/inst/indices/indices.csv new file mode 100644 index 0000000..26f0c21 --- /dev/null +++ b/inst/indices/indices.csv @@ -0,0 +1,4 @@ +index,formula,roles,description +ndvi,"(nir - red) / (nir + red)","nir,red",Normalized Difference Vegetation Index +kndvi,"tanh(pow((nir - red) / (nir + red), 2))","nir,red",Kernel NDVI (tanh of NDVI squared) +ndmi,"(nir - swir16) / (nir + swir16)","nir,swir16",Normalized Difference Moisture Index diff --git a/inst/testdata/neexdzii_break.rds b/inst/testdata/neexdzii_break.rds new file mode 100644 index 0000000..8f4b149 Binary files /dev/null and b/inst/testdata/neexdzii_break.rds differ diff --git a/man/dft_index_expr.Rd b/man/dft_index_expr.Rd new file mode 100644 index 0000000..5efbb47 --- /dev/null +++ b/man/dft_index_expr.Rd @@ -0,0 +1,68 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/dft_index_expr.R +\name{dft_index_expr} +\alias{dft_index_expr} +\title{Apply a spectral index to a data cube} +\usage{ +dft_index_expr( + cube, + index = "kndvi", + source = "sentinel-2-l2a", + roles = NULL, + scale = NULL, + offset = NULL +) +} +\arguments{ +\item{cube}{A \code{gdalcubes} data cube (e.g. the lazy cube built inside +\code{\link[=dft_stac_cube]{dft_stac_cube()}}) whose bands are the source's assets.} + +\item{index}{Character. An index name present in \code{\link[=dft_index_table]{dft_index_table()}} +(default \code{"kndvi"}).} + +\item{source}{Character. Source name passed to \code{\link[=dft_stac_config]{dft_stac_config()}} to resolve +the role→asset map and reflectance scale/offset (default +\code{"sentinel-2-l2a"}).} + +\item{roles}{Named list mapping roles to asset names. When \code{NULL}, taken from +\code{dft_stac_config(source)$roles}.} + +\item{scale, offset}{Numeric reflectance affine transform. When \code{NULL}, taken +from the source config (falling back to \code{1} / \code{0}).} +} +\value{ +A single-band \code{gdalcubes} cube with the band named \code{index}. +} +\description{ +Resolve a named spectral index (e.g. \code{"kndvi"}) into a per-pixel arithmetic +expression over a source's band roles and apply it to a \code{gdalcubes} data +cube, returning a single-band cube named after the index. +} +\details{ +Index formulas are stored in a shipped registry (\code{\link[=dft_index_table]{dft_index_table()}}) written +over band \strong{roles} (\code{red}, \code{nir}, \code{swir16}), not literal asset names. The +roles are resolved to per-source asset names via \code{\link[=dft_stac_config]{dft_stac_config()}}, so the +same \code{"kndvi"} works against Sentinel-2 (\code{B04}/\code{B08}) or any future +reflectance source without changing the formula. + +Reflectance \code{scale}/\code{offset} are folded \strong{into} the expression as a per-band +affine transform \code{(asset * scale + offset)}. This matters for ratio indices: +a non-zero offset does not cancel in \code{(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). +} +\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 +} + +} +\seealso{ +\code{\link[=dft_index_table]{dft_index_table()}} for the registry, \code{\link[=dft_stac_cube]{dft_stac_cube()}} for the +caller that builds the input cube. +} diff --git a/man/dft_index_table.Rd b/man/dft_index_table.Rd new file mode 100644 index 0000000..4cb983c --- /dev/null +++ b/man/dft_index_table.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/dft_index_expr.R +\name{dft_index_table} +\alias{dft_index_table} +\title{Load the shipped spectral-index registry} +\usage{ +dft_index_table() +} +\value{ +A tibble with columns \code{index}, \code{formula}, \code{roles} (comma-separated +role names), and \code{description}. +} +\description{ +Reads the CSV index registry bundled with the package. Each row defines one +index as a \code{gdalcubes}/tinyexpr formula written over band roles. +} +\examples{ +dft_index_table() + +} diff --git a/man/dft_rast_break.Rd b/man/dft_rast_break.Rd new file mode 100644 index 0000000..6e4b152 --- /dev/null +++ b/man/dft_rast_break.Rd @@ -0,0 +1,83 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/dft_rast_break.R +\name{dft_rast_break} +\alias{dft_rast_break} +\title{Detect per-pixel index-trajectory breakpoints} +\usage{ +dft_rast_break( + cube, + history = "all", + start = c(2022, 1), + frequency = NULL, + order = 3, + level = 0.01, + min_obs = 6, + cores = NULL +) +} +\arguments{ +\item{cube}{A monthly index \code{SpatRaster} (the return value of +\code{\link[=dft_stac_cube]{dft_stac_cube()}}): one layer per time step, with a time value per layer.} + +\item{history}{Character. \code{bfastmonitor} history-selection method: \code{"all"} +(default), \code{"ROC"}, or \code{"BP"}.} + +\item{start}{Numeric \code{c(year, period)}. Start of the monitoring period, in the +stack's temporal frequency (e.g. \code{c(2022, 1)} = Jan 2022 for a monthly +stack). Everything before it is the stable history.} + +\item{frequency}{Numeric or \code{NULL}. Seasonal frequency of the time series +(12 for monthly, 1 for annual). When \code{NULL}, derived from the layer time +spacing; when supplied, it must agree with that spacing or the call errors.} + +\item{order}{Integer. Harmonic order of the season-trend model passed to +\code{\link[bfast:bfastmonitor]{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 +\code{\link[=dft_stac_cube]{dft_stac_cube()}} \code{months}), where a high order overfits sparse seasonal +coverage.} + +\item{level}{Numeric. Significance level passed to \code{\link[bfast:bfastmonitor]{bfast::bfastmonitor()}} +(default 0.01).} + +\item{min_obs}{Integer. Minimum non-\code{NA} observations required to attempt a +fit; pixels with fewer return \code{NA} (default 6).} + +\item{cores}{Integer or \code{NULL}. Forked workers for the per-pixel reduction. +When \code{NULL}, uses one fewer than the detected cores.} +} +\value{ +A two-band \link[terra:SpatRaster-class]{terra::SpatRaster} with layers \code{break_date} (decimal year +or \code{NA}) and \code{break_mag} (signed index change; negative = index drop). +} +\description{ +Reduce a monthly index stack (from \code{\link[=dft_stac_cube]{dft_stac_cube()}}) over time with +\code{\link[bfast:bfastmonitor]{bfast::bfastmonitor()}}, returning a two-band \code{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: \emph{when} did the pixel's spectral history break, and by how much? +} +\details{ +\code{bfastmonitor} fits a season-trend model to a stable \emph{history} period, then +watches the \emph{monitoring} period (from \code{start} onward) for a structural break. +The returned \code{break_mag} is the median monitoring-period residual: \strong{negative +means the index dropped} (e.g. vegetation loss / channel scour), positive +means it rose (establishment). \code{break_date} is a decimal year (e.g. +\code{2022.42}) or \code{NA} where no break was detected. + +Pixels are reduced in parallel with \code{parallel::mclapply()} (forked workers, so +the per-pixel logic and its parameters are inherited directly). Pixels with +fewer than \code{min_obs} valid observations short-circuit to \code{NA}. +} +\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 +} + +} +\seealso{ +\code{\link[=dft_stac_cube]{dft_stac_cube()}} (builds the input stack), \code{\link[=dft_index_expr]{dft_index_expr()}}. +} diff --git a/man/dft_rast_trend.Rd b/man/dft_rast_trend.Rd new file mode 100644 index 0000000..1cb5c8d --- /dev/null +++ b/man/dft_rast_trend.Rd @@ -0,0 +1,48 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/dft_rast_trend.R +\name{dft_rast_trend} +\alias{dft_rast_trend} +\title{Per-pixel index trend over time} +\usage{ +dft_rast_trend(cube, min_obs = 6, cores = NULL) +} +\arguments{ +\item{cube}{A monthly index \code{SpatRaster} (the return value of +\code{\link[=dft_stac_cube]{dft_stac_cube()}}): one layer per time step, with a time value per layer.} + +\item{min_obs}{Integer. Minimum non-\code{NA} observations required; pixels with +fewer return \code{NA} (default 6).} + +\item{cores}{Integer or \code{NULL}. Forked workers for the per-pixel reduction. +When \code{NULL}, uses one fewer than the detected cores.} +} +\value{ +A two-band \link[terra:SpatRaster-class]{terra::SpatRaster} with layers \code{trend} (index change per +year — negative = declining greenness, positive = recovering) and \code{trend_p} +(Mann-Kendall two-sided p-value; small = a significant monotonic trend). +} +\description{ +Reduce an index stack (from \code{\link[=dft_stac_cube]{dft_stac_cube()}}) to a per-pixel \strong{trend} — the +rate at which the index rises or falls across the whole record. Where +\code{\link[=dft_rast_break]{dft_rast_break()}} finds an abrupt structural break, this measures the slow, +monotonic direction: gradual degradation (a declining stand) or recovery (a +restoration planting greening up) that annual categorical labels cannot show. +} +\details{ +The slope is a \strong{Theil-Sen} estimate — the median of all pairwise slopes — +which is resistant to a single anomalous season (e.g. a smoke/drought year), +unlike an ordinary least-squares fit. Significance is the non-parametric +\strong{Mann-Kendall} test for a monotonic trend. +} +\examples{ +\dontrun{ +aoi <- sf::st_read(system.file("extdata", "example_aoi.gpkg", package = "drift")) +cube <- dft_stac_cube(aoi, index = "kndvi", datetime = "2017-01-01/2023-12-31") +trend <- dft_rast_trend(cube) +terra::plot(trend[["trend"]]) # negative (red) = declining, positive = recovering +} + +} +\seealso{ +\code{\link[=dft_rast_break]{dft_rast_break()}} (abrupt breaks), \code{\link[=dft_stac_cube]{dft_stac_cube()}}. +} diff --git a/man/dft_stac_config.Rd b/man/dft_stac_config.Rd index 1790af8..1e5d112 100644 --- a/man/dft_stac_config.Rd +++ b/man/dft_stac_config.Rd @@ -2,30 +2,43 @@ % Please edit documentation in R/dft_stac_config.R \name{dft_stac_config} \alias{dft_stac_config} -\title{Get STAC configuration for a known land cover source} +\title{Get STAC configuration for a known source} \usage{ -dft_stac_config(source = c("io-lulc", "esa-worldcover")) +dft_stac_config(source = c("io-lulc", "esa-worldcover", "sentinel-2-l2a")) } \arguments{ -\item{source}{Character. One of \code{"io-lulc"} (Esri IO LULC annual v02) -or \code{"esa-worldcover"} (ESA WorldCover).} +\item{source}{Character. One of \code{"io-lulc"} (Esri IO LULC annual v02), +\code{"esa-worldcover"} (ESA WorldCover), or \code{"sentinel-2-l2a"} (Sentinel-2 +L2A surface reflectance, a cube source).} } \value{ -A list with elements: -\describe{ -\item{stac_url}{STAC API endpoint} -\item{collection}{Collection ID} -\item{asset}{Asset name to download} -\item{available_years}{Integer vector of available years} -} +A list. Categorical sources have elements \code{stac_url}, \code{collection}, +\code{asset}, \code{available_years}. Cube sources have \code{stac_url}, \code{collection}, +\code{cube = TRUE}, \code{roles} (a named list mapping \code{red}/\code{nir}/\code{swir16}/\code{mask} +to asset names), \code{mask_values} (integer mask classes to exclude), +\code{scale}/\code{offset} (DN → reflectance affine transform), and +\code{available_datetime} (an ISO 8601 interval string). The \code{cube} field is +absent (not \code{FALSE}) for categorical sources; test with \code{isTRUE(cfg$cube)}. } \description{ Returns connection details for pre-configured STAC collections. -Used as a convenience wrapper around \code{\link[=dft_stac_fetch]{dft_stac_fetch()}} so users -don't need to remember STAC URLs and collection IDs. +Used as a convenience wrapper around \code{\link[=dft_stac_fetch]{dft_stac_fetch()}} (categorical +sources) and \code{\link[=dft_stac_cube]{dft_stac_cube()}} (continuous index-trajectory sources) so +users don't need to remember STAC URLs, collection IDs, and band names. +} +\details{ +Sources are of two kinds. \strong{Categorical} sources (\code{"io-lulc"}, +\code{"esa-worldcover"}) host single-band classified rasters and carry a flat +\code{asset} name for \code{\link[=dft_stac_fetch]{dft_stac_fetch()}}. \strong{Cube} sources (\code{"sentinel-2-l2a"}) +host multi-band reflectance imagery and instead carry a role-based band map +(\code{red}/\code{nir}/\code{swir16}/\code{mask}), mask values, and reflectance scale/offset for +\code{\link[=dft_stac_cube]{dft_stac_cube()}}; they are marked with \code{cube = TRUE}. The role-based schema +means a new reflectance source (e.g. Landsat C2 L2) drops in with no API +change — only the role→asset map and scale/offset differ. } \examples{ dft_stac_config("io-lulc") dft_stac_config("esa-worldcover") +dft_stac_config("sentinel-2-l2a") } diff --git a/man/dft_stac_cube.Rd b/man/dft_stac_cube.Rd new file mode 100644 index 0000000..df566fe --- /dev/null +++ b/man/dft_stac_cube.Rd @@ -0,0 +1,127 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/dft_stac_cube.R +\name{dft_stac_cube} +\alias{dft_stac_cube} +\title{Fetch a masked spectral-index cube from a STAC catalog} +\usage{ +dft_stac_cube( + aoi, + source = "sentinel-2-l2a", + index = "kndvi", + datetime = NULL, + res = 10, + crs = NULL, + dt = "P1M", + aggregation = "median", + resampling = "bilinear", + cloud_cover_max = 60, + months = NULL, + mask_values = NULL, + cache_dir = NULL, + force = FALSE, + sign_fn = rstac::sign_planetary_computer() +) +} +\arguments{ +\item{aoi}{An \code{sf} polygon defining the area of interest.} + +\item{source}{Character. A cube source name for \code{\link[=dft_stac_config]{dft_stac_config()}} (default +\code{"sentinel-2-l2a"}). Must be a source with \code{cube = TRUE}.} + +\item{index}{Character. Spectral index from \code{\link[=dft_index_table]{dft_index_table()}} (default +\code{"kndvi"}). Determines which band roles (and thus assets) are fetched.} + +\item{datetime}{Character. ISO 8601 interval \code{"start/end"}. When \code{NULL}, uses +\code{available_datetime} from \code{\link[=dft_stac_config]{dft_stac_config()}}.} + +\item{res}{Numeric. Output pixel size in CRS units (default 10).} + +\item{crs}{Character. Target CRS as an EPSG string. When \code{NULL}, auto-detected +from the AOI centroid's UTM zone.} + +\item{dt}{Character. ISO 8601 duration for the temporal aggregation window +(default \code{"P1M"}, monthly). The cadence \code{\link[=dft_rast_break]{dft_rast_break()}}'s \code{frequency} +must agree with.} + +\item{aggregation}{Character. Temporal aggregation for multiple scenes in one +\code{dt} window (default \code{"median"}).} + +\item{resampling}{Character. Spatial resampling (default \code{"bilinear"}).} + +\item{cloud_cover_max}{Numeric. Scene-level \code{eo:cloud_cover} maximum percent +for the STAC pre-filter (default 60).} + +\item{months}{Integer vector of calendar months (1-12) to keep, or \code{NULL} +(default) for all. Restricting to the growing season (e.g. \code{6:9}) both +sharpens the vegetation signal — snow and low-sun winter scenes carry no +vegetation information — and cuts the number of scenes streamed. Months with +no retained scenes become \code{NA} in the monthly cube, so the per-pixel series +stays regular at \code{frequency = 12} for \code{\link[=dft_rast_break]{dft_rast_break()}}. Prefer a longer +\code{datetime} window when using this, so enough growing-season history remains +to fit a stable BFAST baseline.} + +\item{mask_values}{Integer vector of mask-band classes to exclude. When +\code{NULL}, uses \code{mask_values} from \code{\link[=dft_stac_config]{dft_stac_config()}} (e.g. Sentinel-2 SCL +cloud / shadow / cirrus classes).} + +\item{cache_dir}{Character. Cache directory. When \code{NULL}, uses +\code{\link[=dft_cache_path]{dft_cache_path()}}.} + +\item{force}{Logical. Re-fetch even if cached, overwriting the cached raster +(default \code{FALSE}).} + +\item{sign_fn}{A signing function for STAC assets. Default is +\code{\link[rstac:items_sign_planetary_computer]{rstac::sign_planetary_computer()}}.} +} +\value{ +A \link[terra:SpatRaster-class]{terra::SpatRaster} index stack — one layer per time step, with a +time value per layer — cached as a GeoTIFF. The stack spans the AOI +\strong{bounding box} (cloud-masked but not clipped to the AOI polygon); clip the +reduced raster from \code{\link[=dft_rast_break]{dft_rast_break()}} with \code{terra::mask()} if a tight AOI is +needed. For sources with a reflectance-offset baseline boundary (Sentinel-2), +items are split at the boundary and offset-corrected per side, so a series +crossing it carries no artificial index step. +} +\description{ +Sibling of \code{\link[=dft_stac_fetch]{dft_stac_fetch()}} for continuous change detection. Where +\code{\link[=dft_stac_fetch]{dft_stac_fetch()}} materializes one categorical raster per year, this builds +a sub-annual reflectance cube, masks clouds, computes a spectral index over +band roles, and returns the index time series as a \code{SpatRaster} (one layer +per time step) — the input to \code{\link[=dft_rast_break]{dft_rast_break()}} for per-pixel trajectory +breakpoint detection. +} +\details{ +The index stack is materialized once to a GeoTIFF under \code{\link[=dft_cache_path]{dft_cache_path()}} +as \verb{/cube_.tif}, keyed by a hash of the AOI geometry and every +cube-affecting parameter. Because it is invariant to \code{\link[=dft_rast_break]{dft_rast_break()}}'s +parameters, caching it here makes bfast parameter sweeps cheap — they re-read +the local raster instead of re-streaming COGs. + +Three STAC-query specifics distinguish cube mode from \code{\link[=dft_stac_fetch]{dft_stac_fetch()}}: +pagination via \code{\link[rstac:items_functions]{rstac::items_fetch()}} is mandatory (a monthly multi-year query +returns hundreds of items; a single page silently truncates); the query uses +\code{intersects} with the AOI geometry, not a bounding box (floodplain polygons +are highly non-rectangular); and a scene-level \code{eo:cloud_cover} pre-filter +shrinks the collection before any pixel is read, complementing per-pixel mask +filtering. +} +\examples{ +\dontrun{ +# Monthly kNDVI cube for a floodplain reach (requires network + gdalcubes) +aoi <- sf::st_read(system.file("extdata", "example_aoi.gpkg", package = "drift")) +cube <- dft_stac_cube( + aoi, + source = "sentinel-2-l2a", + index = "kndvi", + datetime = "2019-01-01/2023-12-31", + dt = "P1M" +) +breaks <- dft_rast_break(cube, start = c(2022, 1)) +} + +} +\seealso{ +\code{\link[=dft_rast_break]{dft_rast_break()}} (the reducer that consumes this cube), +\code{\link[=dft_index_expr]{dft_index_expr()}} (the index applied), \code{\link[=dft_stac_fetch]{dft_stac_fetch()}} (categorical +sibling). +} diff --git a/planning/active/findings.md b/planning/active/findings.md new file mode 100644 index 0000000..14cc650 --- /dev/null +++ b/planning/active/findings.md @@ -0,0 +1,82 @@ +# Findings — Trajectory vignette redo (#30, continued) + +**Continues** `planning/archive/2026-07-issue-30-index-trajectory/findings.md` — read that for the +full empirical journey. Condensed carry-forward below, then this cycle's findings. + +## Journey so far (carried from the archived #30 cycle) + +- **`gdalcubes::filter_geom()` segfaults / all-NA cube** (0.7.3) → cube spans the AOI bbox; clip the + reduced raster with `terra::mask()`. (drift#32) +- **`reduce_time` R-callback runs in spawned workers** where closures fail → moot after the terra pivot. +- **gdalcubes can't read a terra-written NetCDF** → `dft_stac_cube()` returns a terra SpatRaster stack; + `dft_rast_break()` reduces it via forked `parallel::mclapply` (fast: 102k px ~8 s). +- **Sentinel-2 +1000 DN offset flips at 2022-01-25** (baseline 04.00). A uniform offset makes a false + whole-AOI break at 2022 (kNDVI `tanh` hides it; 99% of pixels "broke" at 2022.42). Fixed by a + **baseline-conditional item split** (`offset_boundary`/`offset_before`), coalesced with `terra::cover`. + Validated: per-year kNDVI aligns across 2022; breaks fell to 25% of pixels. +- **Growing-season focus** (`months = 6:9` + snow mask + `order` knob) sharpens signal and cuts data ~3×. +- Shipped v0.3.0, `R CMD check` 0/0/0, PR #33 open. Follow-ups #31 (label breaks), #32 (AOI-polygon clip). + +## Why the redo (review of the shipped vignette) + +- The 2018–2023 / monitor-from-2022 example is dominated by a **region-wide 2023 kNDVI decline** + (~96% of breaks negative, most dated 2023.4–2023.7; per-year mean 2022≈0.53 → 2023≈0.43). Most + plausibly a **2023 dry/wildfire-smoke year** — smoke depresses kNDVI and is NOT an SCL class, so it + leaks through the cloud mask. This swamps localized change. +- The reach is **mostly rangeland with a riparian tree corridor** (from the LULC vignette), so the + broad dip is largely rangeland; **logging = the subset where the green tree corridor drops sharply**. +- User ground truth: the LULC vignette shows **Trees → Rangeland** transitions; a cut is visible at the + **north confluence** changing 2020→2023; **some cuts are 2022–2023** (inside our monitoring window). +- Key method note carried into this cycle: `bfastmonitor` treats everything before `start` as stable + history, so a pre-`start` cut is invisible; 2022–2023 cuts with `start = 2022` ARE catchable. + +## This cycle's findings + +**2017–2023 fetch** (2026-07-09): offset split fired 87 pre / 85 post; 84 monthly layers; +per-year gs kNDVI 2017=0.476 … 2021=0.524, 2022=0.529, 2023=0.429 (aligns across 2022, real 2023 +dip). Cube built in **11 min**, cached in `data-raw/.break_cache`. Breaks: 13,112 finite (vs 25,550 +on the 2018-start run — the longer 2017–2021 history is a more robust baseline, fewer false breaks), +88% negative, dated 2022.42–2023.67. + +**LULC validation** (bfast break map vs IO LULC Trees→Rangeland/Bare, 2020→2023): tree-loss pixels +(n=2132) break at **25% vs 13% background** (2×) and with **median mag −0.054 vs −0.020** (2.7×). +So the trajectory method genuinely picks up LULC tree loss. Break dates on tree-loss pixels cluster +at **median 2023.58 (91% in 2023)** — consistent with "cuts 2022–2023", and the value-add: annual +LULC can only bracket "2020–2023", bfast dates it to summer 2023. + +**Caveats (honest):** the signal is statistical, not a dramatic clearcut. (a) LULC "tree-loss" pixels +sit ~0.1 BELOW intact forest throughout 2017–2021 — IO LULC "Trees" at 10 m is generous (open / +riparian / edge canopy), so they don't start from a high forest baseline; (b) they end ~0.35 in 2023 +(thinned / regrowth, not bare-soil ~0.15); (c) 75% of LULC tree-loss pixels register no bfast break +(below significance in noisy monthly data, or cut pre-2022 = in history). Intact forest (n=4301) stays +~0.5–0.6 and dips only slightly in 2023, so the tree-loss-below-intact divergence IS visible in 2023. + +**Verdict:** the 2017–2023 redo is clearly better than the shipped 2018-start example — more selective +breaks, a real (if modest) LULC agreement, and honest timing. Recommended vignette framing: "the +trajectory method flags the same stands LULC maps as lost, ~2× the rate and stronger, and dates it to +summer 2023 — timing annual LULC cannot give," with the regional 2023 dip acknowledged, not hidden. + +**Reducer decision:** keep `bfastmonitor(start = c(2022,1))` — cuts are 2022–2023 (in-window), and the +2017–2021 history is a strong baseline. Full `bfast()` not needed for this example. + +## Pivot: categorical-vs-continuous reconciliation → the real value (QA + trend) + +User pushed on why bfast's gentle signal disagrees with LULC's dramatic Trees→Rangeland *inside* the +floodplain. Diagnosed with the cached cube: the LULC "tree-loss" pixels had kNDVI baseline **0.43** +(vs 0.50 intact), only **14%** forest-like (>0.5), dropped a median of **0.037** by 2023, only **2%** +a genuine high→low crash. So on this **deciduous-riparian** reach, IO LULC "Trees" is a loose label +and its "Trees→Rangeland" is largely **borderline pixels flipping category on a tiny greenness +change** — the categorical product **overstates** forest loss. In peak summer, deciduous riparian +trees and grass are both green (~0.43), so kNDVI/bfast structurally cannot separate them the way the +structure-based classifier can → the trajectory tools are **weak** on this landscape for tree-vs-grass, +but valuable as a **QA check** on the categorical change and for **gradual degradation/recovery**. + +**Added `dft_rast_trend()`** (Theil-Sen slope + Mann-Kendall p, robust to the 2023 dip). AOI-restricted +trend by group: tree-loss median **~0** (6% sig-declining), intact forest **+0.0026** (2%), open +background **+0.013** (17% recovering). → floodplain is **stable-to-recovering**, and the mapped +"forest loss" is **not backed by a real greenness decline** — an evidence-based "leave it standing" +result and the honest, goal-aligned framing the vignette now uses. + +**Palette bug caught in self-review:** `hcl.colors("Blue-Red 3")` puts red at the HIGH end, so +red = recovering by default; captions said red = declining. Fixed with `rev = TRUE` on both maps so +red = declining/loss (intuitive), consistent across the trend and break maps. diff --git a/planning/active/progress.md b/planning/active/progress.md new file mode 100644 index 0000000..309e85f --- /dev/null +++ b/planning/active/progress.md @@ -0,0 +1,18 @@ +# Progress — Trajectory vignette redo (#30, continued) + +**Continues** `planning/archive/2026-07-issue-30-index-trajectory/progress.md`. The #30 pipeline was +built, validated, and shipped as v0.3.0 (PR #33 open) in the prior cycle — see that archive for the +full log. This cycle refines the trajectory *vignette example* on the same branch / same PR. + +## Session 2026-07-09 + +- Reopened PWF as an explicit continuation of the archived #30 cycle (see task_plan.md / findings.md + headers). Not a new branch or issue — commits land on PR #33. +- Trigger: user review of the shipped vignette. It used 2018–2023 / monitor-from-2022 and the result + was a region-wide 2023 dip (dry/smoke year), not the localized change; the intended scour-vs-stable + contrast was muddy and a known 2022–2023 logging cut didn't clearly show. +- Plan approved: redo the example at **2017–2023** (apples-to-apples with the LULC vignette), validate + the bfast break map against the LULC Trees→Rangeland transition at the north-confluence cut, and + reframe the vignette as LULC (what / roughly when) vs trajectory (exactly when). +- Next: Phase 1 — re-fetch the 2017–2023 growing-season cube (offset split; ~20 min one-time, + then cached). diff --git a/planning/active/task_plan.md b/planning/active/task_plan.md new file mode 100644 index 0000000..df7cdbe --- /dev/null +++ b/planning/active/task_plan.md @@ -0,0 +1,65 @@ +# Task: Trajectory vignette redo — 2017–2023, LULC-comparable (#30, continued) + +**Continuation of** `planning/archive/2026-07-issue-30-index-trajectory/` (v0.3.0 released, +PR #33 open). That cycle built and shipped the pipeline (`dft_stac_cube` / `dft_index_expr` / +`dft_rast_break`) and a first trajectory vignette. On review the vignette example was weak: it +used a 2018–2023 / monitor-from-2022 window, and the result was dominated by a region-wide 2023 +kNDVI dip (dry/smoke year) rather than the localized change we care about — the intended +"scour vs stable" contrast was muddy and a known **logging cut (2022–2023)** did not clearly show. + +**This cycle** refines the *example* (not the pipeline) so the trajectory vignette is apples-to-apples +with the existing LULC vignette (same Neexdzii Kwa reach, **2017–2023**) and actually surfaces the +logging. Commits land on the **same branch / PR #33**. + +## Phase 1: Re-fetch 2017–2023 growing-season cube +- [x] `dft_stac_cube(datetime = "2017-01-01/2023-12-31", months = 6:9)` — offset split fired 87 pre / 85 post; 84 monthly layers; per-year kNDVI aligns across 2022. Built in **11 min**, cached + +## Phase 2: Reduce + reducer decision +- [x] `dft_rast_break(start = c(2022, 1))`, 2017–2021 history → 13,112 finite breaks (vs 25,550 on the 2018 start — longer baseline halves false breaks) +- [x] Kept `bfastmonitor` — cuts are 2022–2023 (in-window); full `bfast()` not needed + +## Phase 3: Validate against LULC ground truth +- [x] LULC Trees→Rangeland/Bare vs bfast breaks: tree-loss pixels break **25% vs 13%** (2×), median mag **−0.054 vs −0.020** (2.7×), dated **median 2023.58** (91% in 2023). Real, modest signal +- [x] Quantified: LULC brackets "2020–2023"; trajectory dates it to summer 2023. Caveats logged (IO LULC "Trees" generous; regional 2023 dip; 75% of tree-loss pixels no break) + +## Phase 4: Rebuild the vignette honestly +- [x] `data-raw/vignette_data_break.R`: 2017–2023 window; precomputes break raster + LULC tree-loss mask + grouped trajectory (tree-loss/intact/background) + agreement stats → 140 KB artifact +- [x] `vignettes/trajectory-break-detection.Rmd`: reframed as complement to LULC (what/roughly-when vs exactly-when); agreement table; break map w/ tree-loss outlined; grouped trajectory; regional-2023-dip acknowledged. **User agreed to framing + trajectory panel** +- [x] Regenerated `inst/testdata/neexdzii_break.rds` + +## Phase 5: Verify + land on PR #33 +- [x] Vignette renders + re-builds clean under `R CMD check` (only a spurious future-timestamp NOTE); `devtools::test()` 286 pass; `lintr` clean on the two new files +- [x] Atomic commit on branch `30-continuous-index-trajectory-change-detec` (joins PR #33) +- [ ] `/planning-archive` (append to / supersede the existing #30 archive) + +## Phase 6: Reframe as QA + degradation/recovery (goal-aligned rebuild) + +Investigating why bfast's gentle signal disagrees with LULC's dramatic Trees→Rangeland +revealed the real value: LULC "tree-loss" pixels had kNDVI baseline 0.43 (vs 0.50 intact), +dropped only 0.037 by 2023, only 2% a real crash. So on this deciduous-riparian reach IO LULC +likely **overstates** forest loss (borderline label flips), and the continuous trajectory's job is +to **question/QA** the categorical change and catch **gradual degradation/recovery** the annual +labels miss — which is exactly the project goal (justify leave-standing / restoration). User +approved the reframe. + +- [x] `R/dft_rast_trend.R`: new export — Theil-Sen slope + Mann-Kendall p, mclapply, `.dft_trend_pixel` helper. 16 tests pass, lint clean +- [x] NEWS 0.3.0 updated to list `dft_rast_trend`; docs regenerated +- [x] Rebuilt `vignettes/trajectory-break-detection.Rmd` around QA + degradation/recovery. **QA result (AOI): LULC tree-loss median trend ~0, only ~6% significantly declining (vs ~2% elsewhere); open floodplain greening (+0.013/yr, 17% recovering); intact forest stable.** So the mapped "forest loss" is mostly not real greenness decline — caution on area numbers. Fixed an inverted trend-map palette (red must = declining) +- [x] Generator precomputes break + trend + LULC-QA table + grouped trajectory (AOI-restricted groups) → artifact +- [x] Render + `R CMD check` (0/0/1 spurious-timestamp) + `lintr` clean + `/code-check` (stats verified correct; added nlyr>=2 + min_obs>=2 guards for clean failure); commit on PR #33 + +## Phase 7: Interactive map + patch-level reframe (user-driven) + +Reconciling with the user's own imagery: 2 real cutblocks in the floodplain align with red +patches in the trend/break maps → the tools DO catch real harvest in spots; some red outside the +LULC polygons = removal LULC missed → complementary. User steered: drop whole-floodplain averaging, +add an interactive map (as close to the LULC vignette's as possible), scientific (non-advocacy) tone. + +- [x] Generator: LULC Trees→Rangeland transition polygons via `dft_rast_transition` + `dft_transition_vectors` (filtered to the actual loss transition, NOT stable "Trees→Trees" which looked like the whole AOI); dropped the aggregate QA table + 3-group averaged trajectory + an auto cut-patch trajectory (couldn't reliably auto-pick a clean cutblock — interactive map is the honest verifier) +- [x] Vignette: rebuilt around patchy/complementary framing (red in/out/absent vs LULC outline); interactive leaflet mirroring `dft_map_interactive` basemaps (Light/Esri/Google/OpenTopoMap) + centering, with break/trend rasters + beige (#e3e2c3) Trees→Rangeland polygons + AOI, toggleable/fullscreen; trend relabeled "per year 2017-2023"; opening paragraph reworded to a neutral methods-first tone +- [x] `/code-check` Clean (contract/CRS/leaflet verified; removed an unused `aoi_v`) + `R CMD check` 0/0/1-spurious + lint clean + commit on PR #33 + +## Validation +- [x] Tests pass (302); check clean (0 err / 0 warn / 1 spurious NOTE) +- [x] Commits on PR #33, not a new branch +- [ ] `/planning-archive` on completion diff --git a/planning/archive/2026-07-issue-30-index-trajectory/.gitkeep b/planning/archive/2026-07-issue-30-index-trajectory/.gitkeep new file mode 100644 index 0000000..e69de29 diff --git a/planning/archive/2026-07-issue-30-index-trajectory/README.md b/planning/archive/2026-07-issue-30-index-trajectory/README.md new file mode 100644 index 0000000..eb1e8fc --- /dev/null +++ b/planning/archive/2026-07-issue-30-index-trajectory/README.md @@ -0,0 +1,45 @@ +# Issue #30 — Continuous index-trajectory change detection (Sentinel-2 + BFAST) + +## Outcome + +Added a continuous change-detection pipeline alongside the categorical one, released as **v0.3.0** +(`R CMD check`: 0 errors / 0 warnings / 0 notes). Four new exports: `dft_stac_cube()` (fetch a +cloud-masked monthly spectral-index stack from Sentinel-2 via a new role-based `sentinel-2-l2a` +config source), `dft_index_expr()` / `dft_index_table()` (a table-driven index registry — NDVI, +kNDVI, NDMI — written over band roles), and `dft_rast_break()` (per-pixel `bfast::bfastmonitor()` +reduction to a two-band break-date / break-magnitude raster). A network-decoupled vignette +(`trajectory-break-detection.Rmd`) demonstrates it on the same Neexdzii Kwa reach as the +land-cover vignette, loading a small committed artifact (132 KB) generated by +`data-raw/vignette_data_break.R`. + +## What the empirical work changed vs. the original plan + +The plan was sound but three things could only be settled by running real code, and each changed +the implementation materially: + +1. **The reduce_time R-callback runs in spawned gdalcubes workers at every parallel setting**, and a + closure over enclosing locals fails there (`object 'band' not found`). First fixed with a + self-contained embed-the-helper callback; ultimately mooted by the terra pivot below. +2. **`gdalcubes::filter_geom()` segfaults / yields an all-NA cube** on the pinned build. Dropped it; + the cube spans the AOI bbox and callers clip the reduced raster (matches `dft_stac_fetch`). The + AOI-polygon clip is deferred to #32. +3. **The Planetary Computer +1000 DN reflectance offset only applies from 2022-01-25** (baseline + 04.00); a uniform offset produced a *false whole-AOI break* at the boundary (99% of pixels broke + at 2022.42 — an artifact `tanh`-kNDVI hid). Fixed with a baseline-conditional item split + (`offset_boundary` / `offset_before` in config), coalesced with `terra::cover`. Validated: per-year + kNDVI now aligns across 2022 and breaks fell to 25% of pixels (real change). A peer validation + session independently confirmed the offset finding. + +**Architecture pivot:** because gdalcubes cannot read a terra-written NetCDF, the coalesce had to +live in terra — so `dft_stac_cube()` now returns a **terra SpatRaster stack** and `dft_rast_break()` +reduces it via forked `parallel::mclapply` (simpler and faster — 102k px in 8.3s — than the +gdalcubes-worker path, which was removed). User steered two decisions: growing-season focus +(`months` filter + snow mask + `order` knob; better signal, ~3× less data) and fixing the offset +properly now rather than deferring. + +Follow-ups filed: **#31** (label breaks with from/to land-cover classes, depends on #19), **#32** +(restore AOI-polygon clip once the gdalcubes `filter_geom` issue is resolved). The offset half of the +original #32 was closed here. + +Closed by: branch `30-continuous-index-trajectory-change-detec` (11 commits, `6da759c`..`c145899`); +PR pending. diff --git a/planning/archive/2026-07-issue-30-index-trajectory/findings.md b/planning/archive/2026-07-issue-30-index-trajectory/findings.md new file mode 100644 index 0000000..40592ac --- /dev/null +++ b/planning/archive/2026-07-issue-30-index-trajectory/findings.md @@ -0,0 +1,165 @@ +# Findings — Continuous index-trajectory change detection (Sentinel-2 + BFAST) (#30) + +## Design validation (session 2026-07-08, plan mode) + +Verified directly against the installed toolchain (read-only): gdalcubes 0.7.3, +rstac 1.0.1, terra 1.9.11, zoo 1.8.15. `bfast` / `strucchangeRcpp` NOT installed. + +### gdalcubes `reduce_time` R-callback contract (from `reduce_time.cube.Rd`) +- Signature: `reduce_time(x, expr, ..., FUN, names = NULL, load_pkgs = FALSE, load_env = FALSE)`. +- FUN receives a **2-D array: rows = bands, cols = time steps**. Row names ARE + band names — the shipped example indexes `x["B04", ]`, so `x["kndvi", ]` works + **iff the input cube's band is literally named `kndvi`**. +- FUN must return a numeric vector of length `length(names)` → those become the + output bands. +- **FUN runs in spawned worker R processes** (`load_pkgs`/`load_env`). Therefore: + - `bfast` must be loaded in workers → `reduce_time(..., load_pkgs = "bfast")`. + The main-process `rlang::check_installed("bfast")` gate alone is insufficient. + - `band`/`ts_start`/`frequency`/`start`/`level`/`min_obs` must be + **closure-captured scalars**, never derived from `cube` inside FUN. +- gdalcubes ships `L8NY18` (228 local Landsat TIFs at + `system.file("L8NY18", package="gdalcubes")`) → an **offline** cube fixture for + reduce_time / apply_pixel shape tests (no network). +- `create_image_collection(files, date_time=, band_names=)` builds a collection + from local GeoTIFFs deterministically → synthetic offline cube feasible. + +### Corrections to the issue-#30 sketch (all confirmed) +1. `tanh(ndvi^2)` → gdalcubes `apply_pixel` uses the **tinyexpr C engine**; `^` + is not R exponentiation. Use `tanh(pow((nir-red)/(nir+red), 2))`. tinyexpr + provides `pow`, `tanh`, `sqrt`, `exp`, `log`, etc. +2. `cube_view(extent = col)` sizes the cube to the union bbox of all returned + scenes (full S2 tile), then `filter_geom` masks most to NaN AFTER the work. + Build `extent` explicitly from `bbox_target` + parsed `t0`/`t1`, as + `dft_stac_fetch()` already does. +3. `ext_filter(\`eo:cloud_cover\` <= cloud_cover_max)` uses NSE — a runtime + variable needs `{{cloud_cover_max}}` and the property back-ticked. spacehakr's + example only ever used a literal `<= 10`, so it never hit this. +4. `start_of(cube)` and `terra_from_cube()` are pseudo-code. Real: ts start from + `dimension_values(cube, "M")$t` (main process); cube→terra via + `write_ncdf(overwrite=TRUE)` → `terra::rast()` → `terra::mask()` (repo idiom). +5. `stac_search` cannot take both `bbox` and `intersects` — cube mode uses + `intersects` only. +6. Config restructure is **forced backward-compatible** by the existing + `expect_named(cfg, c("stac_url","collection","asset","available_years"))` + test (exact-set). Cube sources add `cube = TRUE`; categorical sources keep + exactly 4 fields → marker is **absence = FALSE** (`isTRUE(cfg$cube)`). +7. Scale/offset must be applied **inside** the index expression as a per-band + affine token `(asset*scale+offset)`. A non-zero offset does NOT cancel in a + ratio index → NDVI-on-raw-DNs is wrong for Landsat (offset -0.2) and + baseline-04.00 S2 (offset -0.1). Only a pure multiplicative scale cancels. + +### API names verified present with the used signatures +`image_mask(band, min, max, values, bits, invert)`, `filter_geom(cube, geom, srs)`, +`apply_pixel(x, ...)` (generic; `expr`/`names`), `raster_cube(..., mask=)`, +`cube_view`, `stac_image_collection`, `write_ncdf`, `ncdf_cube`, +`stac_search(..., intersects=, bbox=, limit=)`, `ext_filter(q, expr)`, +`post_request`, `items_fetch`, `items_sign(items, sign_fn)`, +`sign_planetary_computer()`, `dimension_values`. `bands(cube)` → data.frame with +`name/type/scale/offset/unit`. + +### bfast::bfastmonitor return / semantics (knowledge; bfast not installed) +- `$breakpoint` = **decimal year** (e.g. 2022.34) or `NA` when no break (NOT an + error). `$magnitude` = median of monitoring-period residuals (observed − + predicted) → **negative = index drop = scour / veg loss; positive = + establishment**. Inherent to the method; must be documented for callers who + threshold on it. +- `start = c(year, period)` = monitoring-period start; `history=`, `level=` pass + through. Failure modes to `tryCatch`: too-few history obs for + `response ~ trend + harmon`, singular fits, gappy series → return `c(NA, NA)`. + +### Open correctness subtlety (documented, follow-up candidate) +- S2 baseline-04.00 (+1000 DN) offset is **not uniform** across a 2017→2024 + archive: scenes processed under baseline ≥ 04.00 (~≥ 2022-01-25) carry the + offset; earlier scenes have offset 0. A single `offset=-0.1` is slightly wrong + for the pre-2022 portion. kNDVI is fairly robust (scale cancels in the ratio; + residual is a small additive bias near the 2022 boundary). Ship documented + accepted-bias; per-scene harmonization from item metadata is a follow-up. + +## Phase 4 empirical gate result (2026-07-08) — DESIGN CHANGE + +Ran the smoke test after installing bfast 1.7.2 / strucchangeRcpp 1.5.4. Findings: +- **Leg (a) passed exactly**: a 72-pt monthly series with a −0.3 step at 2022-06 + → `bfastmonitor(start=c(2022,1), level=0.01)` returns `breakpoint=2022.417`, + `magnitude=-0.300`. Scalar `level` is accepted (default is length-2). Stable + series → `breakpoint=NA` (no error). So `history="all"` + scalar `level` work. +- **CLOSURE CAPTURE DOES NOT WORK.** The reducer FUN runs in a spawned worker for + the R-callback path **at every `parallel` setting, including `parallel=1`**. + A closure (even a `force()`d factory frame) fails in the worker with + `object 'band' not found` — gdalcubes does not restore the closure environment. + Both the issue-#30 sketch and the original agent-B design assumed closures work; + they don't. This is the gate's headline catch. +- **SELF-CONTAINED FUN WORKS** at `parallel` 1 and 2. Two variants proven: + inline scalar literals via `substitute()`, or embed the per-pixel helper as a + literal function OBJECT (env detached to `baseenv()`) + inline scalar params. + Chosen: **embed-the-helper-object** — keeps a single logic source + (`.dft_break_pixel()`, directly unit-testable) AND a worker-safe FUN. + Result on the 2×2 synthetic cube: `nlyr=2`, names `break_date`/`break_mag`, + drop pixel `date=2022.42 mag=-0.3`, stable pixels `NA/0`, all-NA pixel `NA/NA`. +- `load_pkgs = "bfast"` is required so the worker can resolve `bfast::bfastmonitor`. +- Degenerate paths (`all(is.na(v))` / `sum(!is.na(v)) < min_obs`) return `c(NA,NA)` + **before any `bfast::` symbol** → unit-testable with no bfast dependency. + +Design for `R/dft_rast_break.R` (proven): +```r +.dft_break_pixel <- function(v, ts_start, frequency, start, history, 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, level = level) + c(m$breakpoint, m$magnitude) + }, error = function(e) c(NA_real_, NA_real_)) +} +# reducer FUN: embed .dft_break_pixel as a literal object (env=baseenv) + inline +# band/ts_start/frequency/start/history/level/min_obs via substitute(); no free +# vars; reduce_time(..., load_pkgs = "bfast"). +``` +`bfastmonitor` return: `$breakpoint` = decimal year or NA; `$magnitude` = median +monitoring-period residual (neg = index drop = scour; ~0 = stable). `ts_start` +and `frequency` come from `gdalcubes::dimension_values(cube, "M")$t` in the main +process; `stop()` if a caller `frequency` disagrees with the dt cadence. + +## Sentinel-2 baseline offset — the correctness catch (2026-07-08) + +A peer validation session + my own growing-season E2E stats surfaced that PC's +`sentinel-2-l2a` +1000 DN reflectance offset only applies from processing +baseline 04.00 (2022-01-25 on); earlier scenes have offset 0. A **uniform** +`offset = -0.1` is therefore wrong pre-2022 and creates a false whole-AOI index +step at the boundary. kNDVI's `tanh` hid it (bounded 0-1) so the cube looked +valid, but my E2E had **90903/91467 (99%) negative breaks all at 2022.42** — the +boundary, not vegetation. User chose to fix it properly (baseline-conditional +split), not defer. + +**terra re-architecture (validated):** +- gdalcubes CANNOT read a terra-written NetCDF (dim-name mismatch), so coalescing + pre/post cubes at the gdalcubes level is impossible. Pivoted to terra. +- `dft_stac_cube()` returns a **terra SpatRaster stack** (materialized GeoTIFF, + `terra::time` set). If the source has `offset_boundary` and items straddle it, + it builds a pre cube (offset `offset_before`) and post cube (offset) over the + SAME full `cube_view`, and coalesces with `terra::cover(pre, post)`. +- `dft_rast_break()` reduces the stack via `parallel::mclapply` (fork -> closures + and package internals inherited, no gdalcubes-worker serialization). Validated: + 102400 px in 8.3 s. `.dft_break_pixel` unchanged. Dropped build_break_reducer, + break_cache_key, gdalcubes reduce_time. +- `cadence_frequency()` now derives frequency from the median day-spacing of the + layer times (26-32d -> 12, 85-95d -> 4, 360-370d -> 1). + +**Split validation E2E (2018-2023, months=6:9):** split fired 75 pre / 85 post; +per-year growing-season mean kNDVI now ALIGNS across the boundary +(2021=0.524 -> 2022=0.529, all years 0.41-0.53, no jump); finite breaks dropped +from 99% to **25% of usable pixels**, spread across 2022-2024. Offset artifact +gone. Config: `sentinel-2-l2a` gains `offset_boundary="2022-01-25"`, +`offset_before=0`. + +## Reuse map +- `stac_cache_key()` / `auto_utm_epsg()` / sf-coercion / `write_ncdf→terra→mask` + idiom in `R/dft_stac_fetch.R` — reuse for `dft_stac_cube()`. +- `dft_class_table()` CSV-reader in `R/dft_stac_classes.R` — mirror for + `dft_index_table()` (`inst/lulc_classes/` → `inst/indices/`). +- `spacehakr::spk_stac_calc()` — pagination/intersects/cloud-filter proof and QA + chip reader; used ONLY in `data-raw/`, never a DESCRIPTION dependency. + +## Issue context + +Full body: `gh issue view 30 --repo NewGraphEnvironment/drift`. +Relates to #9 (deferred BFAST evaluation), #19 (arbitrary factor rasters — +the optional from-to labelling bolt-on depends on it). diff --git a/planning/archive/2026-07-issue-30-index-trajectory/progress.md b/planning/archive/2026-07-issue-30-index-trajectory/progress.md new file mode 100644 index 0000000..504f30e --- /dev/null +++ b/planning/archive/2026-07-issue-30-index-trajectory/progress.md @@ -0,0 +1,61 @@ +# Progress — Continuous index-trajectory change detection (Sentinel-2 + BFAST) (#30) + +## Session 2026-07-08 + +- Plan-mode exploration — two Plan agents vetted the design against issue #30 and + the codebase (acquisition half: config + index registry + cube; reduction half: + bfast reducer + tests + vignette). Both independently confirmed the issue-sketch + corrections. Phases approved by user (one branch, single 0.3.0 release). +- Empirically verified the gdalcubes/rstac/terra APIs read-only (see findings.md). +- Created branch `30-continuous-index-trajectory-change-detec` off main. +- Scaffolded PWF baseline from issue #30 with approved phases. +- **Phase 1 done** (config role-based schema + sentinel-2-l2a + fetch guard + tests). + Fresh-eyes code-check caught a dangling `[dft_stac_cube()]` roxygen link (R CMD + check WARNING) → downgraded to code spans until Phase 3; hardened the fetch guard + to key on `isTRUE(cfg$cube)` alone. 27 config + 16 fetch assertions green. +- **Phase 2 done** (index registry `inst/indices/indices.csv` + `dft_index_expr()` + + `dft_index_table()` + resolver internals; 15 assertions green; code-check Clean). +- **Phase 3 done** (`dft_stac_cube()` + `stac_cube_cache_key()`; 20 assertions green). + Code-check caught a multi-feature-AOI `intersects` bug (silent NoData holes) → + fixed by unioning the AOI. `eo:cloud_cover` NSE declared in globalVariables. +- Committed Phases 2+3 together (tightly-coupled acquisition pipeline; roxygen + cross-refs only resolve once both exist). Full suite: 261 pass / 0 fail / 3 skip. +- **Phase 4 empirical gate RUN** (bfast 1.7.2 installed): proved closure-capture in + the reduce_time R-callback DOES NOT WORK (worker process, any parallel setting) — + the design must build a self-contained FUN. Embed-helper-object approach works at + parallel 1 and 2. See findings.md. +- **Phase 4 done** (`dft_rast_break()` + `.dft_break_pixel`/`cadence_frequency`/ + `build_break_reducer`/`break_cache_key`; bfast→Suggests). 30 assertions green; + full suite 291 pass. Code-check Clean. +- **Real-S2 E2E surfaced a `dft_stac_cube` bug**: `gdalcubes::filter_geom()` yields + an all-NA cube (and intermittently segfaults the compute worker) on this gdalcubes + build. Proven sibling `dft_stac_fetch` works on the same AOI, and the cube stages + work WITHOUT filter_geom (valid kNDVI 0..1). Fix: drop filter_geom; the cube spans + the AOI bbox (cloud-masked), callers clip the reduced raster with terra::mask() + (matches how dft_stac_fetch masks). Fixed E2E verified (48-mo cube valid, break + reduce 25 s, no segfault); committed (c6953a0). +- **User steer — growing-season tuning.** All-months fetch was ~30 min for 100 ha + (the fetch, not the compute). User asked whether restricting to the growing season + sharpens the veg signal AND cuts data — yes on both. Added `dft_stac_cube(months=)` + calendar-month filter (default NULL); snow (SCL 11) now in the default S2 mask; GDAL + /vsicurl config (on.exit restore); `dft_rast_break(order=)` harmonic knob. Growing- + season E2E (months=6:9, 2018-2023): cube in 12 min (vs 30), 24 clean summer obs/pixel, + breaks dated to summer months (2022.42-2023.67). 294 unit tests pass, lint clean. +- Committed growing-season enhancement (61eaa5d); code-check Clean. +- **Peer session + my E2E stats surfaced a real offset bug.** PC's uniform offset=-0.1 + is wrong for pre-2022-01-25 scenes (the +1000 DN baseline offset flips there). A + series crossing the boundary shows a false whole-AOI break: my own growing-season + E2E had 90903/91467 (99%) negative breaks all at 2022.42 — the offset boundary, not + vegetation. User chose to fix it PROPERLY now (baseline-conditional split). +- **Re-architecture (terra route).** gdalcubes can't read a terra-written NetCDF + (round-trip fails), so coalescing pre/post cubes at the gdalcubes level is out. + Pivoted the reduction to terra: dft_stac_cube now returns a terra SpatRaster stack + (materialized GeoTIFF, time set), splitting items at the offset boundary and + correcting each side with its own offset (terra::cover coalesce). dft_rast_break + reduces the stack via parallel::mclapply (fork -> closures + package internals work, + no gdalcubes-worker serialization; validated 102400 px in 8.3 s). Dropped + build_break_reducer/break_cache_key/gdalcubes reduce_time; .dft_break_pixel unchanged. + Config: sentinel-2-l2a gains offset_boundary="2022-01-25" + offset_before=0. Full + suite 286 pass, lint clean. +- Next: confirm the split kills the fake 2022 step (per-year kNDVI aligns; breaks + sparse), regenerate artifact, Phase 5 (vignette + NEWS + 0.3.0 + archive + PR). diff --git a/planning/archive/2026-07-issue-30-index-trajectory/task_plan.md b/planning/archive/2026-07-issue-30-index-trajectory/task_plan.md new file mode 100644 index 0000000..f2b3f8d --- /dev/null +++ b/planning/archive/2026-07-issue-30-index-trajectory/task_plan.md @@ -0,0 +1,50 @@ +# Task: Continuous index-trajectory change detection (Sentinel-2 + BFAST) for floodplain reaches (#30) + +drift detects change today by differencing pre-classified annual land-cover maps +(`dft_stac_fetch` → `dft_rast_classify` → `dft_rast_transition`). Masked to a +floodplain polygon that is the *weakest* place for the Esri IO LULC 9-class +annual product: it collapses gravel bars / wet-dry channel / sparse seasonal +riparian into `bare`/`rangeland`, and annual compositing smears event timing. + +Issue #30 adds the stronger signal: a **continuous index trajectory** per pixel — +detect *when* a pixel's spectral history breaks and by how much — instead of +comparing two categorical labels. The machinery is already a dependency +(`gdalcubes`, `rstac`, `terra`); this is a new fetch mode + a reducer, not a new +toolchain. Four new exports: `dft_stac_cube()`, `dft_index_expr()` +(+ `dft_index_table()`), `dft_rast_break()`. Plus a config restructure and a +network-decoupled vignette. **Whole pipeline on one branch, single 0.3.0 release.** + +## Phase 1: Config — role-based schema + Sentinel-2 source +- [x] `R/dft_stac_config.R`: expand `match.arg` to `c("io-lulc","esa-worldcover","sentinel-2-l2a")`; add `sentinel-2-l2a` entry (`cube=TRUE`, `roles=list(red="B04",nir="B08",swir16="B11",mask="SCL")`, `mask_values=c(3L,8L,9L,10L)`, `scale=1e-4`, `offset=-0.1`, `available_datetime`); keep `io-lulc`/`esa-worldcover` 4-field **verbatim**; add `landsat-c2-l2` as a commented drop-in template (`red`/`nir08`, `qa_pixel`, scale 2.75e-5 / offset -0.2). NOTE: `landsat-c2-l2` left OUT of `match.arg` (in match.arg but no switch case = silent NULL); it stays a commented template until implemented +- [x] `R/dft_stac_fetch.R`: friendly guard in config-resolve block — if `isTRUE(cfg$cube)`, `cli::cli_abort(... use dft_stac_cube())`; keyed on `cube` alone (not `&& is.null(asset)`) so a cube source never falls through to an opaque `stopifnot(years)`; categorical path untouched +- [x] `tests/testthat/test-dft_stac_config.R`: keep the 3 existing tests byte-for-byte; add cube-source assertions (`isTRUE(cfg$cube)`, role names + `$red=="B04"`/`$nir=="B08"`/`$mask=="SCL"`, `mask_values`, `scale`, `offset`, `available_datetime`); backward-compat assertion that `io-lulc`/`esa-worldcover` have no `cube` field + +## Phase 2: `dft_index_expr()` + table-driven index registry +- [x] `inst/indices/indices.csv`: registry (`index,formula,roles,description`) — `ndvi`, `kndvi=tanh(pow((nir-red)/(nir+red),2))`, `ndmi`; formulas in tinyexpr C syntax over **roles**, self-contained (kNDVI inlines NDVI). Uses `pow()` not `^` (gdalcubes tinyexpr) +- [x] `R/dft_index_expr.R`: `dft_index_table()` (reads the CSV, mirrors `dft_class_table()`); internal `index_resolve_expr(index, roles, scale, offset)` — pure string builder, per-role scaled token `(asset*scale±|offset|)` when scale≠1 or offset≠0 else bare asset, longest-role-first word-boundary substitution, unknown-index/absent-role `cli_abort`; internals `index_row()`/`index_roles()`/`scale_token()`; exported `dft_index_expr(cube, index="kndvi", source="sentinel-2-l2a", roles=NULL, scale=NULL, offset=NULL)` → `gdalcubes::apply_pixel(cube, expr, names = index)` +- [x] `tests/testthat/test-dft_index_expr.R` (pure local, no gdalcubes/network): kNDVI-over-S2 resolves to the exact scaled string; scale/offset on (Landsat) vs off (bare) tokens; `ndmi` resolves `swir16→B11`; unknown-index + absent-role errors. 15 assertions green. **Code-check Clean.** + +## Phase 3: `dft_stac_cube()` +- [x] `R/dft_stac_cube.R`: full signature per plan; `check_installed("gdalcubes")`; `!isTRUE(cfg$cube)` reject; reuse `auto_utm_epsg`/sf-coercion/`dft_cache_path`; derive assets from `index_roles(index)` + `cfg$roles$mask`; STAC query with **`intersects` (unioned AOI)** + `ext_filter(\`eo:cloud_cover\` <= {{cloud_cover_max}})` + `post_request() |> items_fetch() |> items_sign()`; **explicit** `cube_view` extent from `bbox_target` + parsed `t0`/`t1`; `raster_cube(col, v, mask=image_mask(mask_asset, values=mask_values)) |> filter_geom(...)`; `dft_index_expr(...)`; materialize to `/cube_.nc` via `write_ncdf(overwrite=TRUE)`, return `gdalcubes::ncdf_cube(f)`; `force` overwrites. **Code-check fixed a multi-feature-AOI bug: union AOI before `intersects` (else silent NoData holes — the #25 class).** `eo:cloud_cover` NSE declared in `utils::globalVariables()` +- [x] Internal `stac_cube_cache_key()` (leaves `stac_cache_key()` byte-for-byte unchanged): hash AOI-WKB + `res`/`target_crs`/`dt`/`aggregation`/`resampling`/`stac_url`/`collection` + cube additions `band_assets`/`datetime`/`index`/`cloud_cover_max`/`sort(mask_values)`/`scale`/`offset` +- [x] `tests/testthat/test-dft_stac_cube.R`: local gdalcubes-missing guard; categorical-source reject; `stac_cube_cache_key` unit tests (determinism/12-hex, per-param sensitivity, mask-order + res-type normalization, sf-attrs ignored); network E2E gated behind `DRIFT_TEST_NETWORK=true` env var (keeps `devtools::test()` network-free per repo convention). 20 assertions green + +## Phase 4: `dft_rast_break()` — bfast reducer (empirical gate FIRST) +- [x] **Smoke-test gate** (bfast 1.7.2 installed): legs (a) `bfastmonitor` on a synthetic step-drop → `breakpoint=2022.417 magnitude=-0.30`; (b) synthetic local cube through `reduce_time` → 2-band raster, correct names, break/NA pixels correct; (c) `frequency=1` on a P1M cube yields different (garbage) → guard justified. **KEY FINDING: closure capture DOES NOT WORK in the reduce_time R-callback (worker process, any parallel setting) — `object 'band' not found`.** Design changed to a self-contained callback (embed helper object + inline literals), proven at parallel 1 and 2. See findings.md. +- [x] `R/dft_rast_break.R`: internal `.dft_break_pixel(v, ts_start, frequency, start, history, level, min_obs)` — degenerate all-NA / ` its asset", { + expr <- drift:::index_resolve_expr("ndmi", s2_roles, scale = 1, offset = 0) + expect_equal(expr, "(B08 - B11) / (B08 + B11)") +}) + +test_that("index_roles reports the roles an index needs", { + expect_setequal(drift:::index_roles("kndvi"), c("nir", "red")) + expect_setequal(drift:::index_roles("ndmi"), c("nir", "swir16")) +}) + +test_that("unknown index errors with the available set", { + expect_error(drift:::index_resolve_expr("bogus", s2_roles), "Unknown index") + expect_error(drift:::index_roles("bogus"), "Unknown index") +}) + +test_that("an index needing an absent role errors", { + # ndmi needs swir16; a role map without it must fail loudly + expect_error( + drift:::index_resolve_expr("ndmi", list(red = "B04", nir = "B08")), + "role" + ) +}) diff --git a/tests/testthat/test-dft_rast_break.R b/tests/testthat/test-dft_rast_break.R new file mode 100644 index 0000000..85214fc --- /dev/null +++ b/tests/testthat/test-dft_rast_break.R @@ -0,0 +1,97 @@ +# ---- guards ---------------------------------------------------------------- + +test_that("dft_rast_break requires bfast", { + skip_if(requireNamespace("bfast", quietly = TRUE), + "bfast is installed, can't test missing-package path") + expect_error(dft_rast_break(terra::rast(nrows = 1, ncols = 1)), "bfast") +}) + +test_that("dft_rast_break errors on non-SpatRaster input", { + skip_if_not_installed("bfast") + expect_error(dft_rast_break(list()), "SpatRaster") +}) + +test_that("dft_rast_break errors when layers lack times", { + skip_if_not_installed("bfast") + r <- terra::rast(nrows = 2, ncols = 2, nlyrs = 3) + terra::values(r) <- 1 + expect_error(dft_rast_break(r), "time") +}) + +# ---- cadence -> frequency (pure) ------------------------------------------- + +test_that("cadence_frequency maps layer times to seasonal frequency", { + monthly <- seq(as.Date("2020-01-01"), by = "month", length.out = 12) + quarterly <- seq(as.Date("2020-01-01"), by = "3 months", length.out = 8) + annual <- seq(as.Date("2018-01-01"), by = "year", length.out = 5) + expect_equal(drift:::cadence_frequency(monthly), 12) + expect_equal(drift:::cadence_frequency(quarterly), 4) + expect_equal(drift:::cadence_frequency(annual), 1) + expect_true(is.na(drift:::cadence_frequency(as.Date("2020-01-01")))) # length 1 +}) + +# ---- .dft_break_pixel degenerate paths (no bfast needed) ------------------- + +test_that(".dft_break_pixel returns c(NA, NA) on all-NA input", { + expect_equal( + drift:::.dft_break_pixel(rep(NA_real_, 24), c(2018, 1), 12, c(2022, 1), + "all", 3, 0.01, 6), + c(NA_real_, NA_real_) + ) +}) + +test_that(".dft_break_pixel returns c(NA, NA) when fewer than min_obs", { + v <- c(0.5, 0.6, NA, NA, 0.55, rep(NA_real_, 19)) # 3 non-NA < min_obs 6 + expect_equal( + drift:::.dft_break_pixel(v, c(2018, 1), 12, c(2022, 1), "all", 3, 0.01, 6), + c(NA_real_, NA_real_) + ) +}) + +# ---- bfast-gated per-pixel behavior ---------------------------------------- + +test_that(".dft_break_pixel detects an injected step drop (negative magnitude)", { + skip_if_not_installed("bfast") + tt <- seq_len(72) + v <- 0.6 + 0.15 * sin(2 * pi * tt / 12) + v[54:72] <- v[54:72] - 0.3 # step drop at 2022-06 (monitoring period) + out <- drift:::.dft_break_pixel(v, c(2018, 1), 12, c(2022, 1), "all", 1, 0.01, 6) + expect_true(is.finite(out[1])) + expect_gt(out[1], 2022) + expect_lt(out[1], 2022.9) + expect_lt(out[2], 0) +}) + +test_that(".dft_break_pixel returns NA break on a stable series (non-error)", { + skip_if_not_installed("bfast") + tt <- seq_len(72) + v <- 0.6 + 0.15 * sin(2 * pi * tt / 12) + out <- drift:::.dft_break_pixel(v, c(2018, 1), 12, c(2022, 1), "all", 1, 0.01, 6) + expect_true(is.na(out[1])) +}) + +# ---- reduction over a stack (terra + bfast, no network) -------------------- + +test_that("dft_rast_break reduces a synthetic stack to a 2-band raster", { + skip_if_not_installed("bfast") + stk <- synthetic_break_stack() + breaks <- dft_rast_break(stk, start = c(2022, 1), order = 1, cores = 2) + + expect_s4_class(breaks, "SpatRaster") + expect_equal(terra::nlyr(breaks), 2) + expect_equal(names(breaks), c("break_date", "break_mag")) + + bd <- terra::values(breaks[["break_date"]])[, 1] + bm <- terra::values(breaks[["break_mag"]])[, 1] + # the engineered drop pixel yields a finite break with negative magnitude; + # the all-NA pixel yields NA + expect_true(any(is.finite(bd))) + expect_true(any(is.na(bd))) + expect_true(any(is.finite(bd) & bm < 0)) +}) + +test_that("dft_rast_break errors when frequency disagrees with the cadence", { + skip_if_not_installed("bfast") + stk <- synthetic_break_stack() # monthly -> cadence 12 + expect_error(dft_rast_break(stk, frequency = 1), "frequency") +}) diff --git a/tests/testthat/test-dft_rast_trend.R b/tests/testthat/test-dft_rast_trend.R new file mode 100644 index 0000000..2ee0b61 --- /dev/null +++ b/tests/testthat/test-dft_rast_trend.R @@ -0,0 +1,66 @@ +# monthly decimal-year time axis for the pure-helper tests +t24 <- seq_len(24) / 12 + +test_that("dft_rast_trend errors on non-SpatRaster input", { + expect_error(dft_rast_trend(list()), "SpatRaster") +}) + +test_that("dft_rast_trend errors when layers lack times", { + r <- terra::rast(nrows = 2, ncols = 2, nlyrs = 3) + terra::values(r) <- 1 + expect_error(dft_rast_trend(r), "time") +}) + +# ---- .dft_trend_pixel degenerate paths (no stats needed) ------------------- + +test_that(".dft_trend_pixel returns c(NA, NA) on all-NA / too-short input", { + expect_equal(drift:::.dft_trend_pixel(rep(NA_real_, 24), t24, 6), + c(NA_real_, NA_real_)) + expect_equal(drift:::.dft_trend_pixel(c(0.5, 0.6, 0.55), t24[1:3], 6), + c(NA_real_, NA_real_)) +}) + +# ---- Theil-Sen slope + Mann-Kendall p direction ---------------------------- + +test_that(".dft_trend_pixel recovers a rising trend (positive, significant)", { + y <- 0.40 + 0.01 * seq_len(24) + out <- drift:::.dft_trend_pixel(y, t24, 6) + expect_gt(out[1], 0) # rising + expect_lt(out[2], 0.05) # significant +}) + +test_that(".dft_trend_pixel recovers a declining trend (negative, significant)", { + y <- 0.60 - 0.01 * seq_len(24) + out <- drift:::.dft_trend_pixel(y, t24, 6) + expect_lt(out[1], 0) + expect_lt(out[2], 0.05) +}) + +test_that(".dft_trend_pixel reports a flat series as ~0 slope, not significant", { + y <- 0.5 + rep(c(0.01, -0.01), 12) # wiggle, no trend + out <- drift:::.dft_trend_pixel(y, t24, 6) + expect_lt(abs(out[1]), 0.02) + expect_gt(out[2], 0.05) +}) + +test_that("Theil-Sen slope is robust to one anomalous season", { + y <- 0.40 + 0.01 * seq_len(24) # rising + y[19:24] <- y[19:24] - 0.30 # last 'year' crashes (a smoke/drought season) + out <- drift:::.dft_trend_pixel(y, t24, 6) + expect_gt(out[1], 0) # median pairwise slope stays positive +}) + +# ---- integration over a stack (terra, no network) -------------------------- + +test_that("dft_rast_trend reduces a synthetic stack to a 2-band raster", { + stk <- synthetic_break_stack() # 3x3, 72 monthly layers, one all-NA pixel + trend <- dft_rast_trend(stk, cores = 2) + + expect_s4_class(trend, "SpatRaster") + expect_equal(terra::nlyr(trend), 2) + expect_equal(names(trend), c("trend", "trend_p")) + # the injected step-drop pixel (pixel 1) trends down; the all-NA pixel is NA + vals <- terra::values(trend[["trend"]])[, 1] + expect_true(is.finite(vals[1])) + expect_true(is.na(vals[9])) +}) diff --git a/tests/testthat/test-dft_stac_config.R b/tests/testthat/test-dft_stac_config.R index 87e335f..9904a24 100644 --- a/tests/testthat/test-dft_stac_config.R +++ b/tests/testthat/test-dft_stac_config.R @@ -20,3 +20,36 @@ test_that("dft_stac_config returns valid config for esa-worldcover", { test_that("dft_stac_config errors on unknown source", { expect_error(dft_stac_config("bogus")) }) + +test_that("categorical sources carry no cube marker (absence = FALSE)", { + # backward-compat: the flat 4-field shape is pinned by expect_named above; + # consumers branch on isTRUE(cfg$cube), which must be FALSE here + expect_false(isTRUE(dft_stac_config("io-lulc")$cube)) + expect_false(isTRUE(dft_stac_config("esa-worldcover")$cube)) +}) + +test_that("dft_stac_config returns role-based cube config for sentinel-2-l2a", { + cfg <- dft_stac_config("sentinel-2-l2a") + expect_type(cfg, "list") + expect_true(isTRUE(cfg$cube)) + expect_equal(cfg$collection, "sentinel-2-l2a") + expect_match(cfg$stac_url, "planetarycomputer") + + # role -> asset map uses Planetary Computer band names (B04/B08/B11/SCL) + expect_named(cfg$roles, c("red", "nir", "swir16", "mask")) + expect_equal(cfg$roles$red, "B04") + expect_equal(cfg$roles$nir, "B08") + expect_equal(cfg$roles$swir16, "B11") + expect_equal(cfg$roles$mask, "SCL") + + # mask values (incl. 11 snow/ice), reflectance affine transform, temporal extent + expect_equal(cfg$mask_values, c(3L, 8L, 9L, 10L, 11L)) + expect_equal(cfg$scale, 1e-4) + expect_equal(cfg$offset, -0.1) + expect_match(cfg$available_datetime, "^\\d{4}-\\d{2}-\\d{2}/\\d{4}-\\d{2}-\\d{2}$") + + # baseline offset boundary: +1000 DN only from 2022-01-25, so pre-boundary + # scenes carry offset 0 and the fetch splits at the boundary + expect_equal(cfg$offset_boundary, "2022-01-25") + expect_equal(cfg$offset_before, 0) +}) diff --git a/tests/testthat/test-dft_stac_cube.R b/tests/testthat/test-dft_stac_cube.R new file mode 100644 index 0000000..2ec39b3 --- /dev/null +++ b/tests/testthat/test-dft_stac_cube.R @@ -0,0 +1,109 @@ +test_that("dft_stac_cube requires gdalcubes", { + skip_if(requireNamespace("gdalcubes", quietly = TRUE), + "gdalcubes is installed, can't test missing-package path") + aoi <- sf::st_read( + system.file("extdata", "example_aoi.gpkg", package = "drift"), + quiet = TRUE + ) + expect_error(dft_stac_cube(aoi), "gdalcubes") +}) + +test_that("dft_stac_cube rejects categorical sources", { + skip_if_not_installed("gdalcubes") + aoi <- sf::st_read( + system.file("extdata", "example_aoi.gpkg", package = "drift"), + quiet = TRUE + ) + expect_error(dft_stac_cube(aoi, source = "io-lulc"), "not a cube source") +}) + +# helpers for stac_cube_cache_key tests: a unit-square polygon (optionally +# shifted) and a key call with fixed defaults so each test varies one input +square_aoi <- function(dx = 0) { + sf::st_sfc( + sf::st_polygon(list(rbind( + c(0 + dx, 0), c(1 + dx, 0), c(1 + dx, 1), c(0 + dx, 1), c(0 + dx, 0) + ))), + crs = 32609 + ) +} + +cube_key <- function(aoi = square_aoi(), res = 10, target_crs = "EPSG:32609", + dt = "P1M", aggregation = "median", resampling = "bilinear", + stac_url = "https://example.com/stac", + collection = "sentinel-2-l2a", + band_assets = c("B08", "B04"), + datetime = "2019-01-01/2023-12-31", index = "kndvi", + cloud_cover_max = 60, mask_values = c(3, 8, 9, 10, 11), + scale = 1e-4, offset = -0.1, months = NULL, + offset_before = 0) { + drift:::stac_cube_cache_key( + aoi, res, target_crs, dt, aggregation, resampling, stac_url, collection, + band_assets, datetime, index, cloud_cover_max, mask_values, scale, offset, + months, offset_before + ) +} + +test_that("stac_cube_cache_key is deterministic and 12-char hex", { + expect_equal(cube_key(), cube_key()) + expect_match(cube_key(), "^[0-9a-f]{12}$") +}) + +test_that("stac_cube_cache_key changes with each cube-affecting parameter", { + base <- cube_key() + expect_false(cube_key(aoi = square_aoi(dx = 0.5)) == base) + expect_false(cube_key(res = 20) == base) + expect_false(cube_key(target_crs = "EPSG:32610") == base) + expect_false(cube_key(dt = "P1Y") == base) + expect_false(cube_key(aggregation = "mean") == base) + expect_false(cube_key(resampling = "near") == base) + expect_false(cube_key(collection = "landsat-c2-l2") == base) + expect_false(cube_key(band_assets = c("B08", "B11")) == base) + expect_false(cube_key(datetime = "2020-01-01/2020-12-31") == base) + expect_false(cube_key(index = "ndvi") == base) + expect_false(cube_key(cloud_cover_max = 20) == base) + expect_false(cube_key(mask_values = c(8, 9)) == base) + expect_false(cube_key(scale = 2.75e-5) == base) + expect_false(cube_key(offset = -0.2) == base) + expect_false(cube_key(months = 6:9) == base) + expect_false(cube_key(offset_before = -0.1) == base) +}) + +test_that("stac_cube_cache_key normalizes months order", { + expect_equal(cube_key(months = c(6, 7, 8, 9)), cube_key(months = c(9, 8, 7, 6))) +}) + +test_that("stac_cube_cache_key normalizes mask_values order and res type", { + expect_equal(cube_key(mask_values = c(3, 8, 9, 10)), + cube_key(mask_values = c(10, 9, 8, 3))) + expect_equal(cube_key(res = 10L), cube_key(res = 10)) +}) + +test_that("stac_cube_cache_key ignores sf attribute columns", { + bare <- square_aoi() + with_attrs <- sf::st_sf(name = "a", area = 1.5, geometry = bare) + expect_equal(cube_key(with_attrs), cube_key(bare)) +}) + +# Network end-to-end against the Planetary Computer. Opt-in only (env var), so +# the default `devtools::test()` stays network-free per the repo convention. +test_that("dft_stac_cube fetches an index stack end-to-end", { + skip_if(Sys.getenv("DRIFT_TEST_NETWORK") != "true", + "network test — set DRIFT_TEST_NETWORK=true to run") + skip_if_not_installed("gdalcubes") + aoi <- sf::st_read( + system.file("extdata", "example_aoi.gpkg", package = "drift"), + quiet = TRUE + ) + cache <- tempfile("drift_cube_") + dir.create(cache) + cube <- dft_stac_cube(aoi, index = "kndvi", + datetime = "2021-06-01/2021-08-31", dt = "P1M", + cache_dir = cache) + expect_s4_class(cube, "SpatRaster") + expect_equal(terra::nlyr(cube), 3) # 3 monthly layers + expect_false(anyNA(terra::time(cube))) # time set per layer + # second call hits the cache (one cube_.tif under the source dir) + expect_length(list.files(file.path(cache, "sentinel-2-l2a"), + pattern = "^cube_.*\\.tif$"), 1) +}) diff --git a/vignettes/trajectory-break-detection.Rmd b/vignettes/trajectory-break-detection.Rmd new file mode 100644 index 0000000..05773db --- /dev/null +++ b/vignettes/trajectory-break-detection.Rmd @@ -0,0 +1,176 @@ +--- +title: "Trajectories as a Check on Land-Cover Change" +output: bookdown::html_document2 +vignette: > + %\VignetteIndexEntry{Trajectories as a Check on Land-Cover Change} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r setup, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + fig.width = 7, + fig.height = 5 +) +``` + +This vignette applies continuous spectral-trajectory analysis to vegetation +change on the Neexdzii Kwa (Upper Bulkley) floodplain reach, over the same +2017–2023 window as the *Land Cover Change Detection* vignette, as a complement +to categorical land-cover classification. Vegetation change on this reach is +spatially patchy rather than reach-wide, so the analysis is framed at the patch +scale: which areas changed, whether a mapped land-cover transition is matched by +a measurable change in vegetation greenness (kNDVI), and when the change +occurred. + +The *Land Cover Change Detection* vignette answers part of that with categorical +labels: a stand mapped as Trees in 2017 and Rangeland in 2023 is a "Trees to +Rangeland" transition. This vignette adds a **continuous** kNDVI trajectory from +Sentinel-2 on the same reach and 2017–2023 window, and treats the two as +**complementary**. An abrupt **break** ([dft_rast_break()], via BFAST) dates +canopy removal; a gradual **trend** ([dft_rast_trend()]) tracks slow degradation +or recovery. Where they agree with the categorical map, a harvest is confirmed +and dated. Where the continuous signal fires *outside* the mapped polygons, it +flags removal the categorical product missed. + +## The pipeline + +```{r pipeline, eval = FALSE} +library(drift) + +aoi <- sf::st_read( + system.file("extdata", "example_aoi.gpkg", package = "drift") +) + +cube <- dft_stac_cube(aoi, index = "kndvi", + datetime = "2017-01-01/2023-12-31", months = 6:9) +breaks <- dft_rast_break(cube, start = c(2022, 1)) # abrupt: dated canopy removal +trend <- dft_rast_trend(cube) # gradual: which way, how fast +``` + +This vignette loads a small pre-computed artifact (generated by +`data-raw/vignette_data_break.R`) so every chunk runs live, no network: + +```{r load} +library(terra) + +d <- readRDS(system.file("testdata", "neexdzii_break.rds", package = "drift")) +break_mag <- terra::unwrap(d$break_r) +trend <- terra::unwrap(d$trend_r) +tree_trans <- d$tree_trans # LULC Trees -> Rangeland polygons +aoi <- d$aoi +``` + +## Where has canopy been removed? + +`dft_rast_break()` measures the size of an abrupt drop in kNDVI. The strong +losses are coherent blocks — the shape of harvest — not scattered noise. The +beige outlines are the stands the land-cover product maps as Trees to +Rangeland. + +```{r break-map, fig.cap = "Break magnitude (red = a large kNDVI drop). Beige outlines are LULC Trees-to-Rangeland. Some red blocks sit inside the outlines (confirmed, dated harvest); some fall outside them (removal the categorical map missed); some outlines carry no break (a label change without real canopy loss).", fig.height = 5.5} +terra::plot(break_mag, col = hcl.colors(21, "Blue-Red 3", rev = TRUE), + range = c(-0.45, 0.45), main = "kNDVI break magnitude") +plot(sf::st_geometry(sf::st_transform(tree_trans, terra::crs(break_mag))), + add = TRUE, border = "#7d7a52", lwd = 1.4) +``` + +The three cases in that figure are the whole point. Red *inside* an outline is a +harvest both methods agree on — and the break dates it. Red *outside* the +outlines is real canopy loss the categorical product did not flag. An outline +with no red is a mapped transition not backed by a spectral drop — likely a +borderline pixel changing label rather than trees coming off. The two products +catch different real things; used together they are more complete than either +alone. + +## Gradual change, and recovery + +Not all change is a clean step. `dft_rast_trend()` gives a robust (Theil-Sen) +slope per pixel — the average kNDVI change **per year over 2017–2023** — for the +slow signals: a stand thinning, or a bar revegetating. Because a linear slope +spreads an abrupt cut across the whole record, it *understates* the blocks (use +the break magnitude above for those); its value is the gradual story the +categorical labels cannot show. + +```{r trend-map, fig.cap = "kNDVI trend per year, 2017-2023 (red = declining, blue = recovering). Much of the open floodplain is greening; decline concentrates on the harvested blocks and some channel edges.", fig.height = 5.5} +terra::plot(trend[["trend"]], col = hcl.colors(21, "Blue-Red 3", rev = TRUE), + range = c(-0.04, 0.04), main = "kNDVI trend per year (2017-2023)") +``` + +## Explore it + +The interactive map below carries the same basemaps and framing as the +land-cover vignette. Toggle the break magnitude, the trend, and the LULC +Trees-to-Rangeland polygons over Esri or Google satellite imagery, and zoom to +the loss blocks to see the harvest on the ground alongside what each method +detected. + +```{r interactive, fig.height = 6} +library(leaflet) + +aoi4326 <- sf::st_transform(aoi, 4326) +bb <- sf::st_bbox(aoi4326) +break_ll <- terra::project(break_mag, "EPSG:4326") +trend_ll <- terra::project(trend[["trend"]], "EPSG:4326") +pal_bk <- colorNumeric(hcl.colors(11, "Blue-Red 3", rev = TRUE), + domain = c(-0.45, 0.45), na.color = "transparent") +pal_tr <- colorNumeric(hcl.colors(11, "Blue-Red 3", rev = TRUE), + domain = c(-0.04, 0.04), na.color = "transparent") +basemaps <- c("Light" = "CartoDB.Positron", + "Esri Satellite" = "Esri.WorldImagery", + "Google Satellite" = "https://mt1.google.com/vt/lyrs=s&x={x}&y={y}&z={z}", + "OpenTopoMap" = "OpenTopoMap") + +m <- leaflet::leaflet() |> + leaflet::setView(mean(bb[c("xmin", "xmax")]), mean(bb[c("ymin", "ymax")]), + zoom = 14) +for (i in seq_along(basemaps)) { + bm <- basemaps[[i]] + nm <- names(basemaps)[i] + m <- if (grepl("^https?://", bm)) { + leaflet::addTiles(m, urlTemplate = bm, group = nm) + } else { + leaflet::addProviderTiles(m, bm, group = nm) + } +} +overlays <- c("Break magnitude", "Trend (per yr)", "LULC Trees -> Rangeland", "AOI") +m <- m |> + leaflet::addRasterImage(break_ll, colors = pal_bk, opacity = 0.8, + group = "Break magnitude", project = FALSE) |> + leaflet::addRasterImage(trend_ll, colors = pal_tr, opacity = 0.8, + group = "Trend (per yr)", project = FALSE) |> + leaflet::addPolygons(data = sf::st_transform(tree_trans, 4326), + fillColor = "#e3e2c3", fillOpacity = 0.6, + color = "#b3b083", weight = 1, + group = "LULC Trees -> Rangeland") |> + leaflet::addPolygons(data = aoi4326, fill = FALSE, color = "red", weight = 2, + group = "AOI") |> + leaflet::addLayersControl( + baseGroups = names(basemaps), overlayGroups = overlays, + options = leaflet::layersControlOptions(collapsed = FALSE) + ) |> + leaflet::hideGroup("Trend (per yr)") +if (requireNamespace("leaflet.extras", quietly = TRUE)) { + m <- leaflet.extras::addFullscreenControl(m) +} +m +``` + +## When to use which + +The tools are complements, patch by patch: + +- **Categorical** (`dft_stac_fetch()` + `dft_rast_transition()`) — for named + transitions and hectares. Treat its transitions as hypotheses to confirm. +- **Break** (`dft_rast_break()`) — to confirm and **date** canopy removal, and to + catch harvest the categorical map misses. +- **Trend** (`dft_rast_trend()`) — for gradual degradation and recovery the annual + labels cannot show. + +One caveat for this landscape: in peak summer, deciduous riparian trees and grass +are both green, so kNDVI cannot cleanly separate them — the continuous tools are +strongest on clear canopy removal and gradual trends, and weakest at telling +riparian tree from rangeland. Used with the categorical map, they say whether a +mapped change is *real*, which *way* it is going, and *when* it happened.