diff --git a/DESCRIPTION b/DESCRIPTION
index e6bfc84..eb2063a 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: drift
 Title: Detecting Riparian and Inland Floodplain Transitions
-Version: 0.6.0
-Date: 2026-07-09
+Version: 0.7.0
+Date: 2026-07-11
 Authors@R: c(
     person("Allan", "Irvine", , "al@newgraphenvironment.com", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0002-3495-2128")),
diff --git a/NEWS.md b/NEWS.md
index c73f71d..a39c608 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,7 @@
+# drift 0.7.0
+
+- `dft_stac_cube()` gains `tile_size` (default `NULL`), the continuous-path twin of `dft_stac_fetch()`'s `tile_size` (#36): an opt-in that bounds the STAC *read* to the AOI footprint (#38). By default one gdalcubes cube is streamed over the whole AOI bounding box, so the COG streaming — the dominant cost (~10-30 min for a multi-year monthly Sentinel-2 fetch) — scales with the bbox, not the AOI; for a thin, diagonal floodplain corridor the bbox is largely empty. 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 — each carrying the full SCL mask, spectral index, and 2022 baseline-offset split — then mosaicked with `terra::merge()`, so a corridor reads close to its footprint. This is the `filter_geom`-independent path (the polygon clip that would do this in-cube segfaults on the pinned gdalcubes build). The cube always caches a `.tif` and a tiled read keys distinctly, so untiled caches are untouched and `tile_size = NULL` is byte-for-byte the previous behavior. Because the cube resamples with bilinear, a tiled cube faithfully reproduces the untiled cube (bilinear-aligned correlation ~0.997, per-layer means within ~1e-3, no tile seams) but lands on a bbox-anchored grid that is sub-pixel-offset from — not pixel-identical to — the untiled cube; the offset is immaterial to the per-pixel `dft_rast_break()` / `dft_rast_trend()` reducers. With `tile_size` set, `clip = FALSE` returns the AOI-intersecting tile union (with `NA` where empty tiles were skipped), not a gap-free bounding box.
+
 # 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.
diff --git a/R/dft_stac_cube.R b/R/dft_stac_cube.R
index f00d921..06b59f8 100644
--- a/R/dft_stac_cube.R
+++ b/R/dft_stac_cube.R
@@ -9,7 +9,8 @@
 #'
 #' The index stack is materialized once to a GeoTIFF under [dft_cache_path()]
 #' as `<source>/cube_<key>.tif`, keyed by a hash of the AOI geometry and every
-#' cube-affecting parameter. Because it is invariant to [dft_rast_break()]'s
+#' cube-affecting parameter (including `clip` and `tile_size`, so a tiled read
+#' keys apart from an untiled one). Because it is invariant to [dft_rast_break()]'s
 #' parameters, caching it here makes bfast parameter sweeps cheap — they re-read
 #' the local raster instead of re-streaming COGs.
 #'
@@ -40,10 +41,13 @@
 #' @param clip Logical. When `TRUE` (default), clip the returned stack to the AOI
 #'   polygon with `terra::mask()` (cells outside → `NA` on every layer), so
 #'   [dft_rast_break()] / [dft_rast_trend()] reduce only in-polygon pixels. Set
-#'   `FALSE` to keep the full bounding box (e.g. for surrounding context, or to
-#'   mask later with a different polygon). Note this clips the *output* only — the
-#'   full bbox of COGs is still streamed either way (the AOI cannot be pushed into
-#'   the read on the pinned gdalcubes build; see `inst/notes/gdalcubes-pc-gotchas.md`).
+#'   `FALSE` to keep the wider extent (e.g. for surrounding context, or to mask
+#'   later with a different polygon). This clips the *output* — with the default
+#'   `tile_size = NULL` the full bbox of COGs is still streamed either way, so
+#'   `clip = FALSE` returns the full bounding box. When `tile_size` is set the read
+#'   is tiled, so `clip = FALSE` returns the **AOI-intersecting tile union** (a
+#'   stair-stepped superset of the polygon with `NA` where empty tiles were
+#'   skipped), not a gap-free bounding box.
 #' @param cloud_cover_max Numeric. Scene-level `eo:cloud_cover` maximum percent
 #'   for the STAC pre-filter (default 60).
 #' @param months Integer vector of calendar months (1-12) to keep, or `NULL`
@@ -57,6 +61,21 @@
 #' @param mask_values Integer vector of mask-band classes to exclude. When
 #'   `NULL`, uses `mask_values` from [dft_stac_config()] (e.g. Sentinel-2 SCL
 #'   cloud / shadow / cirrus classes).
+#' @param tile_size Numeric or `NULL` (default). Edge length, in CRS units
+#'   (metres for the default UTM CRS), of the read-tiling grid (#38). When `NULL`,
+#'   one cube is streamed over the whole AOI bounding box (the read 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) reads close to its
+#'   footprint. Snapped to a multiple of `res`. Smaller tiles waste less bbox but
+#'   cost more per-tile round trips; there is no auto-tuning. The cube always caches
+#'   a `.tif` either way; a tiled read keys distinctly (see the caching note above),
+#'   so untiled caches are untouched and `tile_size = NULL` is byte-for-byte the
+#'   previous behavior. This is the continuous-path twin of [dft_stac_fetch()]'s
+#'   `tile_size` — the `filter_geom`-independent way to bound the read. Because the
+#'   cube resamples with bilinear, a tiled cube faithfully reproduces the untiled
+#'   cube (the per-pixel reducers are unaffected) but lands on a bbox-anchored grid
+#'   that is sub-pixel-offset from — not pixel-identical to — the untiled cube.
 #' @param cache_dir Character. Cache directory. When `NULL`, uses
 #'   [dft_cache_path()].
 #' @param force Logical. Re-fetch even if cached, overwriting the cached raster
@@ -68,7 +87,8 @@
 #'   time value per layer — cached as a GeoTIFF. By default (`clip = TRUE`) the
 #'   stack is clipped to the AOI polygon (cloud-masked, cells outside the polygon
 #'   `NA`), so the reduced raster from [dft_rast_break()] is already polygon-tight;
-#'   pass `clip = FALSE` for the full AOI **bounding box**. For sources with a
+#'   pass `clip = FALSE` for the full AOI **bounding box** (or, with `tile_size`
+#'   set, the AOI-intersecting tile union). For sources with a
 #'   reflectance-offset baseline boundary (Sentinel-2), items are split at the
 #'   boundary and offset-corrected per side, so a series crossing it carries no
 #'   artificial index step.
@@ -105,6 +125,7 @@ dft_stac_cube <- function(aoi,
                           cloud_cover_max = 60,
                           months = NULL,
                           mask_values = NULL,
+                          tile_size = NULL,
                           cache_dir = NULL,
                           force = FALSE,
                           sign_fn = rstac::sign_planetary_computer()) {
@@ -150,6 +171,10 @@ dft_stac_cube <- function(aoi,
   # truthy-but-non-TRUE clip (e.g. 1 or "TRUE") must not skip the mask yet key as
   # TRUE, which would let a later clip=TRUE read the unclipped cube (#32).
   clip <- isTRUE(as.logical(clip))
+  # normalize tile_size ONCE (snap to a multiple of res) so the path gate
+  # (is.null) and the cache-key append derive from the same scalar (#36/#38).
+  # tile_size_check() is the shared download-tiling helper from #36.
+  if (!is.null(tile_size)) tile_size <- tile_size_check(tile_size, res)
 
   # Ensure aoi is sf
   if (inherits(aoi, "SpatVector")) aoi <- sf::st_as_sf(aoi)
@@ -181,7 +206,8 @@ dft_stac_cube <- function(aoi,
   cache_key <- stac_cube_cache_key(
     aoi_target, res, target_crs, dt, aggregation, resampling,
     cfg$stac_url, cfg$collection, band_assets, datetime, index,
-    cloud_cover_max, mask_values, scale, offset, months, offset_before, clip
+    cloud_cover_max, mask_values, scale, offset, months, offset_before, clip,
+    tile_size
   )
   cache_file <- file.path(cache_source_dir, paste0("cube_", cache_key, ".tif"))
 
@@ -228,24 +254,31 @@ dft_stac_cube <- function(aoi,
   message("  ", n_items, " items returned")
   if (n_items == 0) stop("No STAC items found for ", cfg$collection)
 
-  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 = t0, t1 = t1
-    ),
-    dx = res, dy = res, dt = dt,
-    aggregation = aggregation, resampling = resampling
-  )
+  # Baseline-conditional offset: split items at the boundary so each side is
+  # corrected with its own offset (below). The split is by item date, so it is
+  # the same for every read extent — compute it (and announce it) once, before
+  # assembling any cube.
+  is_pre <- rep(FALSE, length(items$features))
+  if (!is.null(offset_boundary)) {
+    item_date <- as.Date(substr(
+      vapply(items$features, function(f) f$properties$datetime %||% NA_character_, ""),
+      1, 10
+    ))
+    is_pre <- !is.na(item_date) & item_date < as.Date(offset_boundary)
+  }
+  if (any(is_pre) && !all(is_pre)) {
+    message("  offset split at ", offset_boundary, ": ",
+            sum(is_pre), " pre / ", sum(!is_pre), " post")
+  }
 
-  # Build the index cube for one item subset with one offset, materialize it, and
-  # read it back as a terra stack. The cube spans the AOI bounding box:
-  # gdalcubes::filter_geom() to clip to the polygon yields an all-NA cube (and can
-  # crash the compute worker) on the pinned build, so we mask clouds here and clip
-  # the assembled stack to the AOI polygon afterward with terra::mask() (see
-  # stac_cube_clip, #32), as the sibling dft_stac_fetch() does — never filter_geom().
-  build_index_stack <- function(features, offset_use) {
+  # Build the masked index cube for one item subset with one offset over a given
+  # cube_view, materialize it, and read it back as a terra stack. The cube spans
+  # the view's extent: gdalcubes::filter_geom() to clip to the polygon yields an
+  # all-NA cube (and can crash the compute worker) on the pinned build, so we mask
+  # clouds here and clip the assembled stack to the AOI polygon afterward with
+  # terra::mask() (see stac_cube_clip, #32), as dft_stac_fetch() does — never
+  # filter_geom().
+  build_index_stack <- function(features, offset_use, v) {
     img_col <- gdalcubes::stac_image_collection(
       features, asset_names = c(band_assets, mask_asset)
     )
@@ -259,27 +292,55 @@ dft_stac_cube <- function(aoi,
     terra::rast(tmp)
   }
 
-  # Baseline-conditional offset: split items at the boundary and correct each
-  # side with its own offset, then coalesce onto the shared monthly grid. Both
-  # subcubes are built over the full view so their layers align for terra::cover.
-  is_pre <- rep(FALSE, length(items$features))
-  if (!is.null(offset_boundary)) {
-    item_date <- as.Date(substr(
-      vapply(items$features, function(f) f$properties$datetime %||% NA_character_, ""),
-      1, 10
-    ))
-    is_pre <- !is.na(item_date) & item_date < as.Date(offset_boundary)
+  # Assemble the full masked index stack for one space extent: build the cube_view
+  # over that extent (same t0/t1/dt for every extent, so every tile yields the same
+  # nlyr), run the offset split, and coalesce the pre/post subcubes with
+  # terra::cover (both built over the same view so their layers align). A local
+  # closure (not @noRd) because it reads the call's items/offset/index/etc.
+  assemble_index_stack <- function(extent) {
+    v <- gdalcubes::cube_view(
+      srs = target_crs,
+      extent = list(
+        left = extent$left, right = extent$right,
+        bottom = extent$bottom, top = extent$top,
+        t0 = t0, t1 = t1
+      ),
+      dx = res, dy = res, dt = dt,
+      aggregation = aggregation, resampling = resampling
+    )
+    if (any(is_pre) && !all(is_pre)) {
+      terra::cover(
+        build_index_stack(items$features[is_pre], offset_before, v),
+        build_index_stack(items$features[!is_pre], offset, v)
+      )
+    } else {
+      build_index_stack(items$features,
+                        if (all(is_pre)) offset_before else offset, v)
+    }
   }
 
-  if (any(is_pre) && !all(is_pre)) {
-    message("  offset split at ", offset_boundary, ": ",
-            sum(is_pre), " pre / ", sum(!is_pre), " post")
-    stk <- terra::cover(
-      build_index_stack(items$features[is_pre], offset_before),
-      build_index_stack(items$features[!is_pre], offset)
-    )
+  # Untiled (tile_size = NULL): one cube over the AOI bounding box — unchanged
+  # behavior. Tiled (#38): stream only the res-aligned tiles that intersect the
+  # AOI polygon and mosaic them, so a sparse corridor reads near its footprint
+  # instead of the full bbox. tile_grid()/tile_size_check() are the shared
+  # download-tiling helpers from #36 (defined in dft_stac_fetch.R); the GDAL
+  # /vsicurl tuning set at the top of this function already covers the extra
+  # per-tile COG opens.
+  bbox_ext <- list(
+    left = bbox_target[["xmin"]], right = bbox_target[["xmax"]],
+    bottom = bbox_target[["ymin"]], top = bbox_target[["ymax"]]
+  )
+  if (is.null(tile_size)) {
+    stk <- assemble_index_stack(bbox_ext)
   } else {
-    stk <- build_index_stack(items$features, if (all(is_pre)) offset_before else offset)
+    tiles <- tile_grid(aoi_target, tile_size, res)
+    message("  tiling read into ", length(tiles), " tile(s) intersecting the AOI")
+    tile_stacks <- lapply(tiles, assemble_index_stack)
+    # every tile shares t0/t1/dt so nlyr is uniform; guard so a future per-tile
+    # time bound fails legibly here rather than deep inside terra::merge.
+    # terra::nlyr() returns a double, so the vapply template is numeric(1).
+    stopifnot(length(unique(vapply(tile_stacks, terra::nlyr, numeric(1)))) == 1L)
+    stk <- mosaic_stacks(tile_stacks)
   }
 
   # Restore the AOI-polygon clip removed in #30: mask the assembled stack (never
@@ -311,6 +372,24 @@ stac_cube_clip <- function(stk, aoi) {
 }
 
 
+#' Mosaic per-tile index stacks into one raster (in-memory, multi-layer)
+#'
+#' The reassembly step of the tiled read path (#38): each AOI-intersecting tile
+#' is streamed and reduced to its own multi-layer monthly index stack, and this
+#' merges them into one. Tiles are `res`-aligned and non-overlapping (see
+#' `tile_grid()`), so `terra::merge()` reassembles without resampling; it merges
+#' layer-by-layer positionally, so all `nlyr` layers are preserved in order
+#' (layer names/time are set by the caller after the merge). Unlike the fetch
+#' sibling's file-based `mosaic_tiles()`, this takes in-memory `SpatRaster`
+#' stacks (a tile may be the `terra::cover()` of a pre/post offset split) and
+#' returns the merged stack for the caller to clip and write. Skipped-tile gaps
+#' become `NA`, so the mosaic is the tile union, not a gap-free bounding box.
+#' @noRd
+mosaic_stacks <- function(stacks) {
+  terra::merge(terra::sprc(stacks))
+}
+
+
 #' Cache key for one STAC index-cube parameter set
 #'
 #' Cube-mode analogue of `stac_cache_key()` (kept separate so the fetch key
@@ -320,21 +399,28 @@ stac_cube_clip <- function(stk, aoi) {
 #' `scale`/`offset` are included because they change pixel values. `clip` is
 #' included because it changes the written extent (polygon vs bbox), so a
 #' `clip = FALSE` request must not read a clipped cube (or vice versa).
+#' `tile_size` (the download-tiling grid, #38) is appended to the hash ONLY when
+#' non-NULL, so an untiled cube keeps the exact legacy 18-element hash (existing
+#' `cube_<key>.tif` stay valid) while a tiled read keys distinctly. It must
+#' arrive already snapped by the caller (see `tile_size_check()`).
 #' @noRd
 stac_cube_cache_key <- function(aoi_target, res, target_crs, dt, aggregation,
                                 resampling, stac_url, collection, band_assets,
                                 datetime, index, cloud_cover_max, mask_values,
                                 scale, offset, months = NULL,
-                                offset_before = 0, clip = TRUE) {
+                                offset_before = 0, clip = TRUE,
+                                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, band_assets, datetime, index,
-      as.numeric(cloud_cover_max), sort(as.numeric(mask_values)),
-      as.numeric(scale), as.numeric(offset), sort(as.numeric(months)),
-      as.numeric(offset_before), as.logical(clip)
-    )),
-    1, 12
+  parts <- list(
+    geom_wkb, as.numeric(res), target_crs, dt, aggregation, resampling,
+    stac_url, collection, band_assets, datetime, index,
+    as.numeric(cloud_cover_max), sort(as.numeric(mask_values)),
+    as.numeric(scale), as.numeric(offset), sort(as.numeric(months)),
+    as.numeric(offset_before), as.logical(clip)
   )
+  # A tiled read caches the same .tif shape but over the AOI-intersecting tile
+  # union; keying it apart stops a tiled cube being served for an untiled request
+  # (or vice versa). Appending only when non-NULL preserves the legacy key.
+  if (!is.null(tile_size)) parts <- c(parts, list(as.numeric(tile_size)))
+  substr(rlang::hash(parts), 1, 12)
 }
diff --git a/inst/notes/gdalcubes-pc-gotchas.md b/inst/notes/gdalcubes-pc-gotchas.md
index 9b6351c..2d19cc1 100644
--- a/inst/notes/gdalcubes-pc-gotchas.md
+++ b/inst/notes/gdalcubes-pc-gotchas.md
@@ -13,10 +13,11 @@ bfast 1.7.2.
   the AOI polygon client-side (helper `stac_cube_clip()` = `terra::mask(stk,
   terra::vect(aoi))`), so the cube is polygon-tight and
   `dft_rast_break()`/`dft_rast_trend()` skip out-of-AOI pixels via their
-  `rowSums(!is.na) >= min_obs` gate. **Residual:** this clips the *output* only —
-  `cube_view(extent = bbox)` still streams the full bbox of COGs, so fetch time is
-  unchanged; pushing the AOI into the read would need a working `filter_geom` or
-  server-side windowing. `clip = FALSE` keeps the full bbox.
+  `rowSums(!is.na) >= min_obs` gate. **Residual (resolved, #38):** this clips the
+  *output* only — `cube_view(extent = bbox)` streams the full bbox of COGs, so the
+  clip alone does not cut fetch time. `dft_stac_cube(tile_size = <metres>)` now bounds
+  the *read* by tiling the `cube_view` (see the next bullet). `clip = FALSE` keeps the
+  full bbox (or, with `tile_size`, the AOI-intersecting tile union).
 - **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 = <metres>)` splits the AOI bbox into a `res`-aligned
@@ -28,8 +29,26 @@ bfast 1.7.2.
   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**.
+  extensions and keys. The continuous `dft_stac_cube(tile_size = <metres>)` (#38)
+  applies the same technique to the reflectance-cube read (per-tile `cube_view` +
+  SCL mask + the 2022 offset split + `terra::cover`, mosaicked by `mosaic_stacks()`;
+  the cube already caches `.tif`, so no extension split there).
+- **Bilinear tiling is not co-lattice with the untiled cube (#38).** The cube's
+  default `resampling = "bilinear"` makes the tiled read sensitive to grid alignment
+  in a way the categorical `near` path (#36) is not. gdalcubes *enlarges the untiled
+  bbox extent symmetrically* to align with `dx/dy` (observed ~0.5 px on a ~3.3 km
+  reach), while the tiles anchor at the bbox lower-left — so the tiled cube is **not
+  co-lattice** with the untiled cube and is **not pixel-identical** to it. It is still
+  a faithful resampling of the same source: bilinear-aligned correlation ~0.997,
+  per-layer means within ~1e-3, and **no tile seams** (gdalcubes reads the source
+  margin at tile edges, so |diff| at seams == interior). The only difference is a
+  benign sub-pixel grid offset, immaterial to the per-pixel `dft_rast_break()` /
+  `dft_rast_trend()` reducers. Consequence for tests/QA: compare a tiled cube to an
+  untiled one by **bilinear-aligning** one onto the other (`terra::resample(...,
+  method = "bilinear")`) and checking correlation + per-layer means — never
+  pixel-for-pixel. Re-anchoring the tiles to gdalcubes' enlarged origin to force
+  co-lattice was rejected: it would couple `tile_grid()` to gdalcubes internals and
+  change the shared #36 helper for no reducer-visible gain.
 - **`reduce_time()` R-callback runs in spawned worker processes at EVERY parallel
   setting** (incl. `parallel = 1`). A closure over enclosing locals fails there
   (`object 'band' not found`). Options: build a self-contained callback (inline
diff --git a/man/dft_stac_cube.Rd b/man/dft_stac_cube.Rd
index c64e3aa..6ba9911 100644
--- a/man/dft_stac_cube.Rd
+++ b/man/dft_stac_cube.Rd
@@ -18,6 +18,7 @@ dft_stac_cube(
   cloud_cover_max = 60,
   months = NULL,
   mask_values = NULL,
+  tile_size = NULL,
   cache_dir = NULL,
   force = FALSE,
   sign_fn = rstac::sign_planetary_computer()
@@ -52,10 +53,13 @@ must agree with.}
 \item{clip}{Logical. When \code{TRUE} (default), clip the returned stack to the AOI
 polygon with \code{terra::mask()} (cells outside → \code{NA} on every layer), so
 \code{\link[=dft_rast_break]{dft_rast_break()}} / \code{\link[=dft_rast_trend]{dft_rast_trend()}} reduce only in-polygon pixels. Set
-\code{FALSE} to keep the full bounding box (e.g. for surrounding context, or to
-mask later with a different polygon). Note this clips the \emph{output} only — the
-full bbox of COGs is still streamed either way (the AOI cannot be pushed into
-the read on the pinned gdalcubes build; see \code{inst/notes/gdalcubes-pc-gotchas.md}).}
+\code{FALSE} to keep the wider extent (e.g. for surrounding context, or to mask
+later with a different polygon). This clips the \emph{output} — with the default
+\code{tile_size = NULL} the full bbox of COGs is still streamed either way, so
+\code{clip = FALSE} returns the full bounding box. When \code{tile_size} is set the read
+is tiled, so \code{clip = FALSE} returns the \strong{AOI-intersecting tile union} (a
+stair-stepped superset of the polygon with \code{NA} where empty tiles were
+skipped), not a gap-free bounding box.}
 
 \item{cloud_cover_max}{Numeric. Scene-level \code{eo:cloud_cover} maximum percent
 for the STAC pre-filter (default 60).}
@@ -73,6 +77,22 @@ to fit a stable BFAST baseline.}
 \code{NULL}, uses \code{mask_values} from \code{\link[=dft_stac_config]{dft_stac_config()}} (e.g. Sentinel-2 SCL
 cloud / shadow / cirrus classes).}
 
+\item{tile_size}{Numeric or \code{NULL} (default). Edge length, in CRS units
+(metres for the default UTM CRS), of the read-tiling grid (#38). When \code{NULL},
+one cube is streamed over the whole AOI bounding box (the read 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) reads close to its
+footprint. Snapped to a multiple of \code{res}. Smaller tiles waste less bbox but
+cost more per-tile round trips; there is no auto-tuning. The cube always caches
+a \code{.tif} either way; a tiled read keys distinctly (see the caching note above),
+so untiled caches are untouched and \code{tile_size = NULL} is byte-for-byte the
+previous behavior. This is the continuous-path twin of \code{\link[=dft_stac_fetch]{dft_stac_fetch()}}'s
+\code{tile_size} — the \code{filter_geom}-independent way to bound the read. Because the
+cube resamples with bilinear, a tiled cube faithfully reproduces the untiled
+cube (the per-pixel reducers are unaffected) but lands on a bbox-anchored grid
+that is sub-pixel-offset from — not pixel-identical to — the untiled cube.}
+
 \item{cache_dir}{Character. Cache directory. When \code{NULL}, uses
 \code{\link[=dft_cache_path]{dft_cache_path()}}.}
 
@@ -87,7 +107,8 @@ A \link[terra:SpatRaster-class]{terra::SpatRaster} index stack — one layer per
 time value per layer — cached as a GeoTIFF. By default (\code{clip = TRUE}) the
 stack is clipped to the AOI polygon (cloud-masked, cells outside the polygon
 \code{NA}), so the reduced raster from \code{\link[=dft_rast_break]{dft_rast_break()}} is already polygon-tight;
-pass \code{clip = FALSE} for the full AOI \strong{bounding box}. For sources with a
+pass \code{clip = FALSE} for the full AOI \strong{bounding box} (or, with \code{tile_size}
+set, the AOI-intersecting tile union). For sources with a
 reflectance-offset baseline boundary (Sentinel-2), items are split at the
 boundary and offset-corrected per side, so a series crossing it carries no
 artificial index step.
@@ -103,7 +124,8 @@ breakpoint detection.
 \details{
 The index stack is materialized once to a GeoTIFF under \code{\link[=dft_cache_path]{dft_cache_path()}}
 as \verb{<source>/cube_<key>.tif}, keyed by a hash of the AOI geometry and every
-cube-affecting parameter. Because it is invariant to \code{\link[=dft_rast_break]{dft_rast_break()}}'s
+cube-affecting parameter (including \code{clip} and \code{tile_size}, so a tiled read
+keys apart from an untiled one). Because it is invariant to \code{\link[=dft_rast_break]{dft_rast_break()}}'s
 parameters, caching it here makes bfast parameter sweeps cheap — they re-read
 the local raster instead of re-streaming COGs.
 
diff --git a/planning/archive/2026-07-issue-38-tile-stac-cube/.gitkeep b/planning/archive/2026-07-issue-38-tile-stac-cube/.gitkeep
new file mode 100644
index 0000000..e69de29
diff --git a/planning/archive/2026-07-issue-38-tile-stac-cube/README.md b/planning/archive/2026-07-issue-38-tile-stac-cube/README.md
new file mode 100644
index 0000000..7bd12b2
--- /dev/null
+++ b/planning/archive/2026-07-issue-38-tile-stac-cube/README.md
@@ -0,0 +1,48 @@
+## Outcome
+
+Added an opt-in `tile_size` to `dft_stac_cube()` (#38) that bounds the STAC **read**
+to the AOI footprint — the continuous-path twin of #36's `dft_stac_fetch(tile_size=)`.
+By default one gdalcubes cube is streamed over the whole AOI bounding box, so for a
+thin, diagonal floodplain corridor (area/bbox ≈ 0.1) the COG streaming (~10–30 min for a
+multi-year monthly Sentinel-2 fetch) is ~10× larger than the polygon needs. When
+`tile_size` is set, the bbox is split into a `res`-aligned grid and only tiles
+intersecting the AOI polygon are streamed — each carrying the full SCL mask, spectral
+index, and 2022 baseline-offset split — then mosaicked with a new `mosaic_stacks()`. The
+`filter_geom`-independent path (the in-cube polygon clip segfaults on the pinned build).
+
+Delivered tests-first across six atomic phases: the cube cache key gained a conditional
+`tile_size` append guarded by a **frozen golden-hash** (`638a2be11fdf`, byte-for-byte
+preservation of the untiled key so existing `cube_*.tif` stay valid); `mosaic_stacks()`
+with three offline oracles (multi-layer merge, cover-then-merge == merge-then-cover
+commutativity, tile-union extent); the core refactor into an `assemble_index_stack(extent)`
+closure with the untiled path behavior-preserving and the tiled branch guarded on uniform
+nlyr; an opt-in network e2e; docs; and the release bump.
+
+**What was learned / decided (the network e2e earned its keep):**
+- `terra::nlyr()` returns a **double** — the uniform-nlyr guard's `vapply(..., integer(1))`
+  template errored on the live run; fixed to `numeric(1)`. The tiled fetch itself worked
+  (all 12 tiles streamed, offset-split, covered) — only the strict guard tripped.
+- The tiled and untiled cubes are **not co-lattice**: gdalcubes enlarges the untiled bbox
+  extent symmetrically to align with `dx/dy` (~0.5 px), while tiles anchor at the bbox
+  lower-left. So they cannot be compared pixel-for-pixel. Confirmed offline on saved real
+  cubes that this is a **benign grid offset, not a bug**: bilinear-aligned correlation
+  ≈ 0.997, per-layer means within ~1e-3, and **no tile seams** (edge |diff| == interior —
+  gdalcubes reads the source margin at tile edges). The network test was rewritten to
+  assert grid-robust equivalence (bilinear-aligned correlation + per-layer means), with
+  thresholds measured on the real cubes. Offset-split-under-tiling stays covered by the
+  offline commutativity oracle.
+- Tile anchoring stays at bbox-LL (re-anchoring to gdalcubes' enlarged origin would couple
+  `tile_grid()` to gdalcubes internals and change the shared #36 helper for no
+  reducer-visible gain). The bilinear-tiling lesson is recorded in
+  `inst/notes/gdalcubes-pc-gotchas.md`.
+- The cube already caches `.tif` and already sets the GDAL `/vsicurl` config, so #36's
+  `.nc`/`.tif` extension routing and scoped-config machinery were **not** needed here.
+- `tile_size_check()` / `tile_grid()` were reused verbatim from #36 (same package
+  namespace). Note: `object_usage_linter` resolves cross-file internals against the
+  *installed* namespace, so a stale install shows false positives — reinstall to clear.
+
+`devtools::check()` clean (0 errors / 0 warnings / 0 notes); 365 pass / 6 skip.
+Released as **v0.7.0**.
+
+Closed by: commits 90f9d93..d7a5094 on branch
+`38-bound-dft-stac-cube-streaming-to-the-aoi` / PR (Fixes #38) → v0.7.0.
diff --git a/planning/archive/2026-07-issue-38-tile-stac-cube/findings.md b/planning/archive/2026-07-issue-38-tile-stac-cube/findings.md
new file mode 100644
index 0000000..58ba571
--- /dev/null
+++ b/planning/archive/2026-07-issue-38-tile-stac-cube/findings.md
@@ -0,0 +1,108 @@
+# Findings — Bound `dft_stac_cube()` streaming to the AOI (#38)
+
+## Issue context
+
+`dft_stac_cube()` (`R/dft_stac_cube.R:231-240`) builds its gdalcubes cube over the AOI
+**bounding box**, so the COG streaming (~10-30 min for a multi-year monthly Sentinel-2
+fetch) scales with the bbox, not the AOI polygon. Measured area/bbox ≈ 0.105 on the
+packaged example AOI → ~10× streaming overhead for a corridor. #32 restored polygon-tight
+**output** (`clip = TRUE`, `terra::mask()`) but that runs after streaming — it does not
+reduce the read. #38 = reduce the READ by tiling the `cube_view`, mirroring #36 on the
+categorical sibling `dft_stac_fetch()`. `gdalcubes::filter_geom()` (the direct fix)
+segfaults on the pinned build (out of scope).
+
+## Design (approved plan + Plan-agent review)
+
+Reuses #36's `tile_size_check()` / `tile_grid()` verbatim (same package namespace).
+Cube path is harder per tile (SCL mask + index + 2022 offset-split + `terra::cover`;
+multi-layer stacks) but simpler in caching (always `.tif`; GDAL config already
+unconditional at `:116-130`). `mosaic_stacks()` = in-memory multi-layer
+`terra::merge(terra::sprc(...))`; distinct from #36's file-based single-layer
+`mosaic_tiles`. `assemble_index_stack(extent)` is a **local closure** (like the existing
+`build_index_stack`), NOT `@noRd` — it must close over `dft_stac_cube`'s call-locals.
+
+## Plan-agent review — issues incorporated
+
+- **B1 (Blocker):** `assemble_index_stack` must be a local closure, not a top-level
+  `@noRd` function — a top-level fn can't close over `items`/`is_pre`/`offset`/`t0`/… Only
+  `mosaic_stacks` is `@noRd`.
+- **O1 (Ordering):** the cube key has **no** frozen golden literal today (fetch froze
+  `79f67b7b9dae`; cube tests assert only determinism+distinctness). `stac_cube_cache_key`
+  is currently a flat `hash(list(...))`. Freeze the current untiled literal FIRST, then
+  migrate to the append pattern, then confirm unchanged — else every `cube_<key>.tif`
+  orphans and re-streams (10-30 min each).
+- **G1 (Gap):** cover-then-merge == merge-then-cover has no CI test; only the opt-in
+  network test would touch it. Add an offline synthetic commutativity test.
+- **G3 (Gap):** `clip = FALSE` docstring promises "full bounding box", but under tiling
+  the mosaic returns the AOI-intersecting **tile union** (stair-stepped, ⊊ full bbox).
+  Amend the docs; add an offline extent-semantics test.
+- **A1:** assert uniform nlyr across tiles before merge (`stopifnot`) for a legible
+  failure if a future edit derives per-tile time bounds.
+- **G2:** the per-tile temp NetCDFs — accept the pre-existing leak (untiled already leaks
+  1-2; tiling makes it 2×N of the same kind). Do NOT add unlink plumbing: terra's
+  NetCDF-backed rasters are lazy until `writeRaster`, so early unlink corrupts the mosaic.
+- **AC1:** kNDVI is **float** (bilinear + median), not integer classes — the network
+  equality test can't copy #36's integer `expect_equal`; needs an explicit tolerance /
+  robust quantile (~1e-6 FP jitter from differing gdalcubes chunk boundaries; possible
+  edge-cell deltas). And don't `compareGeom`-assert identical extents (tiled = tile-union,
+  untiled = full bbox).
+- **AC2:** the e2e window must **straddle 2022-01-25** so the offset-split-under-tiling
+  path actually runs (the existing cube network test uses a pre-boundary window).
+- **S1:** reuse the two shared helpers in place (zero churn); extraction to
+  `R/stac_tiling.R` is a future option if a third consumer appears.
+- **A2/A3:** multi-layer merge preserves positional layer order (names/time overwritten
+  post-merge); nlyr fixed by the shared `t0/t1/dt` axis (zero-overlap tiles → all-NA
+  layers, not fewer). In-polygon values identical after mask; only extent/dims differ.
+
+## Network e2e findings (Phase 4) — what only the live path surfaced
+
+The opt-in network test (32-min run over the 7-month straddling window) did its job
+and caught two things the offline oracles structurally could not (they never invoke
+`dft_stac_cube`, which is network-bound):
+
+1. **`terra::nlyr()` returns a double, not integer.** The uniform-nlyr guard
+   `vapply(tile_stacks, terra::nlyr, integer(1))` errored ("values must be type
+   'integer' but result is 'double'"). Fixed: `numeric(1)` template. The tiled fetch
+   itself worked — all 12 tiles streamed, offset-split, and covered without error; it
+   only tripped at the strict guard *after* assembly. (The untiled network test passed.)
+
+2. **Tiled and untiled cubes are NOT co-lattice.** gdalcubes enlarges the untiled
+   bbox extent *symmetrically* to align with `dx/dy` ("extent will be enlarged by
+   0.79 [x] / 4.88 [y] at both sides"), so the untiled grid origin is
+   `683391.9, 6027436` while the tile grid (anchored at bbox-LL, each tile a `res`
+   multiple so un-enlarged) is `683392.6, 6027441` — offset ~0.79m (x), ~4.88m (y,
+   nearly half a pixel). For integer LULC + `near` (#36) a sub-pixel shift is absorbed
+   (same source pixel), so #36's exact `expect_equal` held. For float kNDVI + bilinear
+   (#38) the shift changes every cell, and a `near`-resample comparison then shifts one
+   grid against the other → median |diff| 0.009, q99 0.11, max 0.24. This is a
+   **test-methodology** issue (both grids are valid; the tiled cube is not wrong), not a
+   code bug — the tiled cube lands on a bbox-anchored grid, sub-pixel-offset from the
+   untiled grid, immaterial to per-pixel trajectory reducers. Confirming via a
+   bilinear-aligned comparison (diagnose_grid.R): if bilinear-aligned diffs collapse to
+   ~1e-3 the values are correct and it was purely the offset.
+
+**Resolution (confirmed via saved 2021-07/08 cubes, offline):** bilinear alignment gives
+cor **0.997**, median |diff| **3.4e-3**, per-layer means agree to **6e-4**; and mean |diff|
+at tile seams (0.0075) **==** interior (0.0076) → **no seam artifacts** (gdalcubes reads the
+source margin at tile edges). So the tiled cube is a faithful resampling of the source,
+differing from untiled only by the benign sub-pixel grid offset. Fix is test-methodology,
+not code: the network test now aligns tiled onto the untiled grid with **bilinear** and
+asserts grid-robust equivalence — `cor > 0.98`, `median |diff| < 0.01`, per-layer mean
+agreement `< 0.01` — thresholds measured on the real cubes with headroom. Window switched
+to growing-season (`2021-07-01/2021-08-31`) for robust coverage; offset-split-under-tiling
+stays covered by the offline commutativity oracle. Code fix carried here: `terra::nlyr()`
+returns double → guard template `numeric(1)`.
+
+**Decision — tile anchoring stays at bbox-LL (not the gdalcubes-enlarged origin).**
+Re-anchoring tiles to match the untiled grid would (a) couple `tile_grid` to gdalcubes'
+internal extent-enlargement, (b) change the shared #36 helper, and (c) still differ at
+edges — for no scientific gain, since the tiled cube is already an equivalent resampling
+and downstream reducers are per-pixel. Documented as a benign grid offset in roxygen/NEWS.
+
+## Key line references
+- `R/dft_stac_cube.R`: cube_view `:231-240`; `build_index_stack` `:248-260`; offset-split +
+  cover `:265-283`; unified tail `:290-295`; `stac_cube_cache_key` `:324-340`.
+- `R/dft_stac_fetch.R`: `tile_size_check` `:233`; `tile_grid` `:271`; `stac_cache_key`
+  append pattern `:210-222`.
+- `tests/testthat/test-dft_stac_cube.R`: `cube_key()` helper `:31-45`; distinctness block
+  `:52-73`; network e2e `:117-140`.
diff --git a/planning/archive/2026-07-issue-38-tile-stac-cube/progress.md b/planning/archive/2026-07-issue-38-tile-stac-cube/progress.md
new file mode 100644
index 0000000..3387820
--- /dev/null
+++ b/planning/archive/2026-07-issue-38-tile-stac-cube/progress.md
@@ -0,0 +1,9 @@
+# Progress — Bound `dft_stac_cube()` streaming to the AOI (#38)
+
+## Session 2026-07-11
+
+- Plan-mode exploration — read `R/dft_stac_cube.R` + `R/dft_stac_fetch.R` (#36 tiling) in full; confirmed cube test patterns, S2 config (offset_boundary 2022-01-25), example AOI
+- Plan-agent review caught B1 (local-closure vs @noRd), O1 (author cube-key golden guardian before refactor), G1/G3 (offline commutativity + clip=FALSE extent tests), A1 (uniform-nlyr guard), AC1/AC2 (float tolerance + 2022-straddling network window) — all folded into the approved plan
+- Created branch `38-bound-dft-stac-cube-streaming-to-the-aoi` off main
+- Scaffolded PWF baseline with approved phases
+- Next: Phase 1 — freeze the cube-key golden guardian, then the conditional `tile_size` append
diff --git a/planning/archive/2026-07-issue-38-tile-stac-cube/task_plan.md b/planning/archive/2026-07-issue-38-tile-stac-cube/task_plan.md
new file mode 100644
index 0000000..bfe40a1
--- /dev/null
+++ b/planning/archive/2026-07-issue-38-tile-stac-cube/task_plan.md
@@ -0,0 +1,56 @@
+# Task: Bound `dft_stac_cube()` streaming to the AOI (#38)
+
+`dft_stac_cube()` builds its gdalcubes cube over the AOI **bounding box**
+(`cube_view(extent = bbox_target)`), so the COG streaming — the dominant cost,
+~10-30 min for a multi-year monthly Sentinel-2 fetch — scales with the bbox, not the
+AOI polygon. For a thin corridor floodplain the bbox is ~90% empty (area/bbox ≈ 0.105
+on the packaged example AOI → ~10× streaming overhead). #32 restored polygon-tight
+**output** via `terra::mask()` but that runs *after* streaming — it does not reduce the
+read. #38 reduces the **read** by tiling the `cube_view` (mirroring #36 for the
+categorical sibling `dft_stac_fetch()`): stream only tiles intersecting the AOI polygon,
+then mosaic. Opt-in `tile_size = NULL` default = today's behavior, byte-for-byte.
+
+Reuses #36's `tile_size_check()` / `tile_grid()` (`R/dft_stac_fetch.R`) in place. The
+cube path is harder (per-tile SCL mask + index + 2022 offset-split + `terra::cover`;
+multi-layer stacks) but simpler in caching (always `.tif`; GDAL config already
+unconditional — no `.nc`/`.tif` routing).
+
+## Phase 1: cube cache key — golden guardian first, then conditional `tile_size` (offline, tests-first)
+- [x] compute the current untiled `stac_cube_cache_key()` for `cube_key()`'s fixed inputs; freeze as `expect_equal(cube_key(), "638a2be11fdf")` — authored **before** touching the function
+- [x] extend the `cube_key()` helper with `tile_size = NULL` (pass through to `stac_cube_cache_key`)
+- [x] distinctness: `cube_key(tile_size = 500) != cube_key()`; distinct sizes → distinct keys; snap-before-key (`504` and `500` at res 10 → same key)
+- [x] refactor `stac_cube_cache_key` to the append-only shape (add trailing `tile_size = NULL`); append only when non-NULL; **confirm the frozen literal is unchanged** (byte-preserving). Existing 18-arg call site keeps working via the default — the `tile_size` param + normalize + call-site pass land in Phase 3 alongside the tiled read (so no distinct-key-but-bbox-read window)
+
+## Phase 2: `mosaic_stacks` + multi-layer merge / commutativity / extent oracles (offline, tests-first)
+- [x] multi-layer merge oracle: reference multi-layer SpatRaster → split into res-aligned non-overlapping tiles → `mosaic_stacks` → exact values per layer; nlyr + layer order preserved
+- [x] **commutativity:** synthetic multi-layer `pre`/`post` stacks with a known NA pattern → `mosaic_stacks(lapply(tiles, \(t) terra::cover(terra::crop(pre,t), terra::crop(post,t))))` equals `terra::cover(pre, post)` cell-for-cell (cover-then-merge == merge-then-cover, in CI)
+- [x] **extent semantics:** `mosaic_stacks` over `tile_grid`-derived synthetic tiles for the example AOI leaves NA gaps where empty tiles were skipped (n_kept 31 < n_full 49) — documents the `clip = FALSE` + `tile_size` narrowing
+- [x] implement `mosaic_stacks()` `@noRd`; tests green (41 pass, 0 fail; lint-clean)
+
+## Phase 3: tile the cube read — refactor + wire `tile_size` end-to-end
+- [x] add `tile_size = NULL` param to `dft_stac_cube`; normalize once via `tile_size_check` at the top; call site passes `tile_size = tile_size`
+- [x] `build_index_stack` gains `v`; add local closure `assemble_index_stack(extent)` (moves cube_view + offset-split + cover inside); untiled routes through it over `bbox_ext` — identical output
+- [x] tiled branch: `tile_grid` → per-tile `assemble_index_stack` → uniform-nlyr `stopifnot` → `mosaic_stacks`; unified clip/time/names/`.tif` tail; reuse #36 helpers in place (comment)
+- [x] roxygen `@param tile_size` + amended clip/read caveat (`clip=FALSE`+`tile_size` = tile-union extent) + cache-doc `tile_size` note
+- [x] `devtools::document()`; `lintr::lint_package()` R/ source clean (reinstalled stale v0.3.0 → object_usage false positives cleared); `devtools::test()` green (41 pass, offline)
+
+## Phase 4: opt-in network e2e (`DRIFT_TEST_NETWORK`) — ran live; caught 2 real issues
+- [x] ran the e2e (32-min straddling fetch) → caught: (1) `terra::nlyr()` returns **double** → guard template fixed to `numeric(1)`; (2) tiled/untiled are **not co-lattice** (gdalcubes enlarges the untiled bbox symmetrically ~0.5px; tiles anchor at bbox-LL) so pixel-identity is the wrong assertion
+- [x] confirmed correctness offline on saved real cubes: bilinear-aligned cor **0.997**, median |diff| **3.4e-3**, per-layer means agree **6e-4**, and **no tile seams** (edge |diff| == interior) → tiled cube is a faithful resampling, benign grid offset only
+- [x] rewrote the network test: grow-season window (`2021-07-01/2021-08-31`); assert `SpatRaster`, equal nlyr, time set, 2 `cube_<key>.tif`, kept-tile count < full grid; equivalence via **bilinear**-aligned `cor > 0.98` + `median |diff| < 0.01` + per-layer mean agreement `< 0.01` (thresholds measured with headroom). Offset-split-under-tiling covered by the offline commutativity oracle
+
+## Phase 5: docs + gotchas + NEWS + version
+- [x] `inst/notes/gdalcubes-pc-gotchas.md`: flipped the #38 residual to resolved-via-cube-tiling; updated the #36 bullet's cross-ref; added the bilinear-tiling not-co-lattice / no-seams lesson
+- [x] `NEWS.md` `# drift 0.7.0` — `dft_stac_cube()` gains `tile_size` (opt-in, default `NULL` = unchanged; bounds the read for sparse AOIs; faithful-but-grid-offset from untiled); Closes #38. Roxygen `@param tile_size` grid-offset note added
+- [x] `DESCRIPTION` `0.6.0 → 0.7.0` + `Date 2026-07-11` (final commit, after `devtools::check`)
+
+## Phase 6: validate, archive, PR, release
+- [x] `devtools::test()` (365 pass / 6 skip) / `lint` / `document` / `check` clean (0/0/0; network tests skip)
+- [ ] `/planning-archive`; `/gh-pr-push` (`Fixes #38`, `Relates to NewGraphEnvironment/sred-2025-2026#16`)
+- [ ] `/gh-pr-merge` → release v0.7.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_cube.R b/tests/testthat/test-dft_stac_cube.R
index efe3ea2..8d05170 100644
--- a/tests/testthat/test-dft_stac_cube.R
+++ b/tests/testthat/test-dft_stac_cube.R
@@ -36,11 +36,16 @@ cube_key <- function(aoi = square_aoi(), res = 10, target_crs = "EPSG:32609",
                      datetime = "2019-01-01/2023-12-31", index = "kndvi",
                      cloud_cover_max = 60, mask_values = c(3, 8, 9, 10, 11),
                      scale = 1e-4, offset = -0.1, months = NULL,
-                     offset_before = 0, clip = TRUE) {
+                     offset_before = 0, clip = TRUE, tile_size = NULL) {
+  # dft_stac_cube snaps tile_size to the pixel grid before it reaches the key;
+  # mirror that here so the snap-before-key test reflects real behavior (#38)
+  if (!is.null(tile_size)) {
+    tile_size <- suppressMessages(drift:::tile_size_check(tile_size, res))
+  }
   drift:::stac_cube_cache_key(
     aoi, res, target_crs, dt, aggregation, resampling, stac_url, collection,
     band_assets, datetime, index, cloud_cover_max, mask_values, scale, offset,
-    months, offset_before, clip
+    months, offset_before, clip, tile_size
   )
 }
 
@@ -49,6 +54,23 @@ test_that("stac_cube_cache_key is deterministic and 12-char hex", {
   expect_match(cube_key(), "^[0-9a-f]{12}$")
 })
 
+test_that("stac_cube_cache_key untiled key is frozen (legacy-cache guardian)", {
+  # Freezes the exact 12-char hash for cube_key()'s fixed inputs so the
+  # tile_size append can't silently perturb the untiled key and orphan every
+  # existing cube_<key>.tif (10-30 min to re-stream). Mirrors the fetch golden
+  # 79f67b7b9dae (#36). If this ever changes, existing cube caches are invalid.
+  expect_equal(cube_key(), "638a2be11fdf")
+})
+
+test_that("stac_cube_cache_key keys tile_size distinctly and after snapping", {
+  base <- cube_key()                                     # tile_size = NULL
+  expect_false(cube_key(tile_size = 500) == base)        # tiled keys apart from untiled
+  expect_false(cube_key(tile_size = 500) ==
+                 cube_key(tile_size = 250))              # distinct sizes -> distinct keys
+  # snapped to the res-lattice before the key: 504 -> 500 (res 10), same key
+  expect_equal(cube_key(tile_size = 504), cube_key(tile_size = 500))
+})
+
 test_that("stac_cube_cache_key changes with each cube-affecting parameter", {
   base <- cube_key()
   expect_false(cube_key(aoi = square_aoi(dx = 0.5)) == base)
@@ -112,6 +134,86 @@ test_that("stac_cube_clip masks cells outside the AOI polygon on every layer", {
   expect_true(all(is.na(vals[!inside, ])))    # outside polygon: NA
 })
 
+# mosaic_stacks(): the in-memory, multi-layer merge that reassembles per-tile
+# index stacks into one raster on the tiled read path (#38). Network-free —
+# synthetic multi-layer rasters split into res-aligned, non-overlapping tiles.
+
+test_that("mosaic_stacks reassembles res-aligned tiles losslessly across layers", {
+  # 3-layer reference on a known lattice; distinct values per cell AND per layer
+  # (values fill layer-major: L1 = 1:64, L2 = 65:128, L3 = 129:192), so a layer
+  # swap or a spatial mis-merge would change the compared values
+  ref <- terra::rast(nrows = 8, ncols = 8, xmin = 0, xmax = 80,
+                     ymin = 0, ymax = 80, crs = "EPSG:32609", nlyrs = 3)
+  terra::values(ref) <- seq_len(terra::ncell(ref) * terra::nlyr(ref))
+  # split into 4 quadrant tiles (40x40), res-aligned, non-overlapping, tiling ref
+  exts <- list(terra::ext(0, 40, 0, 40), terra::ext(40, 80, 0, 40),
+               terra::ext(0, 40, 40, 80), terra::ext(40, 80, 40, 80))
+  merged <- drift:::mosaic_stacks(lapply(exts, function(e) terra::crop(ref, e)))
+
+  expect_s4_class(merged, "SpatRaster")
+  expect_equal(terra::nlyr(merged), 3L)                       # all layers kept
+  expect_true(terra::compareGeom(merged, ref, stopOnError = FALSE))
+  expect_equal(terra::values(merged), terra::values(ref))     # exact, per layer
+})
+
+test_that("tiling commutes with the offset-split cover (per-tile cover + merge == global cover)", {
+  # Mimics the 2022 offset split: pre/post subcubes with complementary NA
+  # patterns, coalesced by terra::cover. Under tiling the cover runs per tile,
+  # then the tiles are merged; that must equal covering the full extent once.
+  base <- function() {
+    terra::rast(nrows = 8, ncols = 8, xmin = 0, xmax = 80,
+                ymin = 0, ymax = 80, crs = "EPSG:32609", nlyrs = 2)
+  }
+  pre <- base()
+  post <- base()
+  n <- terra::ncell(pre)                         # 64
+  odd <- seq_len(n) %% 2 == 1
+  # layer 1: pre holds odd cells (NA on even), post holds even; layer 2 reversed
+  # -> cover must pick pre-then-post per cell, per layer, and tiles partition space
+  terra::values(pre)  <- c(ifelse(odd, seq_len(n), NA),
+                           ifelse(odd, NA, seq_len(n) + 200))
+  terra::values(post) <- c(ifelse(odd, NA, seq_len(n) + 100),
+                           ifelse(odd, seq_len(n) + 300, NA))
+  ref <- terra::cover(pre, post)
+  expect_false(anyNA(terra::values(ref)))        # complementary -> fully coalesced
+
+  exts <- list(terra::ext(0, 40, 0, 40), terra::ext(40, 80, 0, 40),
+               terra::ext(0, 40, 40, 80), terra::ext(40, 80, 40, 80))
+  tiled <- drift:::mosaic_stacks(lapply(exts, function(e) {
+    terra::cover(terra::crop(pre, e), terra::crop(post, e))
+  }))
+
+  expect_true(terra::compareGeom(tiled, ref, stopOnError = FALSE))
+  expect_equal(terra::values(tiled), terra::values(ref))
+})
+
+test_that("mosaic_stacks over tile_grid tiles leaves NA gaps where empty tiles were skipped", {
+  # Documents the clip = FALSE + tile_size contract: the mosaic is the union of
+  # AOI-intersecting tiles, with NA where empty tiles were dropped — NOT a
+  # gap-free bounding box. Uses the packaged diagonal reach (area/bbox ~ 0.105).
+  aoi <- sf::st_read(
+    system.file("extdata", "example_aoi.gpkg", package = "drift"), quiet = TRUE
+  )
+  target_crs <- drift:::auto_utm_epsg(aoi)
+  aoi_t <- sf::st_transform(aoi, as.integer(gsub("EPSG:", "", target_crs)))
+  res <- 10
+  tile_size <- suppressMessages(drift:::tile_size_check(500, res))
+  tiles <- drift:::tile_grid(aoi_t, tile_size, res)
+
+  be <- sf::st_bbox(aoi_t)
+  n_full <- ceiling((be[["xmax"]] - be[["xmin"]]) / tile_size) *
+    ceiling((be[["ymax"]] - be[["ymin"]]) / tile_size)
+  expect_lt(length(tiles), n_full)               # empty tiles dropped (the mechanism)
+
+  # each kept tile filled with 1; merge -> rectangular bbox of kept tiles with
+  # NA in the skipped-tile gaps
+  merged <- drift:::mosaic_stacks(lapply(tiles, function(ext) {
+    terra::rast(xmin = ext$left, xmax = ext$right, ymin = ext$bottom,
+                ymax = ext$top, resolution = res, crs = target_crs, vals = 1)
+  }))
+  expect_true(anyNA(terra::values(merged)))      # gaps -> not a gap-free bbox
+})
+
 # Network end-to-end against the Planetary Computer. Opt-in only (env var), so
 # the default `devtools::test()` stays network-free per the repo convention.
 test_that("dft_stac_cube fetches an index stack end-to-end", {
@@ -138,3 +240,70 @@ test_that("dft_stac_cube fetches an index stack end-to-end", {
   expect_length(list.files(file.path(cache, "sentinel-2-l2a"),
                            pattern = "^cube_.*\\.tif$"), 1)
 })
+
+# Tiled read (#38) vs untiled, end-to-end. Opt-in. This is the only test that
+# exercises the real gdalcubes tiled read against live COGs — that a tile
+# sub-extent reads the same source pixels as the full bbox. The tiled and untiled
+# cubes are NOT co-lattice and can't be compared pixel-for-pixel: gdalcubes
+# enlarges the untiled bbox extent symmetrically to align with dx/dy (~0.5 px),
+# while the tiles are anchored at the bbox lower-left. Both are valid resamplings
+# of the same source, so equivalence is asserted via bilinear alignment
+# (correlation + bulk agreement) and grid-independent per-layer means. Verified
+# offline against saved 2021-07/08 cubes to have no tile seams (edge |diff| ==
+# interior) — the residual is the benign sub-pixel offset, not a seam. Thresholds
+# are measured on those cubes (cor 0.997, median |diff| 3.4e-3, per-layer mean
+# 6e-4) with headroom. A growing-season window is used for robust valid-pixel
+# coverage; the offset split under tiling is proven by the offline commutativity
+# oracle above, not here.
+test_that("dft_stac_cube tiled read reproduces the untiled cube 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
+  )
+  dtwin <- "2021-07-01/2021-08-31"
+  cache <- tempfile("drift_cube_tiled_")
+  dir.create(cache)
+
+  untiled <- dft_stac_cube(aoi, index = "kndvi", datetime = dtwin, dt = "P1M",
+                           cache_dir = cache)
+  tiled <- dft_stac_cube(aoi, index = "kndvi", datetime = dtwin, dt = "P1M",
+                         tile_size = 1000, cache_dir = cache)
+
+  expect_s4_class(tiled, "SpatRaster")
+  expect_equal(terra::nlyr(tiled), terra::nlyr(untiled))   # same monthly axis
+  expect_false(anyNA(terra::time(tiled)))                  # time set per layer
+  # tiled and untiled each cache one cube_<key>.tif, keyed apart (2 files total)
+  expect_length(list.files(file.path(cache, "sentinel-2-l2a"),
+                           pattern = "^cube_.*\\.tif$"), 2)
+
+  # the efficiency claim: for this diagonal reach the tiled read streams fewer
+  # tiles than the full grid (offline-computable, but assert it alongside the fetch)
+  target_crs <- drift:::auto_utm_epsg(aoi)
+  aoi_t <- sf::st_transform(aoi, as.integer(gsub("EPSG:", "", target_crs)))
+  ts <- suppressMessages(drift:::tile_size_check(1000, 10))
+  be <- sf::st_bbox(aoi_t)
+  n_full <- ceiling((be[["xmax"]] - be[["xmin"]]) / ts) *
+    ceiling((be[["ymax"]] - be[["ymin"]]) / ts)
+  expect_lt(length(drift:::tile_grid(aoi_t, ts, 10)), n_full)
+
+  # equivalence over common in-AOI cells: align tiled onto the untiled grid with
+  # BILINEAR (corrects the sub-pixel offset), then compare. A genuinely wrong
+  # tiled read (wrong per-tile offset, scrambled tiles, real seams) would drop the
+  # correlation and shift the means far past these bounds.
+  aligned <- terra::resample(tiled, untiled, method = "bilinear")
+  a <- terra::values(untiled)
+  b <- terra::values(aligned)
+  both <- !is.na(a) & !is.na(b)
+  expect_gt(sum(both), 0)                                        # overlap exists
+  expect_gt(stats::cor(a[both], b[both]), 0.98)                  # spatial pattern reproduced
+  expect_lt(stats::median(abs(a[both] - b[both])), 0.01)         # bulk agreement
+  # grid-independent: per-layer spatial-mean kNDVI agrees closely
+  layer_mean <- function(x) {
+    vapply(seq_len(terra::nlyr(x)),
+           function(i) mean(terra::values(x[[i]]), na.rm = TRUE), 0)
+  }
+  expect_lt(max(abs(layer_mean(tiled) - layer_mean(untiled))), 0.01)
+})