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/R/dft_stac_fetch.R b/R/dft_stac_fetch.R index 392f116..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,39 +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...") - 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 - ) - cube <- gdalcubes::raster_cube(col, v) - gdalcubes::write_ncdf(cube, cache_file, overwrite = TRUE) - 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) @@ -164,18 +202,134 @@ 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) +} + + +#' 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"]]) + }) +} + + +#' 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 +} + + +#' 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 } 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/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/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/archive/2026-07-issue-36-tile-stac-fetch/findings.md b/planning/archive/2026-07-issue-36-tile-stac-fetch/findings.md new file mode 100644 index 0000000..1496c5f --- /dev/null +++ b/planning/archive/2026-07-issue-36-tile-stac-fetch/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/archive/2026-07-issue-36-tile-stac-fetch/progress.md b/planning/archive/2026-07-issue-36-tile-stac-fetch/progress.md new file mode 100644 index 0000000..4e00a84 --- /dev/null +++ b/planning/archive/2026-07-issue-36-tile-stac-fetch/progress.md @@ -0,0 +1,40 @@ +# 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. +- 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. +- 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). +- 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). +- 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. +- 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/archive/2026-07-issue-36-tile-stac-fetch/task_plan.md b/planning/archive/2026-07-issue-36-tile-stac-fetch/task_plan.md new file mode 100644 index 0000000..bd082e2 --- /dev/null +++ b/planning/archive/2026-07-issue-36-tile-stac-fetch/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) +- [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) +- [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) +- [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 +- [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 +- [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 +- [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 + +## Validation +- [ ] Tests pass (`devtools::test()`); network 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_fetch.R b/tests/testthat/test-dft_stac_fetch.R index f0fbe2b..13ca84e 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,8 +80,190 @@ 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) 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) + ) +}) + +# --- 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]) +})