From 127de941688691ad3b4a62c7ccd113c745fadbdc Mon Sep 17 00:00:00 2001 From: almac2022 Date: Thu, 9 Jul 2026 15:26:01 -0700 Subject: [PATCH 1/7] Initialize PWF baseline for #36 Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01DKqnW6Fni3tc44Lo5D5fJz --- planning/active/findings.md | 79 ++++++++++++++++++++++++++++++++++++ planning/active/progress.md | 13 ++++++ planning/active/task_plan.md | 51 +++++++++++++++++++++++ 3 files changed, 143 insertions(+) create mode 100644 planning/active/findings.md create mode 100644 planning/active/progress.md create mode 100644 planning/active/task_plan.md diff --git a/planning/active/findings.md b/planning/active/findings.md new file mode 100644 index 0000000..1496c5f --- /dev/null +++ b/planning/active/findings.md @@ -0,0 +1,79 @@ +# Findings — Tile dft_stac_fetch (#36) + +## Issue context + +`dft_stac_fetch()` builds a single gdalcubes cube over the AOI's **bounding box** +and masks to the polygon afterward. For a thin, spread-out AOI — a floodplain +following a river corridor — the bounding box is largely empty, so the download +scales with the bbox, not the AOI. Large-floodplain fetches are download-bound +(~30 min per group at 10 m × 3 yr). + +Measured on an io-lulc 10 m floodplain fetch (NECR): the grid is 57.0 M cells, +of which 5.73 M (10.1%) fall inside the floodplain — a **~10× download overhead**. +Thinner / more diagonal reaches are worse. + +**Distinct from #32:** #32 restores polygon-tight *compute* in `dft_stac_cube()` +via `gdalcubes::filter_geom()`, which is blocked upstream (segfault in +`gc_exec_worker`). This issue is the `filter_geom`-independent path: cut the +**download** by tiling, on the categorical `dft_stac_fetch()` (io-lulc) path. + +**Proposed:** tile the AOI into a grid, fetch only tiles intersecting the AOI +(skip empty tiles), and mosaic. Peak download approaches the union of intersecting +tile bboxes — near the AOI footprint for a corridor, rather than the full bbox. +Tile size trades request count against per-request waste. + +Related: #32 (`filter_geom`, blocked), #34 (transition memory fixed; fetch now +the dominant cost), #38 (same residual on the cube path). + +## Design notes (from Plan-agent review, 2026-07-09) + +- **Mosaic == untiled over the AOI, not raster-identical.** `st_make_grid` + over-hangs `xmax/ymax` by up to `tile_size − remainder`, and `terra::mask()` + does not crop — so the tiled masked raster keeps a larger extent with a wider + NA margin. Downstream (`dft_rast_classify` on `layer=1`, `dft_rast_summarize` + reduce-by-value) ignores the NA margin. Test oracle must crop both to their + common AOI intersection before comparing values — NOT `all.equal` on raw + dimensions. Inter-tile alignment is exact as long as `tile_size` is snapped to + a multiple of `res` and the grid is anchored at `(xmin, ymin)`; gdalcubes + streams whatever source COG pixels each output window needs (incl. just outside + the tile edge), so even bilinear edges match — no seams. +- **BLOCKER → write the tiled mosaic as `.tif`, not `.nc`.** terra's NetCDF + *write* is fragile on the pinned stack (gotchas note: terra↔NetCDF round-trip, + layer naming, NoData differ). The sibling `dft_stac_cube` writes + `terra::writeRaster(stk, ...)` to `.tif` — proven precedent. Derive the cache + extension from `is.null(tile_size)` so lookup + writer agree. Downstream is + layer-name-agnostic (`dft_rast_classify` uses `layer=1`; `stac_items` attr set + on the list), so `.nc`-vs-`.tif` read differences don't matter. +- **Leave boundary tiles un-trimmed.** Clipping the max-edge tiles to the bbox + would make their span not a multiple of `res`, breaking congruence/alignment. + The `< tile_size` overhang is masked away — bounded waste, worth it. +- **Cache-key conditional-append is safe** iff one normalized `is.null(tile_size)` + predicate drives both the path gate and the key-append (satisfies the #32 + normalize-once convention). Snap `tile_size` BEFORE hashing so `504`/`500` + (res 10) → same key. Add a golden regression freezing `cache_key(NULL)` to the + current hash — that test is the guardian of legacy-cache preservation. +- **Guards:** validate `tile_size` is NULL or one positive finite numeric (abort + on `NA`/`0`/negative/`Inf`/non-scalar/non-numeric); after snapping require + `>= res` (guards `tile_size < res/2` → 0 → degenerate grid). Empty intersecting + set aborts. Single-tile (`tile_size ≥ bbox`) proceeds through the tiled path + (one tile) — do NOT reroute to untiled (would desync key/format). +- **GDAL `/vsicurl` config on the tiled path.** `dft_stac_fetch` doesn't set the + `GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR` / `VSI_CACHE` / HTTP-multiplex config + that `dft_stac_cube` sets. Tiling multiplies per-item COG opens, so set it (with + `on.exit` restore) or small tiles trade data volume for open latency. +- **tile_size units** = CRS units (metres for the default UTM CRS), same as `res`. +- **tempfile hygiene:** unique `tempfile()` per (year × tile); `unlink()` after + each year's merge. terra rasters are disk-backed + merge streams → bounded RAM. + +## Key code locations + +- `R/dft_stac_fetch.R:122-151` — per-year fetch loop (cube_view over bbox_target, + raster_cube, write_ncdf to `_.nc`, read back, mask). The download is + `write_ncdf`. +- `R/dft_stac_fetch.R:169-179` — `stac_cache_key()` (9 hashed params). +- `R/dft_stac_cube.R:231-294` — sibling: cube_view, `terra::cover` assembly, + `terra::writeRaster → .tif` cache, GDAL config, `st_union` query. Reference. +- `tests/testthat/test-dft_stac_fetch.R:43-72` — `cache_key()` helper + "changes + with each fetch-affecting parameter" block to extend. +- `tests/testthat/test-dft_stac_cube.R:115-140` — `DRIFT_TEST_NETWORK` gate to + mirror for the opt-in e2e. diff --git a/planning/active/progress.md b/planning/active/progress.md new file mode 100644 index 0000000..2675db8 --- /dev/null +++ b/planning/active/progress.md @@ -0,0 +1,13 @@ +# Progress — Tile dft_stac_fetch (#36) + +## Session 2026-07-09 + +- Plan-mode exploration — read `dft_stac_fetch.R`, the cube sibling's assembly, + the fetch test conventions, and the STAC config. Phases approved by user. +- Plan-agent design review caught: write tiled mosaic as `.tif` not `.nc` + (terra NetCDF write fragile); test oracle compares over cropped common AOI + cells not raw dimensions; guard `tile_size ≥ res` after snapping; golden hash + regression for `tile_size = NULL`; GDAL `/vsicurl` config on the tiled path. +- Created branch `36-tile-dft-stac-fetch-to-bound-download-over-spars` off main. +- Scaffolded PWF baseline with approved phases. +- Next: start Phase 1 (tests-first for `tile_grid()` + `tile_size` normalization). diff --git a/planning/active/task_plan.md b/planning/active/task_plan.md new file mode 100644 index 0000000..dd3ea54 --- /dev/null +++ b/planning/active/task_plan.md @@ -0,0 +1,51 @@ +# Task: Tile dft_stac_fetch to bound download over sparse-floodplain bounding boxes (#36) + +`dft_stac_fetch()` builds one gdalcubes cube over the AOI's **bounding box** and +masks to the polygon only after the pixels are streamed. For a thin, diagonal +floodplain corridor the bbox is largely empty, so the download scales with the +bbox, not the AOI (~10× overhead measured on a real io-lulc floodplain: 10.1% of +bbox cells inside the polygon). `gdalcubes::filter_geom()` (the polygon clip that +would fix this) segfaults on the pinned build, so it stays blocked. Fix: add an +opt-in `tile_size` that tiles the AOI bbox into a res-aligned grid, fetches only +tiles intersecting the AOI polygon, and mosaics. Default `NULL` = today's +single-cube behavior, existing caches preserved. Sibling cube-path residual is +tracked separately as #38 (out of scope). + +## Phase 1: `tile_grid()` + `tile_size` normalization (offline, tests first) +- [ ] validation aborts: `NA`, `0`, negative, `Inf`, `c(1,2)`, `"500"`, snapped `< res` +- [ ] snapping: non-multiple snaps to nearest multiple of `res`, messages, snapped value used +- [ ] intersecting subset: a diagonal/L-shaped AOI in a known bbox keeps exactly the expected tiles, drops empty ones, kept-count ≪ full grid (the efficiency mechanism, offline) +- [ ] res-alignment: every extent's `left-xmin`, width, height are exact multiples of `res`; first tile lower-left `== (xmin, ymin)`; single-tile when `tile_size ≥ bbox`; empty/degenerate AOI aborts +- [ ] implement `tile_grid()` + normalization; new tests green + +## Phase 2: cache key — conditional `tile_size` append (offline) +- [ ] golden regression: `cache_key(tile_size = NULL)` == frozen current 12-char hash for fixed inputs (guards legacy-cache preservation) +- [ ] `cache_key(tile_size = 500) != cache_key(tile_size = NULL)`; distinct sizes → distinct keys; snap-before-key: `504` and `500` (res 10) → same key +- [ ] `stac_cache_key()` gains `tile_size = NULL`; append only when non-NULL; call site passes `tile_size`; tests green + +## Phase 3: extract `fetch_extent_to()`, refactor untiled path (behavior-preserving) +- [ ] extract helper; untiled path routes through it, writing straight to `_.nc` — identical filename/format +- [ ] existing tests + opt-in untiled network test still green (no observable change) + +## Phase 4: tiled branch — mosaic assembly + wire `tile_size` end-to-end +- [ ] offline merge oracle: reference SpatRaster split into res-lattice tiles → `terra::merge()` → `all.equal(values(merged), values(reference))`; masked mosaic == masked reference over AOI; `.tif` round-trip preserves single-layer integer codes +- [ ] add `tile_size` to `dft_stac_fetch`; tiled branch (per-tile fetch → merge → `.tif` → read → mask); extension from `is.null(tile_size)`; GDAL config + `on.exit`; tempfile `unlink` +- [ ] roxygen `@param tile_size` + `\dontrun` example + cache-doc `.tif` note +- [ ] opt-in network e2e (`DRIFT_TEST_NETWORK`): fetch example AOI untiled and with a small `tile_size`; assert tiled is a per-year list of single-layer SpatRasters with `stac_items` attr, and **tiled == untiled over cropped common AOI cells** (not raw dimensions) +- [ ] `devtools::document()`; `lintr::lint_package()` clean; `devtools::test()` green + +## Phase 5: docs + gotchas note + NEWS + version +- [ ] `inst/notes/gdalcubes-pc-gotchas.md`: tiling entry (fetch download bounded by tiling the cube_view; tiled mosaic cached as `.tif` via terra; #36) +- [ ] `NEWS.md` `# drift 0.6.0` — `tile_size` (opt-in, default `NULL` = unchanged; bounds download for sparse AOIs; tiled fetches cache as `.tif`); Closes #36 +- [ ] `DESCRIPTION` `0.5.0 → 0.6.0` + `Date` (final commit) + +## Phase 6: validate, archive, PR, release +- [ ] `devtools::test()` / `lint` / `document` / `check` clean (network tests skip) +- [ ] `/planning-archive`; `/gh-pr-push` (`Fixes #36`, `Relates to NewGraphEnvironment/sred-2025-2026#16`) +- [ ] `/gh-pr-merge` → release v0.6.0 + +## Validation +- [ ] Tests pass (`devtools::test()`); network tests skip cleanly +- [ ] `/code-check` clean on each commit +- [ ] PWF checkboxes match landed work +- [ ] `/planning-archive` on completion From c4700a06b6e2c1739f844f92a94a49d559aca3df Mon Sep 17 00:00:00 2001 From: almac2022 Date: Thu, 9 Jul 2026 15:33:02 -0700 Subject: [PATCH 2/7] Add tile_grid() + tile_size_check() helpers for download tiling (#36) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Phase 1 of #36: the offline, network-free core of bounding dft_stac_fetch()'s download to the AOI footprint. tile_size_check() validates a tile_size and snaps it to a multiple of res (so tile pixel grids align for a seam-free merge); tile_grid() splits the AOI bbox into a res-aligned grid anchored at the lower- left and keeps only tiles intersecting the AOI polygon — a thin corridor fetches near its footprint, not its full bbox. Boundary tiles are left un-trimmed (the overhang is masked away; trimming would break res-alignment). Tests written first (validation, snapping, intersecting-subset, alignment, single-tile, degenerate-AOI); 40 pass / 1 skip; lint clean; code-check clean. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01DKqnW6Fni3tc44Lo5D5fJz --- R/dft_stac_fetch.R | 67 ++++++++++++++++++++++++ planning/active/progress.md | 6 ++- planning/active/task_plan.md | 10 ++-- tests/testthat/test-dft_stac_fetch.R | 78 ++++++++++++++++++++++++++++ 4 files changed, 155 insertions(+), 6 deletions(-) diff --git a/R/dft_stac_fetch.R b/R/dft_stac_fetch.R index 392f116..7d5fda4 100644 --- a/R/dft_stac_fetch.R +++ b/R/dft_stac_fetch.R @@ -179,6 +179,73 @@ stac_cache_key <- function(aoi_target, res, target_crs, dt, aggregation, } +#' Validate and snap a download `tile_size` to the pixel grid +#' +#' `tile_size` (CRS units) controls the download-tiling grid (#36). It is +#' snapped to a multiple of `res` so every tile's pixel grid aligns to the same +#' `res`-lattice — a prerequisite for a seam-free `terra::merge()` of the tiles. +#' Caller only invokes this for a non-NULL `tile_size`; `NULL` gates the whole +#' tiled path upstream. Returns the snapped size (a single positive numeric). +#' @noRd +tile_size_check <- function(tile_size, res) { + if (!is.numeric(tile_size) || length(tile_size) != 1L || + !is.finite(tile_size) || tile_size <= 0) { + cli::cli_abort(c( + "{.arg tile_size} must be a single positive finite number (CRS units) \\ + or {.code NULL}.", + "x" = "Got {.obj_type_friendly {tile_size}}." + )) + } + snapped <- round(tile_size / res) * res + if (snapped < res) { + cli::cli_abort(c( + "{.arg tile_size} ({tile_size}) snaps to {snapped}, smaller than \\ + {.arg res} ({res}).", + "i" = "Choose a {.arg tile_size} at least as large as {.arg res}." + )) + } + if (!isTRUE(all.equal(snapped, tile_size))) { + cli::cli_inform( + "{.arg tile_size} snapped from {tile_size} to {snapped} \\ + (a multiple of {.arg res} = {res})." + ) + } + snapped +} + + +#' Build the res-aligned download tiles that intersect the AOI +#' +#' Splits the AOI bounding box into a grid of `tile_size`-square cells anchored +#' at the bbox lower-left (the same origin as the single-cube extent), and keeps +#' only cells that intersect the AOI polygon — so a thin corridor fetches near +#' its footprint, not its full bbox (#36). Boundary cells are left un-trimmed +#' past the bbox: trimming the max edge would break `res`-alignment, and the +#' `< tile_size` overhang is dropped by the final `terra::mask()` anyway. +#' `tile_size` must already be snapped (see [tile_size_check()]). +#' @return A list of `list(left, right, bottom, top)` extents for [gdalcubes::cube_view()]. +#' @noRd +tile_grid <- function(aoi_target, tile_size, res) { + bbox <- sf::st_bbox(aoi_target) + grid <- sf::st_make_grid( + sf::st_as_sfc(bbox), + cellsize = tile_size, + offset = c(bbox[["xmin"]], bbox[["ymin"]]) + ) + aoi_union <- sf::st_union(sf::st_geometry(aoi_target)) + grid <- grid[lengths(sf::st_intersects(grid, aoi_union)) > 0] + if (length(grid) == 0) { + cli::cli_abort("No download tiles intersect the AOI \\ + (is the AOI geometry valid and non-empty?).") + } + lapply(grid, function(cell) { + b <- sf::st_bbox(cell) + list(left = b[["xmin"]], right = b[["xmax"]], + bottom = b[["ymin"]], top = b[["ymax"]]) + }) +} + + #' Auto-detect UTM EPSG code from sf geometry #' @noRd auto_utm_epsg <- function(x) { diff --git a/planning/active/progress.md b/planning/active/progress.md index 2675db8..b46b4a3 100644 --- a/planning/active/progress.md +++ b/planning/active/progress.md @@ -10,4 +10,8 @@ regression for `tile_size = NULL`; GDAL `/vsicurl` config on the tiled path. - Created branch `36-tile-dft-stac-fetch-to-bound-download-over-spars` off main. - Scaffolded PWF baseline with approved phases. -- Next: start Phase 1 (tests-first for `tile_grid()` + `tile_size` normalization). +- Phase 1 done: `tile_size_check()` (validate + snap to res multiple) and + `tile_grid()` (res-aligned tiles intersecting the AOI) added to + `dft_stac_fetch.R`, with 6 offline test blocks written first (confirmed red, + then green — 40 pass / 1 skip). Lint clean. +- Next: Phase 2 (cache-key conditional `tile_size` append + golden regression). diff --git a/planning/active/task_plan.md b/planning/active/task_plan.md index dd3ea54..d1fec93 100644 --- a/planning/active/task_plan.md +++ b/planning/active/task_plan.md @@ -12,11 +12,11 @@ single-cube behavior, existing caches preserved. Sibling cube-path residual is tracked separately as #38 (out of scope). ## Phase 1: `tile_grid()` + `tile_size` normalization (offline, tests first) -- [ ] validation aborts: `NA`, `0`, negative, `Inf`, `c(1,2)`, `"500"`, snapped `< res` -- [ ] snapping: non-multiple snaps to nearest multiple of `res`, messages, snapped value used -- [ ] intersecting subset: a diagonal/L-shaped AOI in a known bbox keeps exactly the expected tiles, drops empty ones, kept-count ≪ full grid (the efficiency mechanism, offline) -- [ ] res-alignment: every extent's `left-xmin`, width, height are exact multiples of `res`; first tile lower-left `== (xmin, ymin)`; single-tile when `tile_size ≥ bbox`; empty/degenerate AOI aborts -- [ ] implement `tile_grid()` + normalization; new tests green +- [x] validation aborts: `NA`, `0`, negative, `Inf`, `c(1,2)`, `"500"`, snapped `< res` +- [x] snapping: non-multiple snaps to nearest multiple of `res`, messages, snapped value used +- [x] intersecting subset: a diagonal/L-shaped AOI in a known bbox keeps exactly the expected tiles, drops empty ones, kept-count ≪ full grid (the efficiency mechanism, offline) +- [x] res-alignment: every extent's `left-xmin`, width, height are exact multiples of `res`; first tile lower-left `== (xmin, ymin)`; single-tile when `tile_size ≥ bbox`; empty/degenerate AOI aborts +- [x] implement `tile_grid()` + normalization; new tests green ## Phase 2: cache key — conditional `tile_size` append (offline) - [ ] golden regression: `cache_key(tile_size = NULL)` == frozen current 12-char hash for fixed inputs (guards legacy-cache preservation) diff --git a/tests/testthat/test-dft_stac_fetch.R b/tests/testthat/test-dft_stac_fetch.R index f0fbe2b..d39cec6 100644 --- a/tests/testthat/test-dft_stac_fetch.R +++ b/tests/testthat/test-dft_stac_fetch.R @@ -80,3 +80,81 @@ test_that("stac_cache_key ignores sf attribute columns", { with_attrs <- sf::st_sf(name = "a", area = 1.5, geometry = bare) expect_equal(cache_key(with_attrs), cache_key(bare)) }) + +# --- tile_size_check(): validate + snap tile_size to a multiple of res ------- +# Offline; the download-tiling normalization (#36). NULL is handled by the +# caller (it gates the whole tiled path); this helper only sees non-NULL input. +test_that("tile_size_check aborts on non-positive / non-finite / non-scalar input", { + expect_error(drift:::tile_size_check(NA, 10), "positive") + expect_error(drift:::tile_size_check(0, 10), "positive") + expect_error(drift:::tile_size_check(-5, 10), "positive") + expect_error(drift:::tile_size_check(Inf, 10), "positive") + expect_error(drift:::tile_size_check(c(1, 2), 10), "positive") + expect_error(drift:::tile_size_check("500", 10), "positive") +}) + +test_that("tile_size_check aborts when the snapped size is smaller than res", { + # 4 snaps to round(4/10)*10 = 0, which is < res + expect_error(drift:::tile_size_check(4, 10), "res") +}) + +test_that("tile_size_check snaps to the nearest multiple of res and returns it", { + expect_equal(drift:::tile_size_check(500, 10), 500) # already aligned + expect_equal(drift:::tile_size_check(504, 10), 500) # rounds down + expect_equal(drift:::tile_size_check(506, 10), 510) # rounds up + expect_message(drift:::tile_size_check(504, 10), "snap") +}) + +# --- tile_grid(): res-aligned tiles intersecting the AOI (offline) ----------- +# A rectangular AOI filling a bbox (all tiles kept) and a thin diagonal corridor +# (most bbox tiles dropped — the download-saving mechanism, tested without network). +rect_aoi <- function(xmin = 0, ymin = 0, xmax = 1000, ymax = 1000, crs = 32609) { + sf::st_sfc( + sf::st_polygon(list(rbind( + c(xmin, ymin), c(xmax, ymin), c(xmax, ymax), c(xmin, ymax), c(xmin, ymin) + ))), + crs = crs + ) +} + +test_that("tile_grid returns res-aligned tiles anchored at (xmin, ymin)", { + aoi <- rect_aoi(0, 0, 1000, 1000) # 2x2 tiles at tile_size 500 + tiles <- drift:::tile_grid(aoi, tile_size = 500, res = 10) + expect_length(tiles, 4) + lefts <- vapply(tiles, `[[`, numeric(1), "left") + bottoms <- vapply(tiles, `[[`, numeric(1), "bottom") + widths <- vapply(tiles, function(t) t$right - t$left, numeric(1)) + heights <- vapply(tiles, function(t) t$top - t$bottom, numeric(1)) + # anchored at (0, 0): every left/bottom is a multiple of tile_size from origin + expect_setequal(lefts, c(0, 500)) + expect_setequal(bottoms, c(0, 500)) + # every tile is tile_size (a multiple of res) wide and tall + expect_true(all(widths == 500)) + expect_true(all(heights == 500)) + # each edge lands on the res-lattice anchored at the bbox lower-left + expect_true(all(lefts %% 10 == 0)) + expect_true(all(bottoms %% 10 == 0)) +}) + +test_that("tile_grid drops bbox tiles that miss the AOI (diagonal corridor)", { + line <- sf::st_sfc(sf::st_linestring(rbind(c(0, 0), c(1000, 1000))), crs = 32609) + aoi <- sf::st_buffer(line, 20) # thin diagonal corridor + tiles <- drift:::tile_grid(aoi, tile_size = 500, res = 10) + # full grid over the buffered bbox is 3x3 = 9; the diagonal keeps a strict subset + expect_gt(length(tiles), 0) + expect_lt(length(tiles), 9) +}) + +test_that("tile_grid yields a single tile when tile_size covers the bbox", { + tiles <- drift:::tile_grid(rect_aoi(0, 0, 400, 400), tile_size = 500, res = 10) + expect_length(tiles, 1) + expect_equal(tiles[[1]]$left, 0) + expect_equal(tiles[[1]]$bottom, 0) +}) + +test_that("tile_grid errors on a degenerate (empty) AOI", { + expect_error( + drift:::tile_grid(sf::st_sfc(sf::st_polygon(), crs = 32609), + tile_size = 500, res = 10) + ) +}) From cf04bfe217e2a5816aaa5c194e35612ac3a3150d Mon Sep 17 00:00:00 2001 From: almac2022 Date: Thu, 9 Jul 2026 15:36:47 -0700 Subject: [PATCH 3/7] Key tiled fetches distinctly, preserving legacy untiled cache keys (#36) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Phase 2 of #36: stac_cache_key() gains an optional tile_size, appended to the hash only when non-NULL. An untiled fetch (tile_size = NULL) reproduces the exact pre-tiling 9-element hash, so existing io-lulc caches stay valid on upgrade; a tiled fetch keys distinctly, so its terra .tif mosaic is never served as an untiled gdalcubes .nc (or vice versa). A frozen golden-hash regression guards the legacy-preservation invariant. tile_size arrives already snapped by the caller (single normalization site — no gate/key desync). 45 pass / 1 skip; lint + code-check clean. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01DKqnW6Fni3tc44Lo5D5fJz --- R/dft_stac_fetch.R | 21 ++++++++++++-------- planning/active/progress.md | 9 ++++++++- planning/active/task_plan.md | 6 +++--- tests/testthat/test-dft_stac_fetch.R | 29 ++++++++++++++++++++++++++-- 4 files changed, 51 insertions(+), 14 deletions(-) diff --git a/R/dft_stac_fetch.R b/R/dft_stac_fetch.R index 7d5fda4..ad85952 100644 --- a/R/dft_stac_fetch.R +++ b/R/dft_stac_fetch.R @@ -164,18 +164,23 @@ dft_stac_fetch <- function(aoi, #' representation differences can't change the key; the CRS enters separately #' as `target_crs`. `res` is coerced to double so `10L` and `10` key alike. #' Callers must pass post-resolution `stac_url`/`collection`/`asset`, never -#' the raw possibly-NULL arguments. +#' the raw possibly-NULL arguments. `tile_size` (the download-tiling grid, #36) +#' is appended to the hash ONLY when non-NULL, so an untiled fetch keeps the +#' exact legacy 9-element hash (existing caches stay valid) while a tiled fetch +#' keys distinctly. It must arrive already snapped by the caller. #' @noRd stac_cache_key <- function(aoi_target, res, target_crs, dt, aggregation, - resampling, stac_url, collection, asset) { + resampling, stac_url, collection, asset, + tile_size = NULL) { 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, asset - )), - 1, 12 + parts <- list( + geom_wkb, as.numeric(res), target_crs, dt, aggregation, + resampling, stac_url, collection, asset ) + # A tiled fetch caches a terra .tif mosaic; an untiled fetch caches a + # gdalcubes .nc. Keying them apart stops one being served as the other. + if (!is.null(tile_size)) parts <- c(parts, list(as.numeric(tile_size))) + substr(rlang::hash(parts), 1, 12) } diff --git a/planning/active/progress.md b/planning/active/progress.md index b46b4a3..27fc15b 100644 --- a/planning/active/progress.md +++ b/planning/active/progress.md @@ -14,4 +14,11 @@ `tile_grid()` (res-aligned tiles intersecting the AOI) added to `dft_stac_fetch.R`, with 6 offline test blocks written first (confirmed red, then green — 40 pass / 1 skip). Lint clean. -- Next: Phase 2 (cache-key conditional `tile_size` append + golden regression). +- Phase 2 done: `stac_cache_key()` gains optional `tile_size`, appended to the + hash only when non-NULL — so an untiled fetch reproduces the exact legacy hash + (frozen `79f67b7b9dae` golden test guards this) and a tiled fetch keys + distinctly. Test helper snaps via `tile_size_check` to mirror production + (504≡500). 45 pass / 1 skip; lint + code-check clean. Call-site wiring + deferred to Phase 4 (arrives with the fetch param). +- Next: Phase 3 (extract `fetch_extent_to()`; refactor untiled path, no + behavior change). diff --git a/planning/active/task_plan.md b/planning/active/task_plan.md index d1fec93..db7e635 100644 --- a/planning/active/task_plan.md +++ b/planning/active/task_plan.md @@ -19,9 +19,9 @@ tracked separately as #38 (out of scope). - [x] implement `tile_grid()` + normalization; new tests green ## Phase 2: cache key — conditional `tile_size` append (offline) -- [ ] golden regression: `cache_key(tile_size = NULL)` == frozen current 12-char hash for fixed inputs (guards legacy-cache preservation) -- [ ] `cache_key(tile_size = 500) != cache_key(tile_size = NULL)`; distinct sizes → distinct keys; snap-before-key: `504` and `500` (res 10) → same key -- [ ] `stac_cache_key()` gains `tile_size = NULL`; append only when non-NULL; call site passes `tile_size`; tests green +- [x] golden regression: `cache_key(tile_size = NULL)` == frozen `79f67b7b9dae` for fixed inputs (guards legacy-cache preservation) +- [x] `cache_key(tile_size = 500) != cache_key(tile_size = NULL)`; distinct sizes → distinct keys; snap-before-key: `504` and `500` (res 10) → same key +- [x] `stac_cache_key()` gains `tile_size = NULL`; appends only when non-NULL (real call-site wiring lands in Phase 4 with the fetch param); tests green ## Phase 3: extract `fetch_extent_to()`, refactor untiled path (behavior-preserving) - [ ] extract helper; untiled path routes through it, writing straight to `_.nc` — identical filename/format diff --git a/tests/testthat/test-dft_stac_fetch.R b/tests/testthat/test-dft_stac_fetch.R index d39cec6..c835398 100644 --- a/tests/testthat/test-dft_stac_fetch.R +++ b/tests/testthat/test-dft_stac_fetch.R @@ -43,9 +43,14 @@ square_aoi <- function(dx = 0) { cache_key <- function(aoi = square_aoi(), res = 10, target_crs = "EPSG:32609", dt = "P1Y", aggregation = "first", resampling = "near", stac_url = "https://example.com/stac", - collection = "test-collection", asset = "data") { + collection = "test-collection", asset = "data", + tile_size = NULL) { + # mirror production: dft_stac_fetch() snaps tile_size once (tile_size_check) + # before it reaches both the tile grid and the cache key + ts <- if (is.null(tile_size)) NULL else + suppressMessages(drift:::tile_size_check(tile_size, res)) drift:::stac_cache_key(aoi, res, target_crs, dt, aggregation, resampling, - stac_url, collection, asset) + stac_url, collection, asset, tile_size = ts) } test_that("stac_cache_key is deterministic and 12-char hex", { @@ -75,6 +80,26 @@ test_that("stac_cache_key treats integer and double res alike", { expect_equal(cache_key(res = 10L), cache_key(res = 10)) }) +test_that("stac_cache_key(tile_size = NULL) reproduces the legacy pre-tiling hash", { + # Frozen guardian of legacy-cache preservation (#36): adding tile_size must + # NOT change the key for an untiled fetch, or every cached io-lulc fetch + # silently re-downloads on upgrade. If this literal must change, that is a + # deliberate cache-format break — flag it, don't just re-freeze. + expect_equal(cache_key(), "79f67b7b9dae") +}) + +test_that("stac_cache_key keys a tiled fetch distinctly from an untiled one", { + base <- cache_key() + expect_false(cache_key(tile_size = 500) == base) + expect_false(cache_key(tile_size = 1000) == base) + expect_false(cache_key(tile_size = 500) == cache_key(tile_size = 1000)) +}) + +test_that("stac_cache_key snaps tile_size before hashing (aligned sizes key alike)", { + # 504 and 500 both snap to 500 (res 10), so they must hit the same cache + expect_equal(cache_key(tile_size = 504), cache_key(tile_size = 500)) +}) + test_that("stac_cache_key ignores sf attribute columns", { bare <- square_aoi() with_attrs <- sf::st_sf(name = "a", area = 1.5, geometry = bare) From b7bdd592180316670d86f0cb4d165b516a7ae61f Mon Sep 17 00:00:00 2001 From: almac2022 Date: Thu, 9 Jul 2026 15:38:26 -0700 Subject: [PATCH 4/7] Extract fetch_extent_to() shared by untiled and (upcoming) tiled fetch (#36) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Phase 3 of #36: pull the cube_view + raster_cube + write_ncdf block out of the per-year loop into fetch_extent_to(). The untiled path now routes through it, writing straight to the existing _.nc cache — identical filename, format, and I/O to before (faithful code-motion; cube_view built identically). Sharing this primitive is what will guarantee each tile fetches identically to the corresponding slice of the untiled cube when the tiled branch lands in Phase 4. 45 pass / 1 skip; lint clean. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01DKqnW6Fni3tc44Lo5D5fJz --- R/dft_stac_fetch.R | 48 ++++++++++++++++++++++++------------ planning/active/progress.md | 8 ++++-- planning/active/task_plan.md | 4 +-- 3 files changed, 40 insertions(+), 20 deletions(-) diff --git a/R/dft_stac_fetch.R b/R/dft_stac_fetch.R index ad85952..9589e4f 100644 --- a/R/dft_stac_fetch.R +++ b/R/dft_stac_fetch.R @@ -127,23 +127,12 @@ dft_stac_fetch <- function(aoi, r <- terra::rast(cache_file) } else { message(" ", yr, ": fetching...") - v <- gdalcubes::cube_view( - srs = target_crs, - extent = list( - left = bbox_target["xmin"], - right = bbox_target["xmax"], - bottom = bbox_target["ymin"], - top = bbox_target["ymax"], - t0 = paste0(yr, "-01-01"), - t1 = paste0(yr, "-12-31") - ), - dx = res, dy = res, - dt = dt, - aggregation = aggregation, - resampling = resampling + ext <- list( + left = bbox_target[["xmin"]], right = bbox_target[["xmax"]], + bottom = bbox_target[["ymin"]], top = bbox_target[["ymax"]] ) - cube <- gdalcubes::raster_cube(col, v) - gdalcubes::write_ncdf(cube, cache_file, overwrite = TRUE) + fetch_extent_to(col, ext, paste0(yr, "-01-01"), paste0(yr, "-12-31"), + target_crs, res, dt, aggregation, resampling, cache_file) r <- terra::rast(cache_file) } @@ -251,6 +240,33 @@ tile_grid <- function(aoi_target, tile_size, res) { } +#' Fetch one gdalcubes cube over a single space+time extent to a NetCDF file +#' +#' The `cube_view` + `raster_cube` + `write_ncdf` block shared by the untiled +#' fetch (one call over the AOI bbox) and the tiled fetch (one call per tile, +#' #36). Sharing this primitive is what guarantees a tile fetches identically to +#' the corresponding slice of the untiled cube. `ext` is a list with +#' `left`/`right`/`bottom`/`top`; `t0`/`t1` bound the year. Writes to `out_nc` +#' and returns it (the caller reads it back with [terra::rast()]). +#' @noRd +fetch_extent_to <- function(col, ext, t0, t1, target_crs, res, dt, + aggregation, resampling, out_nc) { + v <- gdalcubes::cube_view( + srs = target_crs, + extent = list( + left = ext$left, right = ext$right, + bottom = ext$bottom, top = ext$top, + t0 = t0, t1 = t1 + ), + dx = res, dy = res, dt = dt, + aggregation = aggregation, resampling = resampling + ) + cube <- gdalcubes::raster_cube(col, v) + gdalcubes::write_ncdf(cube, out_nc, overwrite = TRUE) + out_nc +} + + #' Auto-detect UTM EPSG code from sf geometry #' @noRd auto_utm_epsg <- function(x) { diff --git a/planning/active/progress.md b/planning/active/progress.md index 27fc15b..4820719 100644 --- a/planning/active/progress.md +++ b/planning/active/progress.md @@ -20,5 +20,9 @@ distinctly. Test helper snaps via `tile_size_check` to mirror production (504≡500). 45 pass / 1 skip; lint + code-check clean. Call-site wiring deferred to Phase 4 (arrives with the fetch param). -- Next: Phase 3 (extract `fetch_extent_to()`; refactor untiled path, no - behavior change). +- Phase 3 done: extracted `fetch_extent_to()` (cube_view + raster_cube + + write_ncdf); untiled path routes through it writing straight to the existing + `_.nc` — faithful code-motion, cube_view built identically, 45 pass / + 1 skip, lint clean. Shared primitive de-risks Phase 4 (tiles reuse it). +- Next: Phase 4 (tiled branch — per-tile fetch → terra::merge → .tif → mask; + wire `tile_size` end-to-end; offline merge oracle + opt-in network e2e). diff --git a/planning/active/task_plan.md b/planning/active/task_plan.md index db7e635..d43f296 100644 --- a/planning/active/task_plan.md +++ b/planning/active/task_plan.md @@ -24,8 +24,8 @@ tracked separately as #38 (out of scope). - [x] `stac_cache_key()` gains `tile_size = NULL`; appends only when non-NULL (real call-site wiring lands in Phase 4 with the fetch param); tests green ## Phase 3: extract `fetch_extent_to()`, refactor untiled path (behavior-preserving) -- [ ] extract helper; untiled path routes through it, writing straight to `_.nc` — identical filename/format -- [ ] existing tests + opt-in untiled network test still green (no observable change) +- [x] extract helper; untiled path routes through it, writing straight to `_.nc` — identical filename/format +- [x] existing offline tests green (45/1 skip); faithful code-motion (cube_view built identically); the Phase 4 network e2e exercises the shared primitive end-to-end ## Phase 4: tiled branch — mosaic assembly + wire `tile_size` end-to-end - [ ] offline merge oracle: reference SpatRaster split into res-lattice tiles → `terra::merge()` → `all.equal(values(merged), values(reference))`; masked mosaic == masked reference over AOI; `.tif` round-trip preserves single-layer integer codes From 6dd67c9f44cc936932051c6ed636a2aed591e0f7 Mon Sep 17 00:00:00 2001 From: almac2022 Date: Thu, 9 Jul 2026 15:50:21 -0700 Subject: [PATCH 5/7] Add opt-in tile_size to bound dft_stac_fetch download to the AOI (#36) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Phase 4 of #36: tiled fetch end-to-end. When tile_size is set, dft_stac_fetch() splits the AOI bbox into a res-aligned grid, streams only tiles that intersect the AOI polygon (via the shared fetch_extent_to()), mosaics them with terra::merge(terra::sprc(...)) (new mosaic_tiles() helper), and caches a terra GeoTIFF — so a thin corridor fetches near its footprint instead of its full bounding box. Extension routes on is.null(tile_size) (.nc untiled / .tif tiled), matching the distinct cache key. tile_size is normalized once (snap to res multiple) so the path gate and cache key agree; GDAL /vsicurl config (mirroring dft_stac_cube) is set with an on.exit restore, scoped to the tiled path, since tiling multiplies per-item COG opens. Default tile_size = NULL is unchanged. Offline merge oracle (byte-for-byte reassembly, single-tile, .tif round-trip) and an opt-in DRIFT_TEST_NETWORK e2e (tiled == untiled over the in-AOI overlap, extension routing). 352 pass / 5 skip; document + lint + code-check clean. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01DKqnW6Fni3tc44Lo5D5fJz --- R/dft_stac_fetch.R | 102 ++++++++++++++++++++++----- man/dft_stac_fetch.Rd | 22 ++++-- planning/active/progress.md | 12 +++- planning/active/task_plan.md | 10 +-- tests/testthat/test-dft_stac_fetch.R | 84 ++++++++++++++++++++++ 5 files changed, 200 insertions(+), 30 deletions(-) diff --git a/R/dft_stac_fetch.R b/R/dft_stac_fetch.R index 9589e4f..e5f0424 100644 --- a/R/dft_stac_fetch.R +++ b/R/dft_stac_fetch.R @@ -6,11 +6,11 @@ #' custom COGs). #' #' Fetched rasters are cached under [dft_cache_path()] as -#' `/_.nc`, where `key` is a hash of the AOI geometry and -#' every fetch parameter that affects the output (`res`, `crs`, `dt`, -#' `aggregation`, `resampling`, `stac_url`, `collection`, `asset`). Repeat -#' calls with the same AOI and parameters reuse the cache; changing any of -#' them re-fetches. +#' `/_.nc` (or `.tif` when `tile_size` is set — see below), +#' where `key` is a hash of the AOI geometry and every fetch parameter that +#' affects the output (`res`, `crs`, `dt`, `aggregation`, `resampling`, +#' `stac_url`, `collection`, `asset`, and `tile_size`). Repeat calls with the +#' same AOI and parameters reuse the cache; changing any of them re-fetches. #' #' @param aoi An `sf` polygon defining the area of interest. #' @param source Character. A known source name passed to [dft_stac_config()]. @@ -29,6 +29,16 @@ #' `"first"`). Use `"median"` for multi-scene composites. #' @param resampling Character. Spatial resampling method (default `"near"` #' for categorical data). +#' @param tile_size Numeric or `NULL` (default). Edge length, in CRS units +#' (metres for the default UTM CRS), of the download-tiling grid. When `NULL`, +#' one cube is streamed over the whole AOI bounding box (the download scales +#' with the bbox, not the AOI). When set, the bbox is split into a grid of +#' `tile_size`-square tiles and only tiles that intersect the AOI polygon are +#' streamed, then mosaicked — so a thin, diagonal AOI (e.g. a floodplain +#' corridor) fetches close to its footprint instead of its full bounding box. +#' Snapped to a multiple of `res`. Smaller tiles waste less bbox but cost more +#' per-tile round trips; there is no auto-tuning. Tiled fetches cache a terra +#' GeoTIFF (`.tif`) rather than a gdalcubes NetCDF (`.nc`). #' @param cache_dir Character. Cache directory path. When `NULL`, uses #' [dft_cache_path()]. #' @param force Logical. Re-fetch even if cached, overwriting the cached file @@ -53,11 +63,34 @@ dft_stac_fetch <- function(aoi, dt = "P1Y", aggregation = "first", resampling = "near", + tile_size = NULL, cache_dir = NULL, force = FALSE, sign_fn = rstac::sign_planetary_computer()) { rlang::check_installed("gdalcubes", reason = "to fetch STAC rasters") + # Normalize tile_size ONCE so the path gate (is.null) and the cache key derive + # from the same snapped scalar. When tiling, tune GDAL for the many extra + # per-item COG opens (restored on exit so the caller's session is untouched). + if (!is.null(tile_size)) { + tile_size <- tile_size_check(tile_size, res) + gdal_cfg <- c( + GDAL_DISABLE_READDIR_ON_OPEN = "EMPTY_DIR", + GDAL_HTTP_MULTIPLEX = "YES", + GDAL_HTTP_VERSION = "2", + VSI_CACHE = "TRUE", + CPL_VSIL_CURL_ALLOWED_EXTENSIONS = ".tif" + ) + old_cfg <- Sys.getenv(names(gdal_cfg), unset = NA) + do.call(Sys.setenv, as.list(gdal_cfg)) + on.exit({ + set_again <- old_cfg[!is.na(old_cfg)] + if (length(set_again)) do.call(Sys.setenv, as.list(set_again)) + unset <- names(old_cfg)[is.na(old_cfg)] + if (length(unset)) Sys.unsetenv(unset) + }, add = TRUE) + } + # Resolve config if (is.null(stac_url) || is.null(collection) || is.null(asset)) { cfg <- dft_stac_config(source) @@ -115,28 +148,44 @@ dft_stac_fetch <- function(aoi, dir.create(cache_source_dir, recursive = TRUE, showWarnings = FALSE) cache_key <- stac_cache_key( aoi_target, res, target_crs, dt, aggregation, resampling, - stac_url, collection, asset + stac_url, collection, asset, tile_size = tile_size ) + # A tiled fetch mosaics per-tile cubes with terra and caches a GeoTIFF; an + # untiled fetch writes a single gdalcubes NetCDF. The grid and full-bbox extent + # are constant across years, so build them once. + ext_out <- if (is.null(tile_size)) "nc" else "tif" + bbox_ext <- list( + left = bbox_target[["xmin"]], right = bbox_target[["xmax"]], + bottom = bbox_target[["ymin"]], top = bbox_target[["ymax"]] + ) + tiles <- if (is.null(tile_size)) NULL else tile_grid(aoi_target, tile_size, res) + # Fetch per year result <- lapply(years, function(yr) { - cache_file <- file.path(cache_source_dir, paste0(yr, "_", cache_key, ".nc")) + cache_file <- file.path(cache_source_dir, + paste0(yr, "_", cache_key, ".", ext_out)) + t0 <- paste0(yr, "-01-01") + t1 <- paste0(yr, "-12-31") if (!force && file.exists(cache_file)) { message(" ", yr, ": cached") - r <- terra::rast(cache_file) - } else { + } else if (is.null(tile_size)) { message(" ", yr, ": fetching...") - ext <- list( - left = bbox_target[["xmin"]], right = bbox_target[["xmax"]], - bottom = bbox_target[["ymin"]], top = bbox_target[["ymax"]] - ) - fetch_extent_to(col, ext, paste0(yr, "-01-01"), paste0(yr, "-12-31"), - target_crs, res, dt, aggregation, resampling, cache_file) - r <- terra::rast(cache_file) + fetch_extent_to(col, bbox_ext, t0, t1, target_crs, res, dt, + aggregation, resampling, cache_file) + } else { + message(" ", yr, ": fetching ", length(tiles), " tile(s)...") + tile_files <- vapply(seq_along(tiles), function(i) { + fetch_extent_to(col, tiles[[i]], t0, t1, target_crs, res, dt, + aggregation, resampling, + tempfile(sprintf("drift_tile%d_", i), fileext = ".nc")) + }, character(1)) + mosaic_tiles(tile_files, cache_file) + unlink(tile_files) } - terra::mask(r, terra::vect(aoi_target)) + terra::mask(terra::rast(cache_file), terra::vect(aoi_target)) }) names(result) <- as.character(years) @@ -216,7 +265,7 @@ tile_size_check <- function(tile_size, res) { #' its footprint, not its full bbox (#36). Boundary cells are left un-trimmed #' past the bbox: trimming the max edge would break `res`-alignment, and the #' `< tile_size` overhang is dropped by the final `terra::mask()` anyway. -#' `tile_size` must already be snapped (see [tile_size_check()]). +#' `tile_size` must already be snapped (see `tile_size_check()`). #' @return A list of `list(left, right, bottom, top)` extents for [gdalcubes::cube_view()]. #' @noRd tile_grid <- function(aoi_target, tile_size, res) { @@ -267,6 +316,23 @@ fetch_extent_to <- function(col, ext, t0, t1, target_crs, res, dt, } +#' Mosaic per-tile fetch outputs into one cache raster +#' +#' Reads each per-tile NetCDF, merges them (tiles are res-aligned and +#' non-overlapping, so `terra::merge()` reassembles without resampling), and +#' writes the mosaic to `out_file`. Written with `terra::writeRaster()` to a +#' GeoTIFF — not `gdalcubes::write_ncdf()` — because terra's own NetCDF *write* +#' is fragile on the pinned stack (see inst/notes/gdalcubes-pc-gotchas.md); this +#' mirrors the `dft_stac_cube()` cache. Returns `out_file`. +#' @noRd +mosaic_tiles <- function(tile_files, out_file) { + rasters <- lapply(tile_files, terra::rast) + merged <- terra::merge(terra::sprc(rasters)) + terra::writeRaster(merged, out_file, overwrite = TRUE) + out_file +} + + #' Auto-detect UTM EPSG code from sf geometry #' @noRd auto_utm_epsg <- function(x) { diff --git a/man/dft_stac_fetch.Rd b/man/dft_stac_fetch.Rd index 372f9b7..36391e1 100644 --- a/man/dft_stac_fetch.Rd +++ b/man/dft_stac_fetch.Rd @@ -16,6 +16,7 @@ dft_stac_fetch( dt = "P1Y", aggregation = "first", resampling = "near", + tile_size = NULL, cache_dir = NULL, force = FALSE, sign_fn = rstac::sign_planetary_computer() @@ -50,6 +51,17 @@ When \code{NULL}, auto-detected from the AOI centroid's UTM zone.} \item{resampling}{Character. Spatial resampling method (default \code{"near"} for categorical data).} +\item{tile_size}{Numeric or \code{NULL} (default). Edge length, in CRS units +(metres for the default UTM CRS), of the download-tiling grid. When \code{NULL}, +one cube is streamed over the whole AOI bounding box (the download scales +with the bbox, not the AOI). When set, the bbox is split into a grid of +\code{tile_size}-square tiles and only tiles that intersect the AOI polygon are +streamed, then mosaicked — so a thin, diagonal AOI (e.g. a floodplain +corridor) fetches close to its footprint instead of its full bounding box. +Snapped to a multiple of \code{res}. Smaller tiles waste less bbox but cost more +per-tile round trips; there is no auto-tuning. Tiled fetches cache a terra +GeoTIFF (\code{.tif}) rather than a gdalcubes NetCDF (\code{.nc}).} + \item{cache_dir}{Character. Cache directory path. When \code{NULL}, uses \code{\link[=dft_cache_path]{dft_cache_path()}}.} @@ -74,9 +86,9 @@ custom COGs). } \details{ Fetched rasters are cached under \code{\link[=dft_cache_path]{dft_cache_path()}} as -\verb{/_.nc}, where \code{key} is a hash of the AOI geometry and -every fetch parameter that affects the output (\code{res}, \code{crs}, \code{dt}, -\code{aggregation}, \code{resampling}, \code{stac_url}, \code{collection}, \code{asset}). Repeat -calls with the same AOI and parameters reuse the cache; changing any of -them re-fetches. +\verb{/_.nc} (or \code{.tif} when \code{tile_size} is set — see below), +where \code{key} is a hash of the AOI geometry and every fetch parameter that +affects the output (\code{res}, \code{crs}, \code{dt}, \code{aggregation}, \code{resampling}, +\code{stac_url}, \code{collection}, \code{asset}, and \code{tile_size}). Repeat calls with the +same AOI and parameters reuse the cache; changing any of them re-fetches. } diff --git a/planning/active/progress.md b/planning/active/progress.md index 4820719..86017ff 100644 --- a/planning/active/progress.md +++ b/planning/active/progress.md @@ -24,5 +24,13 @@ write_ncdf); untiled path routes through it writing straight to the existing `_.nc` — faithful code-motion, cube_view built identically, 45 pass / 1 skip, lint clean. Shared primitive de-risks Phase 4 (tiles reuse it). -- Next: Phase 4 (tiled branch — per-tile fetch → terra::merge → .tif → mask; - wire `tile_size` end-to-end; offline merge oracle + opt-in network e2e). +- Phase 4 done: `mosaic_tiles()` helper + `tile_size` wired end-to-end. Tiled + branch fetches each intersecting tile via the shared `fetch_extent_to()`, + `terra::merge(terra::sprc(...))`, writes a `.tif` mosaic (extension routed by + `is.null(tile_size)`), unlinks temps, then masks. GDAL `/vsicurl` config with + on.exit restore, scoped to the tiled path. Offline merge oracle (byte-for-byte + reassembly + single-tile + .tif round-trip) + opt-in network e2e (tiled == + untiled, resampled onto the untiled grid). Roxygen `@param tile_size` + cache + doc. 352 pass / 5 skip; document + lint + code-check clean. Code-check flagged + the network test's sub-pixel-offset fragility → hardened via near-resample. +- Next: Phase 5 (gotchas note + NEWS 0.6.0 + DESCRIPTION bump). diff --git a/planning/active/task_plan.md b/planning/active/task_plan.md index d43f296..7a59153 100644 --- a/planning/active/task_plan.md +++ b/planning/active/task_plan.md @@ -28,11 +28,11 @@ tracked separately as #38 (out of scope). - [x] existing offline tests green (45/1 skip); faithful code-motion (cube_view built identically); the Phase 4 network e2e exercises the shared primitive end-to-end ## Phase 4: tiled branch — mosaic assembly + wire `tile_size` end-to-end -- [ ] offline merge oracle: reference SpatRaster split into res-lattice tiles → `terra::merge()` → `all.equal(values(merged), values(reference))`; masked mosaic == masked reference over AOI; `.tif` round-trip preserves single-layer integer codes -- [ ] add `tile_size` to `dft_stac_fetch`; tiled branch (per-tile fetch → merge → `.tif` → read → mask); extension from `is.null(tile_size)`; GDAL config + `on.exit`; tempfile `unlink` -- [ ] roxygen `@param tile_size` + `\dontrun` example + cache-doc `.tif` note -- [ ] opt-in network e2e (`DRIFT_TEST_NETWORK`): fetch example AOI untiled and with a small `tile_size`; assert tiled is a per-year list of single-layer SpatRasters with `stac_items` attr, and **tiled == untiled over cropped common AOI cells** (not raw dimensions) -- [ ] `devtools::document()`; `lintr::lint_package()` clean; `devtools::test()` green +- [x] offline merge oracle: `mosaic_tiles()` reassembles res-lattice tiles into the reference grid byte-for-byte (`terra::merge(terra::sprc(...))`); single-tile case; `.tif` round-trip preserves single-layer integer codes +- [x] add `tile_size` to `dft_stac_fetch`; tiled branch (per-tile fetch → merge → `.tif` → read → mask); extension from `is.null(tile_size)`; GDAL config + `on.exit`; tempfile `unlink` +- [x] roxygen `@param tile_size` + cache-doc `.tif` note (no `@examples` — neither `dft_stac_fetch` nor the cube sibling carries one; both network-bound) +- [x] opt-in network e2e (`DRIFT_TEST_NETWORK`): fetch example AOI untiled and with a small `tile_size`; assert per-year single-layer SpatRasters + `stac_items` attr + `.nc`/`.tif` extension routing, and **tiled == untiled** (tiled resampled onto the untiled grid, non-NA overlap) +- [x] `devtools::document()`; `lint` clean; `devtools::test()` 352 pass / 5 skip / 0 fail; code-check clean ## Phase 5: docs + gotchas note + NEWS + version - [ ] `inst/notes/gdalcubes-pc-gotchas.md`: tiling entry (fetch download bounded by tiling the cube_view; tiled mosaic cached as `.tif` via terra; #36) diff --git a/tests/testthat/test-dft_stac_fetch.R b/tests/testthat/test-dft_stac_fetch.R index c835398..13ca84e 100644 --- a/tests/testthat/test-dft_stac_fetch.R +++ b/tests/testthat/test-dft_stac_fetch.R @@ -183,3 +183,87 @@ test_that("tile_grid errors on a degenerate (empty) AOI", { tile_size = 500, res = 10) ) }) + +# --- mosaic_tiles(): reassemble per-tile rasters into one cache raster -------- +# Offline oracle for the tiled fetch (#36): res-aligned tiles that partition a +# reference grid must merge back into that grid, byte-for-byte. +test_that("mosaic_tiles merges res-aligned tiles back into the reference raster", { + ref <- terra::rast(nrows = 20, ncols = 20, xmin = 0, xmax = 200, + ymin = 0, ymax = 200, crs = "EPSG:32609") + terra::values(ref) <- seq_len(terra::ncell(ref)) # distinct code per cell + quads <- list(c(0, 100, 0, 100), c(100, 200, 0, 100), + c(0, 100, 100, 200), c(100, 200, 100, 200)) + tile_files <- vapply(quads, function(e) { + f <- tempfile(fileext = ".tif") + terra::writeRaster(terra::crop(ref, terra::ext(e[1], e[2], e[3], e[4])), f) + f + }, character(1)) + out <- tempfile(fileext = ".tif") + + drift:::mosaic_tiles(tile_files, out) + merged <- terra::rast(out) + + expect_equal(terra::nlyr(merged), 1L) + expect_equal( + c(terra::xmin(merged), terra::xmax(merged), + terra::ymin(merged), terra::ymax(merged)), + c(terra::xmin(ref), terra::xmax(ref), terra::ymin(ref), terra::ymax(ref)) + ) + expect_equal(terra::values(merged), terra::values(ref)) # exact reassembly + unlink(c(tile_files, out)) +}) + +test_that("mosaic_tiles handles a single tile", { + ref <- terra::rast(nrows = 5, ncols = 5, xmin = 0, xmax = 50, + ymin = 0, ymax = 50, crs = "EPSG:32609") + terra::values(ref) <- seq_len(25) + f <- tempfile(fileext = ".tif") + terra::writeRaster(ref, f) + out <- tempfile(fileext = ".tif") + + drift:::mosaic_tiles(f, out) + + expect_equal(terra::values(terra::rast(out)), terra::values(ref)) + unlink(c(f, out)) +}) + +# 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_fetch tiled result matches untiled over the AOI", { + skip_if(Sys.getenv("DRIFT_TEST_NETWORK") != "true", + "network test — set DRIFT_TEST_NETWORK=true to run") + skip_if_not_installed("gdalcubes") + aoi <- sf::st_read( + system.file("extdata", "example_aoi.gpkg", package = "drift"), + quiet = TRUE + ) + cache <- tempfile("drift_fetch_") + dir.create(cache) + + untiled <- dft_stac_fetch(aoi, source = "io-lulc", years = 2020, + cache_dir = cache)[["2020"]] + # small tile_size relative to the AOI bbox → several tiles, most bbox-only + # tiles dropped (the download-saving mechanism) + tiled_list <- dft_stac_fetch(aoi, source = "io-lulc", years = 2020, + tile_size = 500, cache_dir = cache) + tiled <- tiled_list[["2020"]] + + expect_false(is.null(attr(tiled_list, "stac_items"))) + expect_s4_class(tiled, "SpatRaster") + expect_equal(terra::nlyr(tiled), 1L) + # extension routing: untiled caches a gdalcubes .nc, tiled a terra .tif + expect_length(list.files(file.path(cache, "io-lulc"), + pattern = "^2020_.*\\.nc$"), 1) + expect_length(list.files(file.path(cache, "io-lulc"), + pattern = "^2020_.*\\.tif$"), 1) + # tiled == untiled over their common in-AOI cells: tiling changes only which + # bbox pixels are streamed, not the classification. Put the tiled mosaic onto + # the untiled grid (nearest — a no-op where the lattices coincide, robust to + # any sub-pixel offset gdalcubes gives the non-divisible untiled bbox) and + # compare where both are non-NA (the in-AOI overlap). + a <- terra::values(terra::resample(tiled, untiled, method = "near")) + b <- terra::values(untiled) + both <- !is.na(a) & !is.na(b) + expect_gt(sum(both), 0) + expect_equal(a[both], b[both]) +}) From 135e84ed8a38fb665b79cf734161066ae2b387a4 Mon Sep 17 00:00:00 2001 From: almac2022 Date: Thu, 9 Jul 2026 15:52:41 -0700 Subject: [PATCH 6/7] Release prep for #36: gotchas note, NEWS 0.6.0, version bump Phase 5 of #36. Document the download-tiling workaround in the gdalcubes gotchas note (the filter_geom-independent way to bound the read, cross-referencing the cube-path residual #38); add the NEWS 0.6.0 entry for tile_size; bump DESCRIPTION 0.5.0 -> 0.6.0. devtools::check() clean: 0 errors / 0 warnings / 0 notes. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01DKqnW6Fni3tc44Lo5D5fJz --- DESCRIPTION | 2 +- NEWS.md | 4 ++++ inst/notes/gdalcubes-pc-gotchas.md | 13 +++++++++++++ planning/active/progress.md | 6 +++++- planning/active/task_plan.md | 8 ++++---- 5 files changed, 27 insertions(+), 6 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 7c1561a..e6bfc84 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: drift Title: Detecting Riparian and Inland Floodplain Transitions -Version: 0.5.0 +Version: 0.6.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 3ebeae9..c73f71d 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,7 @@ +# drift 0.6.0 + +- `dft_stac_fetch()` gains `tile_size` (default `NULL`), an opt-in that bounds the STAC download to the AOI footprint (#36). By default a single cube is streamed over the whole AOI bounding box, so for a thin, diagonal floodplain corridor (measured ~10% of the bbox inside the polygon) roughly 10× more pixels are downloaded than the AOI needs. When `tile_size` (CRS units — metres for the default UTM CRS) is set, the bbox is split into a `res`-aligned grid and only tiles that intersect the AOI polygon are streamed, then mosaicked with `terra::merge()` — so a corridor fetches close to its footprint. Smaller tiles waste less bbox but cost more per-tile round trips (no auto-tuning). This is the `filter_geom`-independent path (the polygon-clip that would do this in the cube pipeline segfaults on the pinned gdalcubes build). Tiled fetches cache a terra GeoTIFF (`.tif`) rather than a gdalcubes NetCDF (`.nc`) and key distinctly, so existing untiled caches are untouched; `tile_size = NULL` is byte-for-byte the previous behavior. The same read residual on the continuous `dft_stac_cube()` path is tracked as #38. + # 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). diff --git a/inst/notes/gdalcubes-pc-gotchas.md b/inst/notes/gdalcubes-pc-gotchas.md index ef8e9e8..9b6351c 100644 --- a/inst/notes/gdalcubes-pc-gotchas.md +++ b/inst/notes/gdalcubes-pc-gotchas.md @@ -17,6 +17,19 @@ bfast 1.7.2. `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. +- **Download-side workaround without `filter_geom`: tile the `cube_view` (#36).** + Since `filter_geom` can't push the AOI into the read, the categorical + `dft_stac_fetch(tile_size = )` splits the AOI bbox into a `res`-aligned + grid and streams only tiles that intersect the AOI polygon (skipping the empty + bbox corners), then mosaics with `terra::merge()`. For a thin, diagonal + floodplain corridor (measured ~10% of the bbox inside the polygon) this fetches + near the AOI footprint instead of the full bbox. Tiles must be snapped to a + multiple of `res` and anchored at the bbox lower-left so their pixel grids are + co-lattice — otherwise the merge seams. The tiled mosaic is written with + `terra::writeRaster()` to a **`.tif`** (terra's NetCDF *write* is fragile — see + the round-trip bullet below), so tiled and untiled fetches cache under different + extensions and keys. The same read residual on the continuous `dft_stac_cube()` + path (its `cube_view` still streams the full bbox) is tracked as **#38**. - **`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/planning/active/progress.md b/planning/active/progress.md index 86017ff..4e00a84 100644 --- a/planning/active/progress.md +++ b/planning/active/progress.md @@ -33,4 +33,8 @@ untiled, resampled onto the untiled grid). Roxygen `@param tile_size` + cache doc. 352 pass / 5 skip; document + lint + code-check clean. Code-check flagged the network test's sub-pixel-offset fragility → hardened via near-resample. -- Next: Phase 5 (gotchas note + NEWS 0.6.0 + DESCRIPTION bump). +- Phase 5 done: gotchas-note tiling entry (download-side workaround for the + blocked filter_geom; .tif mosaic; cross-ref #38), NEWS 0.6.0 entry, DESCRIPTION + 0.5.0 → 0.6.0. `devtools::check()` clean: 0 errors / 0 warnings / 0 notes, + vignettes rebuild OK. +- Next: Phase 6 (/planning-archive → /gh-pr-push → /gh-pr-merge release v0.6.0). diff --git a/planning/active/task_plan.md b/planning/active/task_plan.md index 7a59153..bd082e2 100644 --- a/planning/active/task_plan.md +++ b/planning/active/task_plan.md @@ -35,12 +35,12 @@ tracked separately as #38 (out of scope). - [x] `devtools::document()`; `lint` clean; `devtools::test()` 352 pass / 5 skip / 0 fail; code-check clean ## Phase 5: docs + gotchas note + NEWS + version -- [ ] `inst/notes/gdalcubes-pc-gotchas.md`: tiling entry (fetch download bounded by tiling the cube_view; tiled mosaic cached as `.tif` via terra; #36) -- [ ] `NEWS.md` `# drift 0.6.0` — `tile_size` (opt-in, default `NULL` = unchanged; bounds download for sparse AOIs; tiled fetches cache as `.tif`); Closes #36 -- [ ] `DESCRIPTION` `0.5.0 → 0.6.0` + `Date` (final commit) +- [x] `inst/notes/gdalcubes-pc-gotchas.md`: tiling entry (fetch download bounded by tiling the cube_view; tiled mosaic cached as `.tif` via terra; #36) +- [x] `NEWS.md` `# drift 0.6.0` — `tile_size` (opt-in, default `NULL` = unchanged; bounds download for sparse AOIs; tiled fetches cache as `.tif`); Closes #36 +- [x] `DESCRIPTION` `0.5.0 → 0.6.0` + `Date` (final commit) ## Phase 6: validate, archive, PR, release -- [ ] `devtools::test()` / `lint` / `document` / `check` clean (network tests skip) +- [x] `devtools::test()` / `lint` / `document` / `check` clean (network tests skip) — check: 0 errors / 0 warnings / 0 notes - [ ] `/planning-archive`; `/gh-pr-push` (`Fixes #36`, `Relates to NewGraphEnvironment/sred-2025-2026#16`) - [ ] `/gh-pr-merge` → release v0.6.0 From acbbe99a56caf7698b5cb26f93525b5f917ac211 Mon Sep 17 00:00:00 2001 From: almac2022 Date: Thu, 9 Jul 2026 15:54:04 -0700 Subject: [PATCH 7/7] Archive planning files for issue #36 Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01DKqnW6Fni3tc44Lo5D5fJz --- .../2026-07-issue-36-tile-stac-fetch/.gitkeep | 0 .../README.md | 42 +++++++++++++++++++ .../findings.md | 0 .../progress.md | 0 .../task_plan.md | 0 5 files changed, 42 insertions(+) create mode 100644 planning/archive/2026-07-issue-36-tile-stac-fetch/.gitkeep create mode 100644 planning/archive/2026-07-issue-36-tile-stac-fetch/README.md rename planning/{active => archive/2026-07-issue-36-tile-stac-fetch}/findings.md (100%) rename planning/{active => archive/2026-07-issue-36-tile-stac-fetch}/progress.md (100%) rename planning/{active => archive/2026-07-issue-36-tile-stac-fetch}/task_plan.md (100%) diff --git a/planning/archive/2026-07-issue-36-tile-stac-fetch/.gitkeep b/planning/archive/2026-07-issue-36-tile-stac-fetch/.gitkeep new file mode 100644 index 0000000..e69de29 diff --git a/planning/archive/2026-07-issue-36-tile-stac-fetch/README.md b/planning/archive/2026-07-issue-36-tile-stac-fetch/README.md new file mode 100644 index 0000000..ec0fd14 --- /dev/null +++ b/planning/archive/2026-07-issue-36-tile-stac-fetch/README.md @@ -0,0 +1,42 @@ +## Outcome + +Added an opt-in `tile_size` to `dft_stac_fetch()` (#36) that bounds the STAC +download to the AOI footprint. By default `dft_stac_fetch()` streams one +gdalcubes cube over the whole AOI bounding box, so a thin, diagonal floodplain +corridor downloads ~10× more pixels than the polygon needs (measured 10.1% of +bbox cells inside a real io-lulc floodplain). When `tile_size` is set, the bbox +is split into a `res`-aligned grid, only tiles that intersect the AOI polygon +are streamed (via the shared `fetch_extent_to()` primitive), and the results are +mosaicked with `terra::merge(terra::sprc(...))` into a `.tif` cache. This is the +`filter_geom`-independent path — the polygon clip that would do this in the cube +pipeline segfaults on the pinned gdalcubes build (#32). + +Delivered tests-first across five atomic phases: `tile_grid()` + `tile_size_check()` +(offline), a conditional cache-key append with a frozen golden-hash regression +(`tile_size = NULL` reproduces the exact legacy key so existing untiled caches +stay valid), extraction of the shared fetch primitive, the tiled branch with a +`.tif`/`.nc` extension split and a GDAL `/vsicurl` config scoped to the tiled +path, and the release docs. + +**What was learned / decided:** +- terra's NetCDF *write* is fragile on the pinned stack, so the tiled mosaic is + written as a GeoTIFF (mirroring `dft_stac_cube()`), not a `.nc`. Cache + extension routes on `is.null(tile_size)`; the cache key already keys the two + apart, so the two formats never collide. +- Tiles must be snapped to a multiple of `res` and anchored at the bbox + lower-left so per-tile pixel grids are co-lattice — otherwise `terra::merge()` + seams. Boundary tiles are left un-trimmed (the overhang is masked away). +- `tile_size` is normalized once up front so the path gate (`is.null`) and the + cache-key append derive from the same snapped scalar (the #32 normalize-once + convention — no gate/key desync). +- Two Plan/code-check agent reviews caught the `.tif`-not-`.nc` blocker and a + sub-pixel-lattice fragility in the opt-in network test (hardened via + near-resample onto the untiled grid). +- Fetch-time streaming on the continuous `dft_stac_cube()` path has the same + residual, tracked separately as **#38** (unchanged by this work). + +`devtools::check()` clean (0 errors / 0 warnings / 0 notes); 352 pass / 5 skip. +Released as **v0.6.0**. + +Closed by: commits 127de94..135e84e on branch +`36-tile-dft-stac-fetch-to-bound-download-over-spars` / PR (Fixes #36) → v0.6.0. diff --git a/planning/active/findings.md b/planning/archive/2026-07-issue-36-tile-stac-fetch/findings.md similarity index 100% rename from planning/active/findings.md rename to planning/archive/2026-07-issue-36-tile-stac-fetch/findings.md diff --git a/planning/active/progress.md b/planning/archive/2026-07-issue-36-tile-stac-fetch/progress.md similarity index 100% rename from planning/active/progress.md rename to planning/archive/2026-07-issue-36-tile-stac-fetch/progress.md diff --git a/planning/active/task_plan.md b/planning/archive/2026-07-issue-36-tile-stac-fetch/task_plan.md similarity index 100% rename from planning/active/task_plan.md rename to planning/archive/2026-07-issue-36-tile-stac-fetch/task_plan.md