Single-pass patch vectorization (fixes OOM on large-extent rasters)

DESCRIPTION

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

NEWS.md

@@ -1,3 +1,8 @@
+# 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.
+- terra dependency floored at `>= 1.8-10`: earlier versions had an edge-wraparound bug in `patches(values = TRUE)` that silently merged patches touching opposite raster edges.
+
 # drift 0.2.3
 
 - Fix silent cross-AOI cache collision in `dft_stac_fetch()` (#25). Cache files were keyed by source + year only, so fetching a second AOI with the same source/year silently returned the first AOI's raster masked to the second AOI's extent. Cache filenames now include a hash of the AOI geometry and all fetch-affecting parameters (`res`, `crs`, `dt`, `aggregation`, `resampling`, `stac_url`, `collection`, `asset`). Existing caches re-fetch on first use after upgrading; `dft_cache_clear()` reclaims the orphaned old-format files.

R/dft_transition_vectors.R

@@ -5,6 +5,10 @@
 #' transition type. Useful for QA in GIS, spatial attribution to management
 #' zones, and patch-level reporting.
 #'
+#' Patches are 8-connected components of same-valued cells, computed in a
+#' single pass over the grid, so large sparse rasters vectorize without
+#' per-class memory cost.
+#'
 #' @param x A factor `SpatRaster` from [dft_rast_transition()] (the `$raster`
 #'   element). Must have a projected CRS.
 #' @param zones Optional `sf` polygon layer for spatial attribution. Any
@@ -15,7 +19,8 @@
 #'   smaller than this are dropped before returning. `NULL` (default) keeps all.
 #'
 #' @return An `sf` data frame (polygon geometry) with columns:
-#'   - `patch_id` (integer) — connected component ID
+#'   - `patch_id` (integer) — connected component ID, numbered in raster
+#'     scan order
 #'   - `transition` (character) — transition label (e.g. "Trees -> Rangeland")
 #'   - `area_ha` (numeric) — patch area in hectares
 #'   - Zone column (if `zones` supplied) — from spatial intersection
@@ -65,62 +70,31 @@ dft_transition_vectors <- function(x,
     }
   }
 
