Merge pull request #37 from NewGraphEnvironment/32-dft-stac-cube-restore-aoi-polygon-clip-f

Restore dft_stac_cube() AOI-polygon clip via terra::mask (#32)

Author: NewGraphEnvironment

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: drift
 Title: Detecting Riparian and Inland Floodplain Transitions
-Version: 0.4.0
+Version: 0.5.0
 Date: 2026-07-09
 Authors@R: c(
     person("Allan", "Irvine", , "al@newgraphenvironment.com", role = c("aut", "cre"),

NEWS.md

@@ -1,3 +1,7 @@
+# drift 0.5.0
+
+- `dft_stac_cube()` gains `clip` (default `TRUE`), restoring AOI-polygon-tight output (#32). The assembled index stack is masked to the AOI polygon with `terra::mask()` — client-side, because `gdalcubes::filter_geom()` segfaults / returns an all-NA cube on the pinned build — so cells outside the polygon are `NA` on every layer. The reduced raster from `dft_rast_break()`/`dft_rast_trend()` is now polygon-tight with no caller-side mask, and those reducers skip out-of-AOI pixels via their valid-observation gate. `clip = FALSE` keeps the full bounding box. This is an output change for callers that relied on the bounding-box extent, and the clip is folded into the cube cache key, so existing cached cubes rebuild once. Note the clip affects the *output* only — the full bbox of COGs is still streamed either way (the AOI cannot be pushed into the read on the pinned gdalcubes build).
+
 # 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.

R/dft_stac_cube.R

@@ -37,6 +37,13 @@
 #' @param aggregation Character. Temporal aggregation for multiple scenes in one
 #'   `dt` window (default `"median"`).
 #' @param resampling Character. Spatial resampling (default `"bilinear"`).
+#' @param clip Logical. When `TRUE` (default), clip the returned stack to the AOI
+#'   polygon with `terra::mask()` (cells outside → `NA` on every layer), so
+#'   [dft_rast_break()] / [dft_rast_trend()] reduce only in-polygon pixels. Set
+#'   `FALSE` to keep the full bounding box (e.g. for surrounding context, or to
+#'   mask later with a different polygon). Note this clips the *output* only — the
+#'   full bbox of COGs is still streamed either way (the AOI cannot be pushed into
+#'   the read on the pinned gdalcubes build; see `inst/notes/gdalcubes-pc-gotchas.md`).
 #' @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`
@@ -58,12 +65,13 @@
 #'   [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.
+#'   time value per layer — cached as a GeoTIFF. By default (`clip = TRUE`) the
+#'   stack is clipped to the AOI polygon (cloud-masked, cells outside the polygon
+#'   `NA`), so the reduced raster from [dft_rast_break()] is already polygon-tight;
+#'   pass `clip = FALSE` for the full AOI **bounding box**. 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
@@ -93,6 +101,7 @@ dft_stac_cube <- function(aoi,
                           dt = "P1M",
                           aggregation = "median",
                           resampling = "bilinear",
+                          clip = TRUE,
                           cloud_cover_max = 60,
                           months = NULL,
                           mask_values = NULL,
@@ -137,6 +146,10 @@ dft_stac_cube <- function(aoi,
   # each side, so a series crossing it has no artificial index step.
   offset_boundary <- cfg$offset_boundary
   offset_before <- cfg$offset_before %||% 0
+  # normalize clip to a single scalar so the mask gate and the cache key agree: a
+  # truthy-but-non-TRUE clip (e.g. 1 or "TRUE") must not skip the mask yet key as
+  # TRUE, which would let a later clip=TRUE read the unclipped cube (#32).
+  clip <- isTRUE(as.logical(clip))
 
   # Ensure aoi is sf
   if (inherits(aoi, "SpatVector")) aoi <- sf::st_as_sf(aoi)
@@ -168,7 +181,7 @@ dft_stac_cube <- function(aoi,
   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
+    cloud_cover_max, mask_values, scale, offset, months, offset_before, clip
   )
   cache_file <- file.path(cache_source_dir, paste0("cube_", cache_key, ".tif"))
 
@@ -190,8 +203,8 @@ dft_stac_cube <- function(aoi,
     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)
+      # cube extent and the terra::mask() 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
@@ -229,8 +242,9 @@ dft_stac_cube <- function(aoi,
   # 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.
+  # crash the compute worker) on the pinned build, so we mask clouds here and clip
+  # the assembled stack to the AOI polygon afterward with terra::mask() (see
+  # stac_cube_clip, #32), as the sibling dft_stac_fetch() does — never filter_geom().
   build_index_stack <- function(features, offset_use) {
     img_col <- gdalcubes::stac_image_collection(
       features, asset_names = c(band_assets, mask_asset)
@@ -268,34 +282,58 @@ dft_stac_cube <- function(aoi,
     stk <- build_index_stack(items$features, if (all(is_pre)) offset_before else offset)
   }
 
+  # Restore the AOI-polygon clip removed in #30: mask the assembled stack (never
+  # gdalcubes::filter_geom, which segfaults on the pinned build). Out-of-polygon
+  # cells become NA on every layer, so dft_rast_break()/dft_rast_trend() skip them
+  # via their `rowSums(!is.na) >= min_obs` gate. `mask` preserves nlyr and time is
+  # set below, so the cached tif — and the cache-read path — need no other change.
+  if (isTRUE(clip)) stk <- stac_cube_clip(stk, aoi_target)
+
   terra::time(stk) <- month_times(terra::nlyr(stk))
   names(stk) <- rep(index, terra::nlyr(stk))
   terra::writeRaster(stk, cache_file, overwrite = TRUE)
   stk
 }
 
 
+#' Clip an index stack to the AOI polygon (client-side terra mask)
+#'
+#' Restores AOI-polygon-tight output without `gdalcubes::filter_geom()`, which
+#' segfaults / returns an all-NA cube on the pinned build (see
+#' `inst/notes/gdalcubes-pc-gotchas.md`, drift#32). Cells whose centre falls
+#' outside the polygon become `NA` on every layer, so [dft_rast_break()] /
+#' [dft_rast_trend()] skip them via their `rowSums(!is.na) >= min_obs` gate. A
+#' multi-feature `aoi` masks to the union. Mirrors the post-read mask in
+#' [dft_stac_fetch()]; `aoi` is already in the stack's CRS.
+#' @noRd
+stac_cube_clip <- function(stk, aoi) {
+  terra::mask(stk, terra::vect(aoi))
+}
+
+
 #' 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.
+#' `scale`/`offset` are included because they change pixel values. `clip` is
+#' included because it changes the written extent (polygon vs bbox), so a
+#' `clip = FALSE` request must not read a clipped cube (or vice versa).
 #' @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) {
+                                offset_before = 0, clip = TRUE) {
   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)
+      as.numeric(offset_before), as.logical(clip)
     )),
     1, 12
   )

inst/notes/gdalcubes-pc-gotchas.md

@@ -8,8 +8,15 @@ bfast 1.7.2.
 - **`gdalcubes::filter_geom()` segfaults / returns an all-NA cube** on this build
   (crashes `gc_exec_worker`, `address 0x120`). Do NOT clip to an AOI polygon
   inside the cube pipeline. Use the AOI bbox in `cube_view(extent=)` and
-  `terra::mask()` the reduced raster afterward (as `dft_stac_fetch()` does).
-  Tracked as drift#32.
+  `terra::mask()` afterward, as `dft_stac_fetch()` does. **Resolved (#32):**
+  `dft_stac_cube(clip = TRUE)` (the default) masks the assembled terra stack to
+  the AOI polygon client-side (helper `stac_cube_clip()` = `terra::mask(stk,
+  terra::vect(aoi))`), so the cube is polygon-tight and
+  `dft_rast_break()`/`dft_rast_trend()` skip out-of-AOI pixels via their
+  `rowSums(!is.na) >= min_obs` gate. **Residual:** this clips the *output* only —
+  `cube_view(extent = bbox)` still streams the full bbox of COGs, so fetch time is
+  unchanged; pushing the AOI into the read would need a working `filter_geom` or
+  server-side windowing. `clip = FALSE` keeps the full bbox.
 - **`reduce_time()` R-callback runs in spawned worker processes at EVERY parallel
   setting** (incl. `parallel = 1`). A closure over enclosing locals fails there
   (`object 'band' not found`). Options: build a self-contained callback (inline

man/dft_stac_cube.Rd

@@ -0,0 +1,26 @@
+## Outcome
+
+Restored `dft_stac_cube()`'s AOI-polygon clip (#32), removed in #30 when
+`gdalcubes::filter_geom()` proved to segfault / return an all-NA cube on the
+pinned gdalcubes 0.7.3 build. The clip is now done client-side: a new
+`clip = TRUE` default masks the assembled terra index stack to the AOI polygon
+with `terra::mask()` (helper `stac_cube_clip()`, mirroring `dft_stac_fetch()`).
+Out-of-polygon cells become `NA` on every layer, so `dft_rast_break()` /
+`dft_rast_trend()` skip them via their existing `rowSums(!is.na) >= min_obs` gate,
+and the reduced raster is polygon-tight with no caller-side mask. `clip = FALSE`
+keeps the full bbox.
+
+Key learnings: the headline cost — streaming the full bbox of COGs via
+`cube_view(extent = bbox)` — is **unchanged** (a post-hoc mask can't push the AOI
+into the read), so the genuine wins are polygon-tight output + a modest reducer
+speedup, not a fetch-time saving; that residual is documented in the gotchas note.
+The load-bearing correctness detail was threading `clip` into **both** the mask
+gate and the cache key, and normalizing `clip <- isTRUE(as.logical(clip))` once so
+a truthy-but-non-`TRUE` input (e.g. `1`, `"TRUE"`) can't skip the mask yet key as
+`TRUE` — which a `/code-check` fresh-eyes pass caught.
+
+Released as **v0.5.0**. `devtools::check()` clean (0E/0W/0N); full suite 319 pass;
+offline `stac_cube_clip()` masking test + `cube_key(clip=FALSE) != base` cover the
+contract network-free.
+
+Closed by: PR (Fixes #32) on branch `32-dft-stac-cube-restore-aoi-polygon-clip-f`.

planning/archive/2026-07-issue-32-stac-cube-aoi-clip/README.md

@@ -0,0 +1,66 @@
+# Findings — dft_stac_cube AOI-polygon clip (#32)
+
+## Issue context
+
+Restore AOI-polygon clipping in `dft_stac_cube()`. It was intended to clip the cube
+to the AOI polygon with `gdalcubes::filter_geom()`; on gdalcubes 0.7.3 that yields an
+entirely-NA cube and can segfault the compute worker, so #30 removed it. The cube now
+spans the AOI **bounding box** and callers clip the reduced raster with
+`terra::mask()`. Cost: `dft_rast_break()` reduces over the whole bbox rather than just
+the floodplain polygon (a few× more pixels for a thin reach). #32 scoped to the AOI
+clip only (the Sentinel-2 baseline-offset half shipped in #30).
+
+## Exploration (2026-07-09)
+
+- **Downstream consumers already skip NA pixels.** `dft_rast_break.R:101-102` and
+  `dft_rast_trend.R:58-59` both do `vals <- terra::values(cube); usable <-
+  which(rowSums(!is.na(vals)) >= min_obs)` and run the expensive per-pixel reducer
+  only on `usable` rows. So masking the cube to the AOI (out-of-polygon → NA on every
+  layer) makes those pixels drop out of the reducer for free — no change to break/trend.
+
+- **Proven in-package clip pattern:** `dft_stac_fetch.R:150` —
+  `terra::mask(r, terra::vect(aoi_target))` — applied client-side after reading. Mirror
+  it at the cube-stack level in `dft_stac_cube()`.
+
+- **`filter_geom` gotcha documented:** `inst/notes/gdalcubes-pc-gotchas.md:8-12`
+  records the segfault (`gc_exec_worker`, `address 0x120`) and tags the fix as #32,
+  recommending bbox `cube_view` + `terra::mask()`.
+
+- **Example AOI is non-rectangular:** single MULTIPOLYGON, 2049 pts, area/bbox ≈ 0.105
+  → a clipped cube has all-NA bbox corners (valid opt-in network assertion).
+
+## Plan-agent review (key corrections adopted)
+
+1. **Compute win is largely illusory — reframe.** `cube_view(extent = bbox_target)`
+   (`dft_stac_cube.R:218-227`) still streams the full bbox of COGs (the ~10-30 min
+   cost per `gotchas:40-41`); `terra::values()` still loads the full bbox matrix (peak
+   memory unchanged). Only the seconds-scale per-pixel reducer loop shrinks. Lead the
+   rationale with **polygon-tight output + no caller-side mask needed**; describe the
+   reducer speedup as a modest secondary effect. Fetch-time pixel savings that
+   `filter_geom` would have given are **not** recoverable — documented residual.
+
+2. **Blocker — thread `clip` into BOTH the mask step and the cache key.** The key fn
+   has a default, so forgetting `clip` at the call site (`:168-172`) makes `clip=TRUE`
+   and `clip=FALSE` hash identically → silent wrong-extent cache hit (no error). Add
+   `expect_false(cube_key(clip = FALSE) == cube_key(clip = TRUE))`.
+
+3. **Simplify the helper** to bare `terra::vect(aoi_target)`, matching fetch — drop the
+   `sf::st_as_sf()` coercion unless fetch's path proves `sfc` needs it. `terra::mask`
+   with a multi-feature SpatVector masks to the union (matches the STAC `st_union`
+   query at `:195`).
+
+4. **Ordering is safe:** mask the `terra::cover(...)` result (offset-split branch) or
+   the single-build result before `time`/`names`/`writeRaster`. `mask` preserves
+   `nlyr`; `time` is re-derived on every read via `month_times(terra::nlyr(r))`
+   (`:183`), so the cache-read path needs no change.
+
+5. **Behavior change + cache churn:** direct cube users now get NA outside the AOI by
+   default; hashing `clip` invalidates every existing 0.4.0 bbox cache (one-time
+   rebuild). Both flagged in NEWS. Considered a `_clip.tif` filename suffix to spare
+   `clip=FALSE` users the re-fetch — rejected as needless branching for a pre-1.0
+   package with ~one cube-cache user.
+
+6. **Vignette shielded:** `vignettes/trajectory-break-detection.Rmd` loads a committed
+   `.rds` (pipeline chunk `eval = FALSE`), so the default change doesn't alter the
+   rendered vignette. `data-raw/vignette_data_break.R`'s `terra::mask` becomes an
+   idempotent no-op under `clip = TRUE` — leave it (defensive).