diff --git a/DESCRIPTION b/DESCRIPTION index 1afb974..449a7da 100644 --- a/DESCRIPTION +++ b/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", , "al@newgraphenvironment.com", 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, diff --git a/NEWS.md b/NEWS.md index 28edc82..a07e71d 100644 --- a/NEWS.md +++ b/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. diff --git a/R/dft_transition_vectors.R b/R/dft_transition_vectors.R index 108a7c7..3fc5c92 100644 --- a/R/dft_transition_vectors.R +++ b/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 diff --git a/man/dft_transition_vectors.Rd b/man/dft_transition_vectors.Rd index 95d568b..d20c0d3 100644 --- a/man/dft_transition_vectors.Rd +++ b/man/dft_transition_vectors.Rd @@ -22,7 +22,8 @@ smaller than this are dropped before returning. \code{NULL} (default) keeps all. \value{ An \code{sf} data frame (polygon geometry) with columns: \itemize{ -\item \code{patch_id} (integer) — connected component ID +\item \code{patch_id} (integer) — connected component ID, numbered in raster +scan order \item \code{transition} (character) — transition label (e.g. "Trees -> Rangeland") \item \code{area_ha} (numeric) — patch area in hectares \item Zone column (if \code{zones} supplied) — from spatial intersection @@ -34,6 +35,11 @@ polygons — one row per connected patch of pixels sharing the same transition type. Useful for QA in GIS, spatial attribution to management zones, and patch-level reporting. } +\details{ +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. +} \examples{ r17 <- terra::rast(system.file("extdata", "example_2017.tif", package = "drift")) r20 <- terra::rast(system.file("extdata", "example_2020.tif", package = "drift")) diff --git a/planning/archive/2026-07-issue-27-transition-vectors-oom/README.md b/planning/archive/2026-07-issue-27-transition-vectors-oom/README.md new file mode 100644 index 0000000..ffbf3ed --- /dev/null +++ b/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`) diff --git a/planning/archive/2026-07-issue-27-transition-vectors-oom/findings.md b/planning/archive/2026-07-issue-27-transition-vectors-oom/findings.md new file mode 100644 index 0000000..0e97ef5 --- /dev/null +++ b/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. diff --git a/planning/archive/2026-07-issue-27-transition-vectors-oom/progress.md b/planning/archive/2026-07-issue-27-transition-vectors-oom/progress.md new file mode 100644 index 0000000..8461bc6 --- /dev/null +++ b/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 diff --git a/planning/archive/2026-07-issue-27-transition-vectors-oom/task_plan.md b/planning/archive/2026-07-issue-27-transition-vectors-oom/task_plan.md new file mode 100644 index 0000000..0545f7d --- /dev/null +++ b/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 diff --git a/tests/testthat/test-dft_transition_vectors.R b/tests/testthat/test-dft_transition_vectors.R index 6bcde52..84ed8e8 100644 --- a/tests/testthat/test-dft_transition_vectors.R +++ b/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_" +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) +})