-  lvls <- terra::cats(x)[[1]]
-  vals <- terra::values(x)[, 1]
-
-  if (all(is.na(vals))) {
-    return(sf::st_sf(
-      patch_id = integer(0), transition = character(0),
-      area_ha = numeric(0), geometry = sf::st_sfc(crs = sf::st_crs(x))
-    ))
-  }
-
-  # Build unique patch IDs per transition type
-  patch_ids <- rep(NA_integer_, terra::ncell(x))
-  patch_transition <- integer(0)
-  offset <- 0L
-
-  for (i in seq_len(nrow(lvls))) {
-    code <- lvls$id[i]
-    mask <- which(vals == code)
-    if (length(mask) == 0) next
+  # 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"
+  polys_sf <- sf::st_as_sf(terra::as.polygons(p))
 
-    r_mask <- terra::rast(x)
-    mask_vals <- rep(NA_integer_, terra::ncell(x))
-    mask_vals[mask] <- 1L
-    terra::values(r_mask) <- mask_vals
-    p <- terra::patches(r_mask, directions = 8)
-    p_vals <- terra::values(p)[, 1]
-
-    valid <- !is.na(p_vals)
-    unique_pids <- sort(unique(p_vals[valid]))
-    for (pid in unique_pids) {
-      offset <- offset + 1L
-      patch_ids[valid & p_vals == pid] <- offset
-      patch_transition <- c(patch_transition, code)
-    }
-  }
-
-  if (offset == 0L) {
+  if (nrow(polys_sf) == 0) {
     return(sf::st_sf(
       patch_id = integer(0), transition = character(0),
       area_ha = numeric(0), geometry = sf::st_sfc(crs = sf::st_crs(x))
     ))
   }
 
-  # Vectorize patch IDs
-  r_pid <- terra::rast(x)
-  names(r_pid) <- "pid"
-  terra::values(r_pid) <- patch_ids
-  polys <- terra::as.polygons(r_pid)
-  polys_sf <- sf::st_as_sf(polys)
-
-  # Map patch IDs to transition labels
-  code_to_label <- stats::setNames(lvls$transition, lvls$id)
-  polys_sf$patch_id <- as.integer(polys_sf$pid)
-  polys_sf$transition <- as.character(
-    code_to_label[as.character(patch_transition[polys_sf$pid])]
-  )
+  # Map each patch to its transition label, touching only the non-NA cells
+  cell_idx <- terra::cells(p)
+  pid_at <- terra::extract(p, cell_idx)[, 1]
+  lab_at <- as.character(terra::extract(x, cell_idx)[, 1])
+  first <- !duplicated(pid_at)
+  lab_map <- stats::setNames(lab_at[first], pid_at[first])
+
+  polys_sf$transition <- unname(lab_map[as.character(polys_sf$pid)])
+  # Cell values with no cats() entry have no label — drop their patches
+  polys_sf <- polys_sf[!is.na(polys_sf$transition), ]
+  polys_sf$patch_id <- seq_len(nrow(polys_sf))
   polys_sf$area_ha <- as.numeric(sf::st_area(polys_sf)) * 1e-4
   polys_sf$pid <- NULL
 

planning/archive/2026-07-issue-27-transition-vectors-oom/README.md

@@ -0,0 +1,21 @@
+# Issue #27 — dft_transition_vectors OOMs on large-extent rasters
+
+## Outcome
+
+Fixed the OOM by replacing the per-class/per-patch loop in `dft_transition_vectors()` with a
+single `terra::patches(x, directions = 8, values = TRUE)` pass (8-connected components of
+same-valued cells) plus a sparse pid→label map via `terra::cells()`/`terra::extract()`. The old
+code scaled as ncell × n_patches — the sleeper was the per-patch remap allocating two full-grid
+logicals per patch — extrapolating to ~18+ GB on the 102.6M-cell UFRA floodplain. Approach chosen
+by empirical benchmark during planning: single-pass was behavior-identical (pinned by a
+regression test captured from the old implementation: 185 patches / 123.11 ha / 57 at
+`patch_area_min = 1000`) and 55× faster (1.9 s vs 122 s at 24M cells); the issue's tiling and
+sparse-window options were behavior-changing or degenerate on mainstem-shaped data, and
+polygonize+explode was disqualified (GDAL is 4-connected). Key learnings: terra ≥ 1.8-10 is a
+correctness floor (`values = TRUE` edge-wraparound bug); `patch_id` ordering became scan-order
+(docs only promised "connected component ID"); the producer `dft_rast_transition()` has the same
+disease (6+ full-grid vectors incl. two character) — filed as #28 rather than scope-creeping.
+Released as v0.2.4.
+
+Closed by: commits 157c66c / baa45c3 / 69f2578, PR pending (branch
+`27-dft-transition-vectors-ooms-on-large-ext`)

planning/archive/2026-07-issue-27-transition-vectors-oom/findings.md

@@ -0,0 +1,93 @@
+# Findings — dft_transition_vectors OOMs on large-extent rasters (#27)
+
+## Issue context
+
+### Summary
+
+`dft_transition_vectors()` processes the **entire raster grid** rather than the ~2% of
+cells that hold data, and does so **once per transition class**. On a large-extent
+floodplain this exhausts memory and the R process is OOM-killed.
+
+### Where
+
+`R/dft_transition_vectors.R`. For each transition class in `terra::cats(x)` it:
+- allocates `rep(NA_integer_, terra::ncell(x))` (a full-grid vector), and
+- runs `terra::patches(r_mask, directions = 8)` over the full grid,
+then finally `terra::as.polygons()` over a full-grid patch-id raster. Nothing is cropped
+to the data extent, and `trim()` does not help when the data follows a mainstem across
+the whole bounding box.
+
+### Repro / evidence
+
+Upper Fraser (UFRA) chinook `ch_ff04` floodplain via a per-area pipeline:
+- transition raster: **8637 x 11875 = 102.6M cells**, extent 119 x 86 km
+- non-NA (actual transitions): **1.87%** of the grid
+- distinct transition classes (loop iterations): **56**
+- `terra::trim(x)` returns **100%** of the grid (data spans the full bbox)
+
+Result: OOM-killed at `dft_transition_vectors()`. The smaller MORR floodplain
+(74.0M-cell grid, fewer classes) completes, so this is grid-cell-count driven, not
+floodplain-area driven.
+
+### Fix options (from issue)
+
+1. **Tile internally**: split into column/row tiles bounded to N cells, vectorize each,
+   `rbind`. (Stop-gap applied in the pipeline driver.)
+2. **Work on sparse cells**: derive patch ids from `which(!is.na(vals))` indices;
+   run `patches()` on a cropped/masked window per class.
+3. **Single `as.polygons()` per class-set** using a combined patch-id raster built from
+   sparse indices.
+
+## Plan-mode exploration (2026-07-07)
+
+### Failure anatomy (beyond the issue)
+
+- The sleeper is the per-patch remap loop (`R/dft_transition_vectors.R:95-101`):
+  `patch_ids[valid & p_vals == pid] <- offset` allocates TWO full-grid logical vectors
+  per patch. Thousands of patches → TBs of allocation churn. The current implementation
+  scales as ncell × n_patches, not just ncell × n_classes.
+- Benchmarked current impl: 122 s / 4.32 GB at 24M cells (48 classes, ~2% non-NA,
+  1,232 patches). Extrapolates to ~18+ GB at the real 102.6M-cell case → matches OOM.
+- `dft_rast_transition()` (producer) has the same disease independently: 6+ full-grid
+  R vectors including two full-grid character vectors (`name_from`/`name_to`,
+  R/dft_rast_transition.R:90-109, ~800 MB each at 102.6M cells) regardless of options,
+  plus a full-grid `patches()` path when `patch_area_min` is set (:118-137, :188-205).
+  Needs its own refactor + correctness harness → follow-up issue, not #27 scope.
+
+### Approach evaluation (empirical, terra 1.9.11, /usr/bin/time -l)
+
+| Approach | 24M cells | Output vs current |
+|---|---|---|
+| current | 122 s / 4.32 GB | reference |
+| **single-pass `patches(values = TRUE)` (chosen)** | **2.2 s / 1.94 GB** | **identical** |
+| sparse + per-class window (issue opt 2) | 54 s / 7.24 GB | identical |
+| `as.polygons(dissolve)` + `st_cast` explode | 2.1 s / 1.13 GB | WRONG (4-connected) |
+
+- `terra::patches(x, directions = 8, values = TRUE)` computes 8-connected components of
+  same-valued cells in one C++ pass — verified: respects class boundaries, merges
+  diagonal same-class cells (incl. crossed diagonals), treats 0 as a real class,
+  factor input fine, all-NA safe.
+- Verified identical to current output on real fixtures (326x314, 19 classes): 185
+  patches, identical per-class counts, sorted areas, total 123.11 ha, union geometry
+  `st_equals`; `patch_area_min = 1000` → 57 patches both ways.
+- Sparse+window degenerates on mainstem-shaped data: every class's bounding window is
+  nearly the full grid, so it does full-grid work × n_classes.
+- Polygonize+explode disqualified: GDAL polygonize is 4-connected (no 8-conn option in
+  `terra::as.polygons`), silently splits diagonal joins (275 vs 185 patches on fixtures)
+  and changes `patch_area_min` semantics.
+
+### Load-bearing details for implementation
+
+- **terra `(>= 1.8-10)` floor is a correctness requirement**: `values = TRUE` added in
+  1.8-5; edge-wraparound bug (patches falsely connected across left/right raster edges,
+  rspatial/terra#1675) fixed in 1.8-10. Older terra would silently corrupt patches.
+- pid→label mapping: `terra::cells(p)` + `terra::extract()` on the factor raster touches
+  only non-NA cells and returns labels straight from `cats()` — no full-grid pull, no
+  label drift.
+- Drop rows whose value has no `cats()` entry after mapping (old code silently skipped
+  them via the `lvls` loop).
+- `patch_id` numbering becomes scan-order rather than class-major — docs only promise
+  "connected component ID"; no test or documented behavior depends on ordering.
+- All-NA early return moves after `as.polygons()` (0 features → same empty sf).
+- Bit-identical output verified under `terraOptions(memmax = 0.3)`. At the real
+  102.6M-cell size expect ~7-8 GB peak (dominated by `as.polygons`) vs OOM before.

planning/archive/2026-07-issue-27-transition-vectors-oom/progress.md

@@ -0,0 +1,15 @@
+# Progress — dft_transition_vectors OOMs on large-extent rasters (#27)
+
+## Session 2026-07-07
+
+- Plan-mode exploration — approaches benchmarked empirically, phases approved by user
+- Created branch `27-dft-transition-vectors-ooms-on-large-ext` off main
+- Scaffolded PWF baseline from issue #27 with approved phases
+- Phase 1 complete: single-pass `patches(values = TRUE)` rewrite + terra floor + 3 new tests
+  (suite: 211 pass, 0 fail; lint clean). Regression guard captured from OLD implementation
+  before the rewrite: 185 patches / 123.11 ha / 57 at patch_area_min = 1000 — new code matches.
+- Final-implementation benchmark: 24M cells, 4,799 patches → 1.9 s (old code: 122 s with only
+  1,232 patches)
+- Phase 2 complete: roxygen single-pass/scan-order notes, NEWS 0.2.4, follow-up issue #28
+  filed for dft_rast_transition's full-grid pipeline, version bump
+- Next: /planning-archive + /gh-pr-push

planning/archive/2026-07-issue-27-transition-vectors-oom/task_plan.md

@@ -0,0 +1,27 @@
+# Task: dft_transition_vectors OOMs on large-extent rasters (processes full grid per class) (#27)
+
+`dft_transition_vectors()` processes the **entire raster grid** rather than the ~2% of
+cells that hold data, and does so **once per transition class**. On a large-extent
+floodplain this exhausts memory and the R process is OOM-killed.
+
+## Phase 1: Single-pass patches(values = TRUE) rewrite
+
+- [x] Replace the per-class/per-patch loop (`R/dft_transition_vectors.R:68-125`) with single-pass `terra::patches(x, directions = 8, values = TRUE)` → `terra::as.polygons()` → sparse pid→label map via `terra::cells()` + `terra::extract()`; keep signature, validation, `patch_area_min` filter, zones block, and return columns unchanged; drop rows whose value has no `cats()` entry (old behavior)
+- [x] Add `terra (>= 1.8-10)` floor to DESCRIPTION Imports (patches(values=TRUE) edge-wraparound bug fixed in 1.8-10 — correctness, not nicety)
+- [x] New tests in `tests/testthat/test-dft_transition_vectors.R`: synthetic decomposition test (~60x80 factor raster with engineered topology — diagonal-only joins, adjacent different classes, crossed diagonals, a code-0 class, holes) asserting equal patch count / per-class counts / sorted areas vs a brute-force per-class reference computed in-test; all-NA raster returns empty sf with correct columns + CRS; fixture regression guard (185 patches / 123.11 total ha)
+- [x] All 10 existing tests stay green unchanged (behavior contract)
+
+## Phase 2: Docs + release
+
+- [x] Roxygen: note single-pass implementation and that `patch_id` numbering is scan-order; `devtools::document()`
+- [x] NEWS.md 0.2.4: OOM fix (grid-cell x n-class x n-patch churn → single pass), identical output guarantee, patch_id ordering note, new terra floor
+- [x] File follow-up issue for `dft_rast_transition()`'s full-grid value pipeline (same OOM class; needs its own refactor + harness) — filed as #28
+- [x] `lintr` clean + full `devtools::test()` pass (211 pass; one pre-existing vignette lint, untouched by this branch)
+- [x] Version bump to 0.2.4 in DESCRIPTION as final commit
+
+## Validation
+
+- [x] Tests pass
+- [x] `/code-check` clean on each commit
+- [x] PWF checkboxes match landed work
+- [x] `/planning-archive` on completion

tests/testthat/test-dft_transition_vectors.R

@@ -107,3 +107,87 @@ test_that("errors when zones supplied without zone_col", {
   expect_error(dft_transition_vectors(result$raster, zones = aoi),
                "zone_col")
 })
+
+# builds a small projected factor raster from an integer matrix; codes become
+# levels labelled "class_<code>"
+transition_test_rast <- function(m) {
+  r <- terra::rast(m)
+  terra::ext(r) <- c(0, ncol(m) * 10, 0, nrow(m) * 10)
+  terra::crs(r) <- "EPSG:32609"
+  codes <- sort(unique(as.vector(m)))
+  terra::set.cats(
+    r, layer = 1,
+    value = data.frame(id = codes, transition = paste0("class_", codes))
+  )
+  r
+}
+
+test_that("patch decomposition matches per-class brute-force reference", {
+  m <- matrix(NA_integer_, nrow = 60, ncol = 80)
+  m[5:7, 5:7] <- 10L    # two class-10 blocks touching only at a corner:
+  m[8:10, 8:10] <- 10L  # 8-connectivity must merge them into one patch
+  m[5:7, 8:10] <- 20L   # class-20 block edge-adjacent to class 10: must stay separate
+  m[20, 20] <- 30L      # crossed diagonals: 30 and 40 each diagonal pairs
+  m[21, 21] <- 30L      # crossing each other — each class merges to one patch
+  m[20, 21] <- 40L
+  m[21, 20] <- 40L
+  m[40:42, 40:45] <- 0L # code 0 is a real class, not background
+  m[30:36, 60:66] <- 50L      # class-50 ring ...
+  m[32:34, 62:64] <- NA_integer_  # ... around a hole
+  m[50, 70] <- 20L      # isolated second class-20 patch
+
+  x <- transition_test_rast(m)
+  patches <- dft_transition_vectors(x)
+
+  # brute-force reference: per-class binary mask -> patches(directions = 8)
+  codes <- sort(unique(as.vector(m)))
+  for (code in codes) {
+    r_class <- terra::classify(transition_test_rast(m), cbind(code, 1),
+                               others = NA)
+    p_ref <- terra::patches(r_class, directions = 8)
+    ref_cells <- sort(terra::freq(p_ref)$count)
+
+    got <- patches[patches$transition == paste0("class_", code), ]
+    got_cells <- sort(round(got$area_ha * 1e4 / 100))  # 10 m cells -> cell counts
+    expect_equal(got_cells, ref_cells, label = paste("class", code))
+  }
+
+  expect_equal(length(unique(patches$patch_id)), nrow(patches))
+  # engineered expectations, independent of the reference implementation
+  n_by_class <- table(patches$transition)
+  expect_equal(unname(n_by_class[["class_10"]]), 1L)  # diagonal merge
+  expect_equal(unname(n_by_class[["class_20"]]), 2L)  # class boundary held
+  expect_equal(unname(n_by_class[["class_30"]]), 1L)  # crossed diagonal
+  expect_equal(unname(n_by_class[["class_40"]]), 1L)
+  expect_equal(unname(n_by_class[["class_0"]]), 1L)   # code 0 not background
+})
+
+test_that("all-NA raster returns empty sf with expected columns and CRS", {
+  m <- matrix(NA_integer_, nrow = 10, ncol = 10)
+  r <- terra::rast(m)
+  terra::ext(r) <- c(0, 100, 0, 100)
+  terra::crs(r) <- "EPSG:32609"
+  terra::set.cats(r, layer = 1,
+                  value = data.frame(id = 1L, transition = "class_1"))
+
+  out <- dft_transition_vectors(r)
+
+  expect_s3_class(out, "sf")
+  expect_equal(nrow(out), 0)
+  expect_true(all(c("patch_id", "transition", "area_ha") %in% names(out)))
+  expect_equal(sf::st_crs(out), sf::st_crs(r))
+})
+
+test_that("fixture decomposition is stable (regression guard)", {
+  r17 <- terra::rast(system.file("extdata", "example_2017.tif", package = "drift"))
+  r20 <- terra::rast(system.file("extdata", "example_2020.tif", package = "drift"))
+  classified <- dft_rast_classify(list("2017" = r17, "2020" = r20), source = "io-lulc")
+  result <- dft_rast_transition(classified, from = "2017", to = "2020")
+
+  patches <- dft_transition_vectors(result$raster)
+
+  expect_equal(nrow(patches), 185L)
+  expect_equal(round(sum(patches$area_ha), 2), 123.11)
+  expect_equal(nrow(dft_transition_vectors(result$raster, patch_area_min = 1000)),
+               57L)
+})