Fix large-floodplain OOM in LULC transition + vectorize (terra-native, changes_only) — v0.4.0

.gitignore

@@ -10,3 +10,6 @@ vignettes/*.html
 
 # local STAC cube cache (generated by data-raw/vignette_data_break.R)
 data-raw/.break_cache/
+
+# terra sidecars regenerated when tests/benchmarks read inst/extdata rasters
+inst/extdata/*.aux.xml

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: drift
 Title: Detecting Riparian and Inland Floodplain Transitions
-Version: 0.3.0
-Date: 2026-07-08
+Version: 0.4.0
+Date: 2026-07-09
 Authors@R: c(
     person("Allan", "Irvine", , "[email protected]", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0002-3495-2128")),

NEWS.md

@@ -1,3 +1,9 @@
+# drift 0.4.0
+
+- Categorical land-cover change detection no longer exhausts memory on large-floodplain AOIs (#34, #28). `dft_rast_transition()` was rewritten to stream entirely through `terra` — transitions are encoded and filtered with raster arithmetic, `terra::subst()`, `patches()`, and a single `terra::freq()`, with no `terra::values()` pull and no full-grid R vectors — so peak memory scales with the number of distinct transitions and patches, not the grid size (producer-only peak at 16M cells dropped from 2.66 GB to 1.63 GB). Output is byte-identical to the previous version, verified by a golden snapshot across the full parameter matrix.
+- `dft_transition_vectors()` gains `changes_only` (default `FALSE`): when `TRUE`, stable (`from == to`) transitions are dropped at the raster level before polygonizing, so `terra::as.polygons()` only builds geometry for actual change patches. On a fragmented floodplain — where the stable mosaic is most of the grid and polygonization dominates memory — this roughly halves peak use (a 9M-cell, 415k-patch benchmark went from 3.83 GB to 1.71 GB). The result equals the default output filtered to change patches. When `patch_area_min` is set, small patches are also dropped before polygonizing, with identical output.
+- `patch_id` in `dft_transition_vectors()` is numbered over the surviving patches when filtering drops any, and an empty result now carries the zone column so per-zone results bind cleanly.
+
 # 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.

R/dft_rast_transition.R

@@ -30,6 +30,12 @@
 #'   - `removed`: A `SpatRaster` of transitions removed by `patch_area_min`
 #'     filtering, or `NULL` when no filtering is applied. Same factor
 #'     encoding as `raster`.
+#'
+#' @details
+#' The transition raster, filters, and patch removal are computed with streamed
+#' `terra` operations (`ifel`, `patches`, `freq`) on the underlying class codes —
+#' no full-grid vectors are pulled into R — so peak memory scales with the number
+#' of distinct transitions and patches, not the grid size.
 #' @export
 #' @examples
 #' r17 <- terra::rast(system.file("extdata", "example_2017.tif", package = "drift"))
@@ -83,65 +89,55 @@ dft_rast_transition <- function(x,
   dft_check_crs(r_from, "dft_rast_transition")
   dft_check_crs(r_to, "dft_rast_transition")
 
-  # Resolve class names to codes
+  # code -> class name lookup (small; used only for factor labels + summary)
   code_lookup <- stats::setNames(class_table$class_name, class_table$code)
 
-  # Get raw integer values (strip factor)
-  v_from <- terra::values(r_from)[, 1]
-  v_to <- terra::values(r_to)[, 1]
-
-  # Map codes to class names
-  name_from <- code_lookup[as.character(v_from)]
-  name_to <- code_lookup[as.character(v_to)]
-
-  # Encode transitions: from_code * 1000 + to_code (supports up to 999 classes)
-  trans_code <- v_from * 1000L + v_to
-
-  # Build mask for filters
-
-  keep <- rep(TRUE, length(trans_code))
-  if (!is.null(from_class)) keep <- keep & (name_from %in% from_class)
-  if (!is.null(to_class)) keep <- keep & (name_to %in% to_class)
-
-  # Also mask where either raster is NA
-  keep <- keep & !is.na(v_from) & !is.na(v_to)
+  # Strip factor to raw integer codes as streamed rasters. `* 1L` returns a new
+  # non-factor raster of codes, preserves NA, and does NOT mutate the inputs
+  # (they are references into the caller's list). Encode each transition as
+  # from_code * 1000 + to_code; NA propagates wherever either input is NA. No
+  # full-grid R vector is materialized.
+  code_from <- r_from * 1L
+  code_to <- r_to * 1L
+  r_trans <- code_from * 1000L + code_to
 
-  trans_code[!keep] <- NA_integer_
+  # from_class / to_class filters as integer code-set membership (streamed)
+  r_trans <- apply_codeset(r_trans, code_from, from_class, class_table)
+  r_trans <- apply_codeset(r_trans, code_to, to_class, class_table)
 
+  cell_area_m2 <- prod(terra::res(r_from))
 
-  # Build transition raster
-  r_trans <- terra::rast(r_from)
-  terra::values(r_trans) <- trans_code
-
-  # Filter small patches of changed pixels
-  removed_codes <- NULL
+  # Filter small patches of *changed* pixels (streamed; no rep()/values())
+  r_removed <- NULL
   if (!is.null(patch_area_min) && patch_area_min > 0) {
-    # Binary raster of actual changes only (from != to), not same-class pixels
-    changed <- !is.na(trans_code) & (v_from != v_to)
-    r_changed <- terra::rast(r_from)
-    r_changed_vals <- rep(NA_integer_, length(trans_code))
-    r_changed_vals[changed] <- 1L
-    terra::values(r_changed) <- r_changed_vals
+    r_changed <- terra::ifel(!is.na(r_trans) & (code_from != code_to), 1L, NA)
     p <- terra::patches(r_changed, directions = 8)
-    cell_area_m2 <- prod(terra::res(r_from))
-    f <- terra::freq(p)
-    f <- f[!is.na(f$value), ]
-    small_ids <- f$value[f$count * cell_area_m2 < patch_area_min]
-    if (length(small_ids) > 0) {
-      p_vals <- terra::values(p)[, 1]
-      mask_small <- !is.na(p_vals) & (p_vals %in% small_ids)
-      removed_codes <- trans_code[mask_small]
-      trans_code[mask_small] <- NA_integer_
-      terra::values(r_trans) <- trans_code
+    f <- tryCatch(terra::freq(p), error = function(e) NULL)   # NULL when no changes
+    if (!is.null(f)) {
+      f <- f[!is.na(f$value), , drop = FALSE]
+      small_ids <- f$value[f$count * cell_area_m2 < patch_area_min]
+      if (length(small_ids) > 0) {
+        # 1 at small-patch cells, NA elsewhere (incl. p == NA / stable cells).
+        # subst() is exact-match and scales; SpatRaster `%in%` is not dispatched
+        # when terra is imported (not attached).
+        sm <- terra::subst(p, small_ids, 1L, others = NA)
+        r_removed <- terra::ifel(!is.na(sm), r_trans, NA)   # capture removed codes first
+        r_trans <- terra::ifel(!is.na(sm), NA, r_trans)     # then drop them
+      }
     }
   }
 
-  # Build factor table from observed transitions
-  valid <- !is.na(trans_code)
-  unique_codes <- sort(unique(trans_code[valid]))
+  # Observed transition codes + counts from a single native freq (value == code,
+  # NA excluded). freq() errors on an all-NA raster -> treat as no transitions.
+  freq_tbl <- tryCatch(terra::freq(r_trans), error = function(e) NULL)
+  if (!is.null(freq_tbl)) {
+    freq_tbl <- freq_tbl[!is.na(freq_tbl$value), , drop = FALSE]
+  }
+
+  m2_to_unit <- switch(unit, "m2" = 1, "ha" = 1e-4, "km2" = 1e-6)
+  cell_area <- cell_area_m2 * m2_to_unit
 
-  if (length(unique_codes) == 0) {
-    # No transitions found — return empty
+  if (is.null(freq_tbl) || nrow(freq_tbl) == 0) {
     terra::set.cats(r_trans, layer = 1,
                     value = data.frame(id = integer(0), transition = character(0)))
     summary_tbl <- tibble::tibble(
@@ -151,58 +147,57 @@ dft_rast_transition <- function(x,
     return(list(raster = r_trans, summary = summary_tbl, removed = NULL))
   }
 
-  from_codes <- unique_codes %/% 1000L
-  to_codes <- unique_codes %% 1000L
+  freq_tbl <- freq_tbl[order(freq_tbl$value), , drop = FALSE]
+  codes <- freq_tbl$value
+  from_codes <- codes %/% 1000L
+  to_codes <- codes %% 1000L
   labels <- paste0(
     code_lookup[as.character(from_codes)], " -> ",
     code_lookup[as.character(to_codes)]
   )
-
-  lvl_df <- data.frame(id = unique_codes, transition = labels)
-  terra::set.cats(r_trans, layer = 1, value = lvl_df)
+  terra::set.cats(r_trans, layer = 1,
+                  value = data.frame(id = codes, transition = labels))
 
   # Summary table
-  m2_to_unit <- switch(unit, "m2" = 1, "ha" = 1e-4, "km2" = 1e-6)
-  res <- terra::res(r_from)
-  cell_area <- res[1] * res[2] * m2_to_unit
-
-  freq_tbl <- terra::freq(r_trans)
-  freq_tbl <- freq_tbl[!is.na(freq_tbl$value), ]
-  total_valid <- sum(!is.na(trans_code))
-
-  # freq on a factor raster returns labels in $value — map back to codes
-  label_to_code <- stats::setNames(lvl_df$id, lvl_df$transition)
-  freq_codes <- label_to_code[as.character(freq_tbl$value)]
-
+  total_valid <- sum(freq_tbl$count)
   summary_tbl <- tibble::tibble(
-    from_class = code_lookup[as.character(freq_codes %/% 1000L)],
-    to_class = code_lookup[as.character(freq_codes %% 1000L)],
+    from_class = code_lookup[as.character(from_codes)],
+    to_class = code_lookup[as.character(to_codes)],
     n_cells = as.integer(freq_tbl$count),
     area = freq_tbl$count * cell_area,
     pct = round(freq_tbl$count / total_valid * 100, 2)
   )
-
   summary_tbl <- summary_tbl[order(summary_tbl$n_cells, decreasing = TRUE), ]
 
-  # Build removed raster when patch filtering was applied
-  r_removed <- NULL
-  if (!is.null(removed_codes)) {
-    r_removed <- terra::rast(r_from)
-    rem_vals <- rep(NA_integer_, terra::ncell(r_from))
-    rem_positions <- which(!is.na(p_vals) & (p_vals %in% small_ids))
-    rem_vals[rem_positions] <- removed_codes
-    terra::values(r_removed) <- rem_vals
-    # Build factor table from removed codes
-    rem_unique <- sort(unique(removed_codes[!is.na(removed_codes)]))
-    rem_from <- rem_unique %/% 1000L
-    rem_to <- rem_unique %% 1000L
-    rem_labels <- paste0(
-      code_lookup[as.character(rem_from)], " -> ",
-      code_lookup[as.character(rem_to)]
-    )
-    terra::set.cats(r_removed, layer = 1,
-                    value = data.frame(id = rem_unique, transition = rem_labels))
+  # Factor levels for the removed raster (only when patch filtering removed cells)
+  if (!is.null(r_removed)) {
+    rf <- terra::freq(r_removed)
+    rf <- rf[!is.na(rf$value), , drop = FALSE]
+    rf <- rf[order(rf$value), , drop = FALSE]
+    terra::set.cats(r_removed, layer = 1, value = data.frame(
+      id = rf$value,
+      transition = paste0(code_lookup[as.character(rf$value %/% 1000L)], " -> ",
+                          code_lookup[as.character(rf$value %% 1000L)])
+    ))
   }
 
   list(raster = r_trans, summary = summary_tbl, removed = r_removed)
 }
+
+#' Mask a transition raster to a from/to class set (streamed)
+#'
+#' Reproduces the old full-grid `name %in% from_class` filter using integer
+#' code-set membership on a code raster, so no full-grid character vector is
+#' built. `code_r %in% keep_codes` is FALSE (not NA) at NA cells, matching the
+#' base-R string filter's NA handling. `NULL` selection is a no-op; an empty
+#' selection (no class name matched) masks everything to NA.
+#' @noRd
+apply_codeset <- function(r_trans, code_r, names_sel, class_table) {
+  if (is.null(names_sel)) return(r_trans)
+  keep_codes <- class_table$code[class_table$class_name %in% names_sel]
+  if (length(keep_codes) == 0L) return(r_trans * NA)   # impossible filter -> all NA
+  # subst(): 1 at in-set codes, NA at out-of-set and NA cells. Avoids SpatRaster
+  # `%in%`, which is not dispatched when terra is imported (not attached).
+  keep_mask <- terra::subst(code_r, keep_codes, 1L, others = NA)
+  terra::ifel(!is.na(keep_mask), r_trans, NA)
+}

R/dft_transition_vectors.R

@@ -17,14 +17,22 @@
 #'   Required when `zones` is supplied.
 #' @param patch_area_min Numeric or `NULL`. Minimum patch area in m². Patches
 #'   smaller than this are dropped before returning. `NULL` (default) keeps all.
+#' @param changes_only Logical. When `TRUE`, drop stable (`from == to`)
+#'   transitions before polygonizing, so only actual change patches are
+#'   vectorized. On a floodplain the stable mosaic is most of the grid, so this
+#'   is the main memory lever for large AOIs. Default `FALSE` keeps every patch
+#'   (including stable ones).
 #'
 #' @return An `sf` data frame (polygon geometry) with columns:
 #'   - `patch_id` (integer) — connected component ID, numbered in raster
-#'     scan order
+#'     scan order over the returned patches
 #'   - `transition` (character) — transition label (e.g. "Trees -> Rangeland")
 #'   - `area_ha` (numeric) — patch area in hectares
 #'   - Zone column (if `zones` supplied) — from spatial intersection
 #'
+#'   When `patch_area_min` or `changes_only` drop patches, `patch_id` is numbered
+#'   over the surviving patches (dense `1..n`), not the pre-filter grid.
+#'
 #' @export
 #' @examples
 #' r17 <- terra::rast(system.file("extdata", "example_2017.tif", package = "drift"))
@@ -39,10 +47,15 @@
 #' # Filter to large patches only
 #' patches_large <- dft_transition_vectors(result$raster, patch_area_min = 1000)
 #' head(patches_large)
+#'
+#' # Only actual changes (drop stable from == to patches before polygonizing)
+#' changes <- dft_transition_vectors(result$raster, changes_only = TRUE)
+#' head(changes)
 dft_transition_vectors <- function(x,
                                    zones = NULL,
                                    zone_col = NULL,
-                                   patch_area_min = NULL) {
+                                   patch_area_min = NULL,
+                                   changes_only = FALSE) {
   if (!inherits(x, "SpatRaster")) {
     stop("`x` must be a SpatRaster (e.g. from dft_rast_transition()$raster).",
          call. = FALSE)
@@ -70,18 +83,56 @@ dft_transition_vectors <- function(x,
     }
   }
 
+  if (!is.logical(changes_only) || length(changes_only) != 1 ||
+        is.na(changes_only)) {
+    stop("`changes_only` must be a single logical (TRUE or FALSE).", call. = FALSE)
+  }
+
+  # Optionally drop stable (from == to) transitions before polygonizing, so
+  # as.polygons() only builds geometry for actual change patches. Stable codes
+  # are those where from_code == to_code (id %/% 1000 == id %% 1000). On a
+  # floodplain the stable mosaic is most of the grid, so this caps the working
+  # set to what callers actually keep.
+  if (isTRUE(changes_only)) {
+    ct <- terra::cats(x)[[1]]
+    stable_ids <- ct$id[(ct$id %/% 1000L) == (ct$id %% 1000L)]
+    if (length(stable_ids) > 0) {
+      # NA at stable cells, 1 at change cells; mask keeps x's factor levels
+      change_mask <- terra::subst(x * 1L, stable_ids, NA, others = 1L)
+      x <- terra::mask(x, change_mask)
+    }
+  }
+
   # Single pass: 8-connected components of same-valued cells.
   # Requires terra >= 1.8-10 (earlier versions falsely connect patches
   # across the left/right raster edges with values = TRUE).
   p <- terra::patches(x, directions = 8, values = TRUE)
   names(p) <- "pid"
+
+  # When filtering by area, drop small patches at the raster level BEFORE
+  # polygonizing (as.polygons is the patch-count-driven hotspot). For axis-
+  # aligned raster polygons st_area == count * cell_area, so this drops a strict
+  # subset of what the trailing st_area filter removes -> identical output.
+  if (!is.null(patch_area_min) && patch_area_min > 0) {
+    fp <- tryCatch(terra::freq(p), error = function(e) NULL)
+    if (!is.null(fp)) {
+      cell_area_m2 <- prod(terra::res(p))
+      small <- fp$value[!is.na(fp$value) & fp$count * cell_area_m2 < patch_area_min]
+      if (length(small) > 0) p <- terra::subst(p, small, NA)
+    }
+  }
+
   polys_sf <- sf::st_as_sf(terra::as.polygons(p))
 
   if (nrow(polys_sf) == 0) {
-    return(sf::st_sf(
+    out <- sf::st_sf(
       patch_id = integer(0), transition = character(0),
       area_ha = numeric(0), geometry = sf::st_sfc(crs = sf::st_crs(x))
-    ))
+    )
+    # match the zone-attributed schema so per-zone results bind cleanly when a
+    # zone has no (or, under changes_only, no change) patches
+    if (!is.null(zones)) out[[zone_col]] <- zones[[zone_col]][0]
+    return(out)
   }
 
   # Map each patch to its transition label, touching only the non-NA cells