diff --git a/DESCRIPTION b/DESCRIPTION index e536cd9..7c1561a 100644 --- a/DESCRIPTION +++ b/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"), diff --git a/NEWS.md b/NEWS.md index 6d1cbfa..3ebeae9 100644 --- a/NEWS.md +++ b/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. diff --git a/R/dft_stac_cube.R b/R/dft_stac_cube.R index 711ebc0..f00d921 100644 --- a/R/dft_stac_cube.R +++ b/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,6 +282,13 @@ 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) @@ -275,19 +296,36 @@ dft_stac_cube <- function(aoi, } +#' 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( @@ -295,7 +333,7 @@ stac_cube_cache_key <- function(aoi_target, res, target_crs, dt, aggregation, 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 ) diff --git a/inst/notes/gdalcubes-pc-gotchas.md b/inst/notes/gdalcubes-pc-gotchas.md index 5cbd73b..ef8e9e8 100644 --- a/inst/notes/gdalcubes-pc-gotchas.md +++ b/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 diff --git a/man/dft_stac_cube.Rd b/man/dft_stac_cube.Rd index df566fe..c64e3aa 100644 --- a/man/dft_stac_cube.Rd +++ b/man/dft_stac_cube.Rd @@ -14,6 +14,7 @@ dft_stac_cube( dt = "P1M", aggregation = "median", resampling = "bilinear", + clip = TRUE, cloud_cover_max = 60, months = NULL, mask_values = NULL, @@ -48,6 +49,14 @@ must agree with.} \item{resampling}{Character. Spatial resampling (default \code{"bilinear"}).} +\item{clip}{Logical. When \code{TRUE} (default), clip the returned stack to the AOI +polygon with \code{terra::mask()} (cells outside → \code{NA} on every layer), so +\code{\link[=dft_rast_break]{dft_rast_break()}} / \code{\link[=dft_rast_trend]{dft_rast_trend()}} reduce only in-polygon pixels. Set +\code{FALSE} to keep the full bounding box (e.g. for surrounding context, or to +mask later with a different polygon). Note this clips the \emph{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 \code{inst/notes/gdalcubes-pc-gotchas.md}).} + \item{cloud_cover_max}{Numeric. Scene-level \code{eo:cloud_cover} maximum percent for the STAC pre-filter (default 60).} @@ -75,12 +84,13 @@ cloud / shadow / cirrus classes).} } \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. +time value per layer — cached as a GeoTIFF. By default (\code{clip = TRUE}) the +stack is clipped to the AOI polygon (cloud-masked, cells outside the polygon +\code{NA}), so the reduced raster from \code{\link[=dft_rast_break]{dft_rast_break()}} is already polygon-tight; +pass \code{clip = FALSE} for the full AOI \strong{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. } \description{ Sibling of \code{\link[=dft_stac_fetch]{dft_stac_fetch()}} for continuous change detection. Where diff --git a/planning/archive/2026-07-issue-32-stac-cube-aoi-clip/README.md b/planning/archive/2026-07-issue-32-stac-cube-aoi-clip/README.md new file mode 100644 index 0000000..ae896b0 --- /dev/null +++ b/planning/archive/2026-07-issue-32-stac-cube-aoi-clip/README.md @@ -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`. diff --git a/planning/archive/2026-07-issue-32-stac-cube-aoi-clip/findings.md b/planning/archive/2026-07-issue-32-stac-cube-aoi-clip/findings.md new file mode 100644 index 0000000..fd8d851 --- /dev/null +++ b/planning/archive/2026-07-issue-32-stac-cube-aoi-clip/findings.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). diff --git a/planning/archive/2026-07-issue-32-stac-cube-aoi-clip/progress.md b/planning/archive/2026-07-issue-32-stac-cube-aoi-clip/progress.md new file mode 100644 index 0000000..b99dfa8 --- /dev/null +++ b/planning/archive/2026-07-issue-32-stac-cube-aoi-clip/progress.md @@ -0,0 +1,35 @@ +# Progress — dft_stac_cube AOI-polygon clip (#32) + +## Session 2026-07-09 + +- Plan-mode exploration — traced `dft_stac_cube` build path, the `filter_geom` + removal (#30, commit c6953a0), the downstream NA-skip gate in + `dft_rast_break`/`dft_rast_trend`, the proven `dft_stac_fetch.R:150` mask pattern, + and the `gotchas` note. Confirmed example AOI non-rectangular (area/bbox ≈ 0.105). +- Plan-agent review — reframed the rationale (compute win illusory; real win is + polygon-tight output + no caller-side mask), flagged the cache-key wiring Blocker, + simplified the helper. Phases approved by user. +- Archived #34 PWF confirmed complete; #34 plan-mode scratch doc left to be overwritten + (user decision — it's redundant with the archived findings/task_plan). +- Created branch `32-dft-stac-cube-restore-aoi-polygon-clip-f` off main. +- Scaffolded PWF baseline with approved phases. +- Phase 1 (commit pending): wrote the test contract — `clip` threaded through the + `cube_key()` helper + `cube_key(clip = FALSE) != base`; offline `stac_cube_clip()` + masking test (synthetic raster + half-covering polygon); network e2e high-NA-fraction + assertion under default clip. Confirmed red against current code (FAIL 6 | SKIP 2 | + PASS 1): 5 cache-key tests error on `unused argument (clip)`, offline test on + `object 'stac_cube_clip' not found`. `/code-check` deferred to the Phase 2 impl diff + (Phase 1 is test-only, self-reviewed). +- Phase 2 (commit pending): confirmed `dft_stac_fetch` normalizes `aoi` identically and + `terra::vect()` accepts `sfc` → bare `terra::vect(aoi_target)` matches the sibling. + Added `clip = TRUE` param, `stac_cube_clip()` `@noRd` helper, the mask step after stk + assembly, `clip` in `stac_cube_cache_key` (hash) + call site, roxygen + `filter_geom` + comment updates. `/code-check` (1 fresh-eyes round) found a fragile mismatch — `isTRUE` + gate vs `as.logical` key disagree for `clip = 1`/`"TRUE"`; fixed by normalizing + `clip <- isTRUE(as.logical(clip))` once up front. Full suite 319 pass, lint clean. +- Phase 3 (commit pending): recorded the #32 resolution in the gotchas note + (clip restored via cube-level `terra::mask()`; residual = fetch-time streaming + still bbox-bound); added the `NEWS.md` 0.5.0 entry (behavior change + cache + rebuild + Closes #32). DESCRIPTION bump deferred to Phase 4 as the final + `Release v0.5.0` commit (r-packages "bump last" rule; matches #34). +- Next: Phase 4 — validate (check/test), Release v0.5.0 bump, archive, PR, merge. diff --git a/planning/archive/2026-07-issue-32-stac-cube-aoi-clip/task_plan.md b/planning/archive/2026-07-issue-32-stac-cube-aoi-clip/task_plan.md new file mode 100644 index 0000000..d062ac5 --- /dev/null +++ b/planning/archive/2026-07-issue-32-stac-cube-aoi-clip/task_plan.md @@ -0,0 +1,45 @@ +# Task: dft_stac_cube — restore AOI-polygon clip (filter_geom segfault workaround) (#32) + +`dft_stac_cube()` originally clipped its STAC index cube to the AOI polygon with +`gdalcubes::filter_geom()`, which segfaults / returns an all-NA cube on the pinned +gdalcubes 0.7.3 build, so #30 removed it — the cube now spans the AOI **bounding +box**. Restore polygon-tight computation *without* `filter_geom` by masking the +assembled terra stack client-side (`terra::mask()`, the proven `dft_stac_fetch.R:150` +pattern) behind a new `clip = TRUE` default. 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. + +Honest scope: the dominant cost (streaming the full bbox of COGs via +`cube_view(extent = bbox_target)`) is unchanged — a post-hoc mask can't push the AOI +into the read. Real wins: polygon-tight output by default + a modest reducer speedup. +Fetch-time streaming stays bbox-bound (documented residual). + +## Phase 1: Test contract (failing first) +- [x] `cube_key()` helper gains `clip`; `cube_key(clip = FALSE) != base` assertion added +- [x] `stac_cube_clip()` offline masking unit test (synthetic raster + polygon) +- [x] opt-in network e2e asserts a bbox-corner cell is all-`NA` under default clip +- [x] confirm the new tests FAIL against current code (FAIL 6 | SKIP 2 | PASS 1) + +## Phase 2: Implement clip +- [x] `clip = TRUE` param + `stac_cube_clip()` helper + apply after stk assembly +- [x] `clip` threaded into `stac_cube_cache_key()` (hash) AND the call site +- [x] roxygen `@param clip` + `@return` rewrite + `filter_geom` comment updates +- [x] `devtools::document()`; `lintr::lint_package()` clean; `devtools::test()` green (319 pass) +- [x] `/code-check`: normalize `clip <- isTRUE(as.logical(clip))` once so the mask gate + and cache key can't disagree for truthy-but-non-TRUE inputs (silent wrong-extent fix) + +## Phase 3: Docs + gotchas note + NEWS +- [x] update `inst/notes/gdalcubes-pc-gotchas.md` filter_geom bullet (#32 resolved + residual) +- [x] `NEWS.md` 0.5.0 entry (behavior change + cache rebuild + Closes #32) + +## Phase 4: Validate, bump, archive, PR, release +- [x] `devtools::test()` / `lint` / `document` / `check` clean (check: 0E/0W/0N; suite 319 pass) +- [x] `DESCRIPTION` 0.4.0 → 0.5.0 + Date, committed as `Release v0.5.0` (final code commit) +- [ ] `/planning-archive`; `/gh-pr-push` (`Fixes #32`, `Relates to NewGraphEnvironment/sred-2025-2026#16`) +- [ ] `/gh-pr-merge` → tag v0.5.0 + +## Validation +- [ ] Tests pass (`devtools::test()`); network/bfast tests skip cleanly +- [ ] `/code-check` clean on each commit +- [ ] PWF checkboxes match landed work +- [ ] `/planning-archive` on completion diff --git a/tests/testthat/test-dft_stac_cube.R b/tests/testthat/test-dft_stac_cube.R index 2ec39b3..efe3ea2 100644 --- a/tests/testthat/test-dft_stac_cube.R +++ b/tests/testthat/test-dft_stac_cube.R @@ -36,11 +36,11 @@ cube_key <- function(aoi = square_aoi(), res = 10, target_crs = "EPSG:32609", 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) { + offset_before = 0, clip = TRUE) { 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 + months, offset_before, clip ) } @@ -67,6 +67,9 @@ test_that("stac_cube_cache_key changes with each cube-affecting parameter", { 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) + # clip must key distinctly, or a clip=FALSE request silently hits a clipped + # (or vice-versa) cached .tif and returns wrong-extent data + expect_false(cube_key(clip = FALSE) == base) }) test_that("stac_cube_cache_key normalizes months order", { @@ -85,6 +88,30 @@ test_that("stac_cube_cache_key ignores sf attribute columns", { expect_equal(cube_key(with_attrs), cube_key(bare)) }) +# stac_cube_clip(): the AOI-polygon clip that replaces gdalcubes::filter_geom +# (#32). Network-free — a synthetic stack + a half-covering polygon. Cells whose +# centre is outside the polygon become NA on every layer; nlyr is preserved. +test_that("stac_cube_clip masks cells outside the AOI polygon on every layer", { + r <- terra::rast(nrows = 10, ncols = 10, xmin = 0, xmax = 10, + ymin = 0, ymax = 10, crs = "EPSG:32609", nlyrs = 2) + terra::values(r) <- 1 # every cell valid on both layers + # AOI covers the left half (x in [0, 5]); no cell centre lands exactly on 5 + aoi <- sf::st_sfc( + sf::st_polygon(list(rbind( + c(0, 0), c(5, 0), c(5, 10), c(0, 10), c(0, 0) + ))), + crs = 32609 + ) + out <- drift:::stac_cube_clip(r, aoi) + + expect_s4_class(out, "SpatRaster") + expect_equal(terra::nlyr(out), 2L) # layers preserved + vals <- terra::values(out) + inside <- terra::xyFromCell(out, seq_len(terra::ncell(out)))[, 1] < 5 + expect_true(all(!is.na(vals[inside, ]))) # inside polygon: retained + expect_true(all(is.na(vals[!inside, ]))) # outside polygon: NA +}) + # 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", { @@ -103,6 +130,10 @@ test_that("dft_stac_cube fetches an index stack end-to-end", { expect_s4_class(cube, "SpatRaster") expect_equal(terra::nlyr(cube), 3) # 3 monthly layers expect_false(anyNA(terra::time(cube))) # time set per layer + # default clip = TRUE clips to the AOI polygon: for this thin reach + # (area / bbox ~= 0.105) most bbox cells are fully NA across all layers + fully_na <- mean(rowSums(!is.na(terra::values(cube))) == 0) + expect_gt(fully_na, 0.5) # 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)