From 6da8fbdd4b5cf41ba434d0e71e86bc130cbcceb1 Mon Sep 17 00:00:00 2001 From: almac2022 Date: Wed, 30 Sep 2026 09:22:29 -0700 Subject: [PATCH 01/10] Initialize PWF baseline for #92 Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01PRhUJsuKABLfpBGktPoiBN --- planning/active/findings.md | 51 +++++++++++++++++++ planning/active/progress.md | 8 +++ planning/active/task_plan.md | 96 ++++++++++++++++++++++++++++++++++++ 3 files changed, 155 insertions(+) create mode 100644 planning/active/findings.md create mode 100644 planning/active/progress.md create mode 100644 planning/active/task_plan.md diff --git a/planning/active/findings.md b/planning/active/findings.md new file mode 100644 index 0000000..01e0f12 --- /dev/null +++ b/planning/active/findings.md @@ -0,0 +1,51 @@ +# Findings — dft_stac_composite(aggregation = "count") silently returns reflectance (#92) + +## Issue context + +## Problem + +`dft_stac_composite(aggregation = "count")` returns **reflectance, not a count**. It raises no error and no warning. + +`aggregation` goes straight into `gdalcubes::cube_view(aggregation = )` (`R/dft_stac_cube.R`, `stac_cube_assemble`), which supports only `"min"`, `"max"`, `"mean"`, `"median"` and `"first"` (gdalcubes 0.7.5, `?cube_view`). An unsupported value is not refused. What comes back looks like surface reflectance. + +This was found by floodplains#93 phase 2. That phase measures the clear-observation windows for reference imagery exactly as drift's docs describe: `aggregation = "count"`, `bands = "red"`, with the SCL mask applied before aggregation. + +## Repro (drift 0.19.0, gdalcubes 0.7.5) + +A 2 km square in the NECR floodplain, July, `res = 20`: + +```r +r <- dft_stac_composite(aoi, years = 2021, months = 7, bands = "red", + aggregation = "count", res = 20, crs = "EPSG:32610")[[1]] +summary(terra::values(r)) +#> Min 0.0077 Median 0.0331 Max 0.1835 # 2021 +#> Min 0.0078 Median 0.0411 Max 0.1802 NA 307 # same call, 2023 +``` + +The values are continuous and sit at typical vegetation red reflectance. The per-pixel maximum for the month should be at most the 17 items the query returned (6 distinct dates). The values are not integers times the band scale (1e-4) either, and 2023 carries no −0.1 offset shift. So this is not a count that was then scaled. `"count"` is simply not being honoured. + +## Asks + +1. **Validate `aggregation`** against what `cube_view` supports, in `dft_stac_composite()` and `dft_stac_cube()`, and refuse anything else. A wrong value should never return plausible reflectance. +2. **Provide clear-observation counts per window.** The obvious way is `dt = "P1D"` in the view, then `gdalcubes::reduce_time(cube, "count()")` over the window, with **no** scale or offset applied to the result. Either a new `aggregation = "count"` path or a separate function. The count is what a caller needs to choose composite windows, and the composite docs already suggest it. + +## Also affected + +floodplains#93's issue body tells callers to use `aggregation = "count"`. That call has to wait for this issue. +3. **Invalidate cached "count" composites.** The cache key hashes `aggregation`, so every `composite_.tif` written by 0.19.0 with `aggregation = "count"` holds reflectance under the key a fixed count would reuse. The fix has to change the key, for example by versioning it, or those stale files will be read back as counts. +4. **Count acquisitions, not items.** Where adjacent MGRS tiles overlap, one acquisition appears as two items. The validation square above spans 3 tiles (09UYV, 10UCE, 10UDE): 17 items but only 6 distinct dates. A count should de-duplicate same-day scenes (`dt = "P1D"` does this if the day is the time step). Also note that snow classes are in the default `mask_values`, so a spring or autumn "clear" count excludes snow as well as cloud. The docs should say so. + + +## Plan-mode exploration (2026-09-30) + +- `gdalcubes::cube_view()` (0.7.5) validates `aggregation` only as `is.character` + length 1; `?cube_view` documents "min", "max", "mean", "median", "first". +- Three callers pass it through unvalidated: `stac_cube_assemble()` (composite + cube) and `fetch_extent_to()` (fetch). +- `stac_cube_assemble()` coalesces the offset split with `terra::cover(pre, post)`. For a count that would discard the post side wherever the pre count is 0, because 0 is not NA. So count must not split. +- `cube_check_nonempty()` tests for `notNA`. An all-zero count passes it, which means a failed read looks the same as "all cloudy". +- The `reduce_time()` closure gotcha (inst/notes/gdalcubes-pc-gotchas.md) applies to R callbacks. The string reducer `"count(B04)"` is built-in C++. +- The #79 BULK floodplain-wide composite: 3 h 16 min, then #88 (the mask of the merged mosaic fails), open. 100 chips: 84.1 min, 0.48 GiB peak. + +## Errors Encountered + +| Error | Resolution | +|-------|------------| diff --git a/planning/active/progress.md b/planning/active/progress.md new file mode 100644 index 0000000..93be213 --- /dev/null +++ b/planning/active/progress.md @@ -0,0 +1,8 @@ +# Progress — dft_stac_composite(aggregation = "count") (#92) + +## Session 2026-09-30 + +- Plan-mode exploration; phases approved by user (count lives in `aggregation = "count"`; count cache keyed under its own tag, `count_.tif`) +- Created branch `92-dft-stac-composite-aggregation-count-si` off main +- Scaffolded PWF baseline from issue #92 with approved phases +- Next: start Phase 1 diff --git a/planning/active/task_plan.md b/planning/active/task_plan.md new file mode 100644 index 0000000..90cb465 --- /dev/null +++ b/planning/active/task_plan.md @@ -0,0 +1,96 @@ +# Task: dft_stac_composite(aggregation = "count") silently returns reflectance, not a clear-observation count (#92) + +## Problem + +`dft_stac_composite(aggregation = "count")` returns **reflectance, not a count**. It raises no error and no warning. + +`aggregation` goes straight into `gdalcubes::cube_view(aggregation = )` (`R/dft_stac_cube.R`, `stac_cube_assemble`), which supports only `"min"`, `"max"`, `"mean"`, `"median"` and `"first"` (gdalcubes 0.7.5, `?cube_view`). An unsupported value is not refused. What comes back looks like surface reflectance. + +This was found by floodplains#93 phase 2. That phase measures the clear-observation windows for reference imagery exactly as drift's docs describe: `aggregation = "count"`, `bands = "red"`, with the SCL mask applied before aggregation. + + +Decisions taken at the gate: + +- **The count is `aggregation = "count"` inside `dft_stac_composite()`**, with no new + export. The count path hashes under **its own tag** and writes `count_.tif`, so + the 0.19.x reflectance-as-count files are orphaned rather than served, and every genuine + composite cache stays valid. + +## Phase 1: Refuse an aggregation gdalcubes does not honour +- [ ] Add `aggregation_check(aggregation, allowed)`, a shared validator in + `R/dft_stac_cube.R`. It requires a single non-NA string in `allowed` and aborts + with a classed error that names the valid set. +- [ ] `.cube_view_aggregations <- c("min", "max", "mean", "median", "first")`, one + constant, commented with its gdalcubes `?cube_view` source. +- [ ] Call the validator in `dft_stac_cube()` and `dft_stac_fetch()` (allowed = the gdalcubes + set) and in `dft_stac_composite()` (allowed = the gdalcubes set plus `"count"`). + Each check runs before any network call or cache lookup. + `dft_stac_fetch()` is the third caller that shares the harness, so it gets the check too. +- [ ] Tests, offline: the check fires in all three functions for `"count"` (except the + composite), `"sum"`, `NA`, `c("median", "mean")` and a non-character value. Confirm + the guard runs before any network call by mocking `stac_cube_items` / the STAC + search to `stop()`. Restore the defect (skip the check) and confirm the tests go red. + +## Phase 2: Clear-observation count path in `dft_stac_composite()` +- [ ] A count window uses `dt = "P1D"` in the `cube_view`, with a NaN-ignoring day-level + aggregation. Same-day items from overlapping MGRS tiles therefore collapse to one + acquisition (ask 4). The pixel function is `gdalcubes::reduce_time(cube, + "count()")` per band, renamed to the band role. It uses the built-in string + reducer, not an R callback, so the gotchas-note closure trap does not apply. No scale + or offset is applied. +- [ ] **No offset split for count.** Hand `stac_cube_assemble()` a `fetched` with `is_pre` + all FALSE, and skip `composite_offset_check()`. Otherwise `terra::cover(pre, post)` + would take the pre-side count wherever it is non-NA, and a count of 0 is non-NA, so + the post side would silently never be counted. A window straddling 2022-01-25 is + valid for a count. +- [ ] Separate cache family: pass a `"count"` tag into `stac_composite_cache_key()` (or + add a family argument) and use the file prefix `count_`. Write with + `datatype = "INT2U"`, so the COG stores integers. +- [ ] Empty-read guard: a window whose count is 0 everywhere is the gdalcubes failed-read + mode wearing the face of "all cloudy". Treat it the way a composite with no clear + pixels is treated: warn and drop the year (`drift_empty_cube`). + `cube_check_nonempty()` cannot see it, because 0 is not NA. +- [ ] Offline tests. Cover the key: count ≠ median over the same inputs, count stays + stable, and composite keys are unchanged, pinned against a literal pre-change key. + Cover the file prefix. Cover the count pixel function on a local + `create_image_collection()` fixture of dated GeoTIFFs: two same-day items count + once, a masked day is not counted, and the result is an integer with no scale. + Cover the straddling-window count, which is not refused. + Restore each defect and confirm its test goes red. + +## Phase 3: Docs +- [ ] `@param aggregation` in `dft_stac_composite()` lists the valid values and describes + `"count"`: distinct clear days per pixel, per band, integer, no scale; the + `cloud_cover_max` pre-filter applies, so the count is of scenes that pass it. + `dft_stac_cube()` and `dft_stac_fetch()` list their valid values. +- [ ] `@param mask_values`: the default SCL classes include **snow**, so a spring or autumn + "clear" count excludes snow as well as cloud (ask 4). +- [ ] Add a count example to `@examples` (inside the existing `\dontrun{}`, since it needs + the network). Run `devtools::document()` and `pkgdown::check_pkgdown()`. + +## Phase 4: Live verification and scale +- [ ] Network test (`DRIFT_TEST_NETWORK`): on the packaged AOI, the count is an integer, + ≥ 0, and ≤ the number of distinct item dates in the window. It must not look like + reflectance, so its maximum must be ≥ 1 and must be a whole number. +- [ ] Reproduce the issue's check on a small square: the old call gave continuous values + around 0.03, and the new max must be ≤ the distinct dates (6 in the issue's July 2021 + square). +- [ ] Scale: count chips on BULK (about 20 points, 300 m buffers, one month) under + `/usr/bin/time -l`, recording peak RSS and wall clock, beside #79's 100-chip median + run (0.48 GiB, 84.1 min). The floodplain-wide read is out of scope while #88 is open. + Numbers go in the PR body and the archive README. + +## Phase 5: Release notes +- [ ] NEWS entry covering: + - the silent reflectance, and which callers now refuse what + - the new count and its semantics + - the orphaned 0.19.x `composite_*.tif` files written under `"count"`, which are never + read and can be reclaimed via `dft_cache_info()` / `dft_cache_clear()` guidance + - snow in the default mask + +## Validation + +- [ ] Tests pass (`devtools::test()`), `lintr::lint_package()` clean +- [ ] `/code-check` clean on each commit +- [ ] PWF checkboxes match landed work +- [ ] `/planning-archive` on completion, then `/gh-pr-push` (PR body: `Relates to NewGraphEnvironment/sred-2025-2026#16`) From 727c5d95d67f5afbb05d782e4cbc629c30b7e155 Mon Sep 17 00:00:00 2001 From: almac2022 Date: Wed, 30 Sep 2026 09:34:39 -0700 Subject: [PATCH 02/10] Refuse an aggregation gdalcubes would not honour (#92) gdalcubes::cube_view() maps an unknown aggregation to "none" without an error, so dft_stac_composite(aggregation = "count") returned reflectance. dft_stac_cube(), dft_stac_fetch() and dft_stac_composite() now validate against the six values measured to survive a cube_view() round trip (min, max, mean, median, first, last), case-insensitively, before any network call or cache lookup. count_values/count_images are held back: they count items and every caller would scale them as reflectance. Relates to #92 Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01PRhUJsuKABLfpBGktPoiBN --- R/dft_stac_composite.R | 1 + R/dft_stac_cube.R | 37 +++++++++++ R/dft_stac_fetch.R | 1 + planning/active/findings.md | 28 +++++++++ planning/active/review-p1-round1.md | 58 ++++++++++++++++++ planning/active/review-p1-round2.md | 35 +++++++++++ planning/active/task_plan.md | 78 +++++++++++++++--------- tests/testthat/test-dft_stac_composite.R | 14 +++++ tests/testthat/test-dft_stac_cube.R | 60 ++++++++++++++++++ 9 files changed, 282 insertions(+), 30 deletions(-) create mode 100644 planning/active/review-p1-round1.md create mode 100644 planning/active/review-p1-round2.md diff --git a/R/dft_stac_composite.R b/R/dft_stac_composite.R index 3fed511..8d897ca 100644 --- a/R/dft_stac_composite.R +++ b/R/dft_stac_composite.R @@ -150,6 +150,7 @@ dft_stac_composite <- function(aoi, "i" = "Use {.fn dft_stac_fetch} for categorical rasters." )) } + aggregation <- aggregation_check(aggregation) years <- composite_years_check(years) months <- composite_months_check(months) band_assets <- composite_band_assets(bands, cfg$roles) diff --git a/R/dft_stac_cube.R b/R/dft_stac_cube.R index 09e4dcb..f952839 100644 --- a/R/dft_stac_cube.R +++ b/R/dft_stac_cube.R @@ -163,6 +163,7 @@ dft_stac_cube <- function(aoi, force = FALSE, sign_fn = rstac::sign_planetary_computer()) { check_gdalcubes("to fetch STAC cubes") + aggregation <- aggregation_check(aggregation) # gdalcubes worker count and GDAL /vsicurl tuning, both restored on exit so # the caller's session is untouched (see stac_cube_session()). @@ -295,6 +296,42 @@ dft_stac_cube <- function(aoi, } +# The aggregations drift passes to gdalcubes::cube_view(). gdalcubes 0.7.5 +# lower-cases the string and maps anything it does not know to "none" without +# error, which is how "count" came back as reflectance (#92). Measured by +# round-tripping each value through cube_view()$aggregation: the six below +# survive. It also honours "count_values" and "count_images", which are held back +# on purpose: they count items, so overlapping MGRS tiles double-count one +# acquisition, and every caller here would then apply a reflectance scale or an +# index expression to a count. "none" is the silent fallback itself. +.cube_view_aggregations <- c("min", "max", "mean", "median", "first", "last") + +#' Refuse an `aggregation` outside `allowed`, before any network call +#' +#' Shared by [dft_stac_fetch()], [dft_stac_cube()] and [dft_stac_composite()], +#' the three callers of `gdalcubes::cube_view()`. An unsupported value is not +#' refused downstream, so it has to be refused here (#92). Matching ignores case, +#' as gdalcubes does, and the value is returned lower-cased so `"Median"` and +#' `"median"` share a cache key. +#' @noRd +aggregation_check <- function(aggregation, allowed = .cube_view_aggregations) { + if (!is.character(aggregation) || length(aggregation) != 1L || + is.na(aggregation) || !tolower(aggregation) %in% allowed) { + cli::cli_abort(c( + "{.arg aggregation} must be one of {.or {.val {allowed}}}.", + "x" = if (is.character(aggregation) && length(aggregation) == 1L && + !is.na(aggregation)) { + "Got {.val {aggregation}}. gdalcubes reads a value it does not know as \\ + {.val none} without an error, so drift passes only these." + } else { + "Got {.obj_type_friendly {aggregation}}." + } + ), class = "drift_bad_aggregation") + } + tolower(aggregation) +} + + #' Set the gdalcubes worker count and GDAL cloud-read tuning for a cube read #' #' Shared by [dft_stac_cube()] and [dft_stac_composite()]. Returns a function diff --git a/R/dft_stac_fetch.R b/R/dft_stac_fetch.R index bf4147b..5a0b3ea 100644 --- a/R/dft_stac_fetch.R +++ b/R/dft_stac_fetch.R @@ -99,6 +99,7 @@ dft_stac_fetch <- function(aoi, force = FALSE, sign_fn = rstac::sign_planetary_computer()) { check_gdalcubes("to fetch STAC rasters") + aggregation <- aggregation_check(aggregation) # 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 diff --git a/planning/active/findings.md b/planning/active/findings.md index 01e0f12..e38e5e0 100644 --- a/planning/active/findings.md +++ b/planning/active/findings.md @@ -45,7 +45,35 @@ floodplains#93's issue body tells callers to use `aggregation = "count"`. That c - The `reduce_time()` closure gotcha (inst/notes/gdalcubes-pc-gotchas.md) applies to R callbacks. The string reducer `"count(B04)"` is built-in C++. - The #79 BULK floodplain-wide composite: 3 h 16 min, then #88 (the mask of the merged mosaic fails), open. 100 chips: 84.1 min, 0.48 GiB peak. +## Count semantics, measured offline (2026-09-30, gdalcubes 0.7.5) + +Probe: 6 local items on a 4x4 grid, EPSG:32609, July 2021. The inputs: + +- 07-03: tiles A and B, both clear +- 07-08: tile A, cloud (SCL 9) +- 07-13: tile A, clear +- 07-18: tile A cloud (SCL 8), tile B clear + +`cube_view(dt = "P1D")` + `image_mask("SCL", c(3, 8, 9, 10, 11))` + `reduce_time(select_bands(cube, "B04"), "count(B04)")`: + +- The count is **3** under `aggregation = "first"`, `"max"` and `"median"` alike. Same-day items collapse to one day, and a masked item does not blank a clear same-day item, because the aggregators skip NaN. +- The output band is named `B04_count`. `rename_bands(rc, B04_count = "red")` renames it. +- **A pixel with zero clear days is NaN, not 0.** That holds both when every day is masked and when the window has no images. Consequences: + - `cube_check_nonempty()` already catches an all-empty count, so the planned extra zero-guard is unnecessary. It is dropped. + - The docs must say that NA means no clear observation (or no coverage), not that a 0 is stored. + - The offset split is still wrong for a count. `terra::cover(pre, post)` keeps the pre-side count wherever it is non-NA, so a pixel with 2 pre days and 3 post days would read 2. + +## Plan review (Plan agent, 2026-09-30) — what was acted on + +- **Blocker: `"count"` must never reach `cube_view`.** The count path uses `"first"` and P1D, and `stac_cube_assemble()` also checks. +- **The "0 everywhere = failed read" premise was wrong.** Measured (probe3, 128x128 fixture): gdalcubes gives an all-NaN chunk empty (NA) and a zero-clear pixel inside a chunk that has clear pixels elsewhere a 0. The same data gave 12,544 zeros under 256 px chunks, and 4,352 zeros plus 8,192 NaN under 64 px chunks. Chunk size follows `parallel`, so the raw output depended on a knob documented as cost-only. **Rule: 0 -> NA**, after which both chunkings are identical. NA was chosen over 0 because a failed chunk read (NaN) cannot be told apart from "all cloudy", so it fails toward NA. +- The first probe's "all masked -> NaN" case used `t0 == t1`. reduce_time passes a single time step through unchanged (reduce_time.cpp:533), so that case proved nothing. Superseded by probe3. +- Accepted: key on `dt = "P1D"`; NEAREST overviews; "count" labels; docs beyond `@param`; NEWS reclaim wording; the note in the gotchas file. +- Not acted on: the misleading "offset split" message from `stac_cube_items()` under count (cosmetic). +- `create_image_collection(one_band_per_file = FALSE)` on 2-band files segfaulted a gdalcubes worker (0.7.5). Separate files per band plus a format JSON work, so the tests use those. + ## Errors Encountered | Error | Resolution | |-------|------------| +| Round 1: `.cube_view_aggregations` pinned to `?cube_view`, which refuses `"last"` (honoured) and case variants | Probed `cube_view()$aggregation`: honours min/max/mean/median/first/last, count_values, count_images; lower-cases; else `"none"`. Set is now the six, case-insensitive; test pins behaviour via round trip, not the Rd | diff --git a/planning/active/review-p1-round1.md b/planning/active/review-p1-round1.md new file mode 100644 index 0000000..1936a7e --- /dev/null +++ b/planning/active/review-p1-round1.md @@ -0,0 +1,58 @@ +# Review: #92 Phase 1 (aggregation validation), round 1 + +## Findings + +- **[bug]** R/dft_stac_cube.R:299-303 (`.cube_view_aggregations`) — the set is the one gdalcubes + *documents*, not the one it *honours*, and the comment says "honours". gdalcubes 0.7.5's + C++ parser (`src/gdalcubes/src/view.cpp`, `aggregation::from_string`, unchanged since + 2024-03, and master is 0.7.5) accepts `none`, `min`, `max`, `mean`, `median`, `first`, + **`last`**, **`count_images`**, **`count_values`**, lower-casing the input first. Anything + else falls through to `AGG_NONE`, which is why `"count"` returned reflectance. Measured + offline by passing each string through `gdalcubes::cube_view()` and reading back + `$aggregation`, which comes from the C++ side: + + ``` + median -> median Median -> median last -> last + count_values -> count_values count_images -> count_images + count -> none sum -> none none -> none + ``` + + This has three effects. + 1. **Regression.** `aggregation = "last"`, a real mode, worked before this commit and is + now refused in all three functions. It is also the one most useful to + `dft_stac_fetch()`'s categorical path, as the counterpart of its default `"first"`. The + same goes for `"Median"` and the other case variants. That one hardly matters, but the + test asserts the refusal as intended behaviour. + 2. **The pin test checks the proxy.** The test "the cube_view aggregation set is the one + gdalcubes documents" (tests/testthat/test-dft_stac_cube.R, new block) parses + `?cube_view`, so it agrees with the constant and can never see the gap above. A test of + what gdalcubes actually does is the round-trip above. Do the check with + `cube_view(..., aggregation = a)$aggregation != "none"`, which is offline and costs + milliseconds. + 3. **Bearing on Phase 2.** gdalcubes already has a native per-pixel count in the view, + `aggregation = "count_values"`, and a per-window image count, `"count_images"`. Phase 2's + plan of `dt = "P1D"` + `reduce_time("count(...)")` may be the right design anyway, + because `count_values` counts items, so it double-counts where MGRS tiles overlap (the + issue's ask 4). Still, decide that on purpose rather than because the option was never + seen. Also: if `count_*` or `last` is added to the allowed set, `count_*` must bypass + the scale/offset in the cube and composite, or the #92 class of bug comes back as + "scaled counts". + + Suggested fix, for the parent to decide: list the set gdalcubes honours, add `"last"` at + least, and hold back the `count_*` pair (and `"none"`) until Phase 2 handles scale. Either + way, correct the comment and replace or add to the Rd-parse test with the round-trip + assertion, so the constant is pinned to behaviour. + +## Checked and fine + +- The guard runs in all three callers of `cube_view()` before any network, cache or + directory side effect. In the composite it follows only `dft_stac_config(source)`, which + is offline. `dft_stac_fetch()` has no path that avoids `cube_view()`, so the check applies + to all of it. +- The guard's short-circuit order is safe for `NULL`, `character(0)`, `NA` and length > 1. + The cli message interpolates in the validator's own frame. `{.obj_type_friendly}` is + already used in the package. +- The mocked bindings `stac_cube_items` and `stac_items_paged` exist, and testthat is + pinned >= 3.2.0. The Rd-parse regex fails loudly, not toward pass, when it does not match + (`[[1]]` on an empty result errors). It gives the same result under the C and UTF-8 + locales. diff --git a/planning/active/review-p1-round2.md b/planning/active/review-p1-round2.md new file mode 100644 index 0000000..b84a26b --- /dev/null +++ b/planning/active/review-p1-round2.md @@ -0,0 +1,35 @@ +# Review: #92 Phase 1 (aggregation validation), round 2 + +## Clean + +No bugs, security issues or fragile paths found in the round-1 fix. + +## Checked and fine + +- **Lower-casing and cache keys.** `aggregation_check()` returns `tolower(aggregation)`, which + is the value `stac_cache_key()` (R/dft_stac_fetch.R:183), `stac_cube_cache_key()` + (R/dft_stac_cube.R:221) and `stac_composite_cache_key()` (R/dft_stac_composite.R:189) hash. + For every lower-case input, including all three defaults (`"first"`, `"median"`), it does + nothing, so no existing cache entry is invalidated. The only keys that move are ones a + caller built with a mixed-case value such as `"Median"`. Those now map to the `"median"` + key, and gdalcubes lower-cased the value itself, so both produced identical content. The + old entries are orphaned but never read as wrong data. The existing key tests + (test-dft_stac_fetch.R:88, test-dft_stac_cube.R:80, test-dft_stac_composite.R:60) all use + lower-case values and call the key functions directly, so the change does not reach them. + `aggregation` is not written into metadata or returned attributes anywhere else, so + lower-casing changes nothing downstream. +- **No bypass.** All three `gdalcubes::cube_view()` calls in `R/` are internal: + `fetch_extent_to()` (dft_stac_fetch.R:608) and `stac_cube_assemble()` (dft_stac_cube.R:593). + Their only callers are `dft_stac_fetch()` (lines 215/221), `dft_stac_cube()` (line 258) + and `dft_stac_composite()` (line 211), and each of those reassigns `aggregation` from the + check before building the key or calling the helper. The remaining `cube_view()` calls are + in `data-raw/` scripts with hard-coded lower-case values, outside the package API. +- **Order.** In all three functions the check runs before any cache-directory creation, + STAC query or cache lookup. In the composite, only `check_gdalcubes()` and the offline + `dft_stac_config(source)` run before it. +- **Messages.** Probed: `NA` gives "Got `NA`." and `"Count"` gives "Got \"Count\", which + gdalcubes would not honour.", with the valid set listed. +- **Tests.** Both changed test files pass when run in a copy with `NOT_CRAN=true`: + test-dft_stac_cube.R gives 141 pass / 0 fail / 3 network skips, and + test-dft_stac_composite.R gives 61 pass / 0 fail / 2 network skips. The round-trip test + runs offline and pins the constant to what gdalcubes does, not to its Rd. diff --git a/planning/active/task_plan.md b/planning/active/task_plan.md index 90cb465..17749a5 100644 --- a/planning/active/task_plan.md +++ b/planning/active/task_plan.md @@ -17,48 +17,65 @@ Decisions taken at the gate: composite cache stays valid. ## Phase 1: Refuse an aggregation gdalcubes does not honour -- [ ] Add `aggregation_check(aggregation, allowed)`, a shared validator in +- [x] Add `aggregation_check(aggregation, allowed)`, a shared validator in `R/dft_stac_cube.R`. It requires a single non-NA string in `allowed` and aborts with a classed error that names the valid set. -- [ ] `.cube_view_aggregations <- c("min", "max", "mean", "median", "first")`, one +- [x] `.cube_view_aggregations <- c("min", "max", "mean", "median", "first")`, one constant, commented with its gdalcubes `?cube_view` source. -- [ ] Call the validator in `dft_stac_cube()` and `dft_stac_fetch()` (allowed = the gdalcubes +- [x] Call the validator in `dft_stac_cube()` and `dft_stac_fetch()` (allowed = the gdalcubes set) and in `dft_stac_composite()` (allowed = the gdalcubes set plus `"count"`). Each check runs before any network call or cache lookup. `dft_stac_fetch()` is the third caller that shares the harness, so it gets the check too. -- [ ] Tests, offline: the check fires in all three functions for `"count"` (except the +- [x] Tests, offline: the check fires in all three functions for `"count"` (except the composite), `"sum"`, `NA`, `c("median", "mean")` and a non-character value. Confirm the guard runs before any network call by mocking `stac_cube_items` / the STAC search to `stop()`. Restore the defect (skip the check) and confirm the tests go red. ## Phase 2: Clear-observation count path in `dft_stac_composite()` -- [ ] A count window uses `dt = "P1D"` in the `cube_view`, with a NaN-ignoring day-level - aggregation. Same-day items from overlapping MGRS tiles therefore collapse to one - acquisition (ask 4). The pixel function is `gdalcubes::reduce_time(cube, - "count()")` per band, renamed to the band role. It uses the built-in string - reducer, not an R callback, so the gotchas-note closure trap does not apply. No scale - or offset is applied. -- [ ] **No offset split for count.** Hand `stac_cube_assemble()` a `fetched` with `is_pre` - all FALSE, and skip `composite_offset_check()`. Otherwise `terra::cover(pre, post)` - would take the pre-side count wherever it is non-NA, and a count of 0 is non-NA, so - the post side would silently never be counted. A window straddling 2022-01-25 is - valid for a count. -- [ ] Separate cache family: pass a `"count"` tag into `stac_composite_cache_key()` (or - add a family argument) and use the file prefix `count_`. Write with - `datatype = "INT2U"`, so the COG stores integers. -- [ ] Empty-read guard: a window whose count is 0 everywhere is the gdalcubes failed-read - mode wearing the face of "all cloudy". Treat it the way a composite with no clear - pixels is treated: warn and drop the year (`drift_empty_cube`). - `cube_check_nonempty()` cannot see it, because 0 is not NA. -- [ ] Offline tests. Cover the key: count ≠ median over the same inputs, count stays - stable, and composite keys are unchanged, pinned against a literal pre-change key. - Cover the file prefix. Cover the count pixel function on a local - `create_image_collection()` fixture of dated GeoTIFFs: two same-day items count - once, a masked day is not counted, and the result is an integer with no scale. - Cover the straddling-window count, which is not refused. - Restore each defect and confirm its test goes red. + +Revised after the plan review (see findings.md, "Plan review" and "Count semantics"). + +- [ ] The count path calls `stac_cube_assemble()` with `dt = "P1D"` and a day + aggregation of `"first"`, never `"count"`. `"count"` must never reach `cube_view`, + because it maps to `AGG_NONE` (the #92 bug). As a last guard, + `aggregation_check()` also runs inside `stac_cube_assemble()`. +- [ ] Add a helper, `composite_count_cube(cube, band_assets, bands)`, that returns + `reduce_time(cube, paste0("count(", band_assets, ")"), names = bands)`. It + applies no scale or offset, and `names = bands` satisfies `composite_layers_order()`. +- [ ] **A pixel with zero clear days is NA, always** (`count_zero_na()`, `0 -> NA`). + gdalcubes returns 0 or NaN for it depending on chunk layout, and chunk size + follows `parallel`. This rule is therefore what keeps `parallel` a cost-only + knob. It fails toward NA because a failed chunk read cannot be told apart from + "all cloudy". An all-NA result is caught by `cube_check_nonempty()` and the year + is dropped (`drift_empty_cube`). +- [ ] **No offset split for count.** Set `is_pre` to all FALSE and skip + `composite_offset_check()`, because `terra::cover(pre, post)` would drop the post + side. A window that straddles 2022-01-25 is valid for a count. +- [ ] Give the count its own cache family: the `"count"` tag in + `stac_composite_cache_key()` (a family argument whose default is `"composite"`, + so the median key is byte-identical, pinned to `03ee8ecc66b832a8`), the key + hashing `dt = "P1D"`, and the file prefix `count_`. The count key must not + equal the stale 0.19 key `08e0c5510e8ae297`. Write the file as INT2U with + `OVERVIEW_RESAMPLING=NEAREST`, and use `"count"` as the label in the + cache-read and hit messages. +- [ ] Offline tests: + - [ ] Test `composite_count_cube()` on a local `create_image_collection()` + fixture of separate B04/SCL files. Cover: a clear same-day item not blanked + by a masked one, a count of distinct days rather than items, integer output + with no scale, and names equal to the roles. + - [ ] Test that the result is chunking-invariant after `count_zero_na()` + (`chunking = c(16, 64, 64)` against `c(16, 256, 256)`). + - [ ] Build `dft_stac_composite(aggregation = "count")` offline with a mocked + `stac_cube_assemble` that captures its arguments. Assert `"first"` and `"P1D"`, + that `is_pre` is all FALSE with no refusal for a straddling window, the + `count_.tif` name, INT2U, and that zeros become NA. + - [ ] Restore each defect and confirm the matching test goes red. +- [ ] Add the `AGG_NONE` root cause and the chunk-dependent 0/NaN behaviour to + `inst/notes/gdalcubes-pc-gotchas.md`. ## Phase 3: Docs +- [ ] Description, `@return` and the Caching section carry the count exception + (integer counts, not reflectance; `count_.tif`). - [ ] `@param aggregation` in `dft_stac_composite()` lists the valid values and describes `"count"`: distinct clear days per pixel, per band, integer, no scale; the `cloud_cover_max` pre-filter applies, so the count is of scenes that pass it. @@ -85,7 +102,8 @@ Decisions taken at the gate: - the silent reflectance, and which callers now refuse what - the new count and its semantics - the orphaned 0.19.x `composite_*.tif` files written under `"count"`, which are never - read and can be reclaimed via `dft_cache_info()` / `dft_cache_clear()` guidance + read. They sit in the current `v2` scheme and cannot be told apart from real + composites, so `scheme = "superseded"` does not reclaim them. State this plainly. - snow in the default mask ## Validation diff --git a/tests/testthat/test-dft_stac_composite.R b/tests/testthat/test-dft_stac_composite.R index 356aafb..157541a 100644 --- a/tests/testthat/test-dft_stac_composite.R +++ b/tests/testthat/test-dft_stac_composite.R @@ -295,3 +295,17 @@ test_that("a built composite leaves only its COG in the cache, no sidecar", { tolerance = 1e-6) expect_equal(terra::time(out[[1]])[1], as.Date("2023-06-01")) }) + +test_that("dft_stac_composite refuses an aggregation gdalcubes would not honour (#92)", { + skip_if_not_installed("gdalcubes") + testthat::local_mocked_bindings( + stac_cube_items = function(...) stop("reached the network") + ) + cache <- withr::local_tempdir() + for (b in list("sum", "mode", NA_character_, c("median", "mean"), 1)) { + expect_error(dft_stac_composite(aoi_pkg(), years = 2021, aggregation = b, + cache_dir = cache), + class = "drift_bad_aggregation") + } + expect_length(list.files(cache, recursive = TRUE), 0L) +}) diff --git a/tests/testthat/test-dft_stac_cube.R b/tests/testthat/test-dft_stac_cube.R index 6557579..37a65d2 100644 --- a/tests/testthat/test-dft_stac_cube.R +++ b/tests/testthat/test-dft_stac_cube.R @@ -754,3 +754,63 @@ test_that("stac_features_resign re-signs features before a read, stubs pass thro list(feat("x"))) expect_identical(drift:::stac_features_resign(fetched, list()), list()) }) + +test_that("aggregation_check refuses anything outside the set drift passes to cube_view (#92)", { + for (ok in c("min", "max", "mean", "median", "first", "last")) { + expect_identical(drift:::aggregation_check(ok), ok) + } + # case-insensitive, as gdalcubes is, and returned lower-cased for the cache key + expect_identical(drift:::aggregation_check("Median"), "median") + expect_identical(drift:::aggregation_check("LAST"), "last") + bad <- list("count", "sum", "none", "count_values", "count_images", "", + NA_character_, NA, c("median", "mean"), character(0), 1, NULL) + for (b in bad) { + expect_error(drift:::aggregation_check(b), class = "drift_bad_aggregation") + } + # the message names the valid set and the value it refused + expect_error(drift:::aggregation_check("count"), "median") + expect_error(drift:::aggregation_check("count"), "count") + # an extended set admits its extra member only + expect_identical(drift:::aggregation_check("count", c("median", "count")), + "count") +}) + +test_that("every allowed aggregation is one gdalcubes::cube_view actually honours", { + # Behaviour, not documentation: gdalcubes maps an unknown aggregation to + # "none" silently, so each allowed value must survive the round trip, and + # "count" must not (it is the #92 value that came back as reflectance). + skip_if_not_installed("gdalcubes") + honoured <- function(a) { + gdalcubes::cube_view( + srs = "EPSG:32609", + extent = list(left = 0, right = 40, bottom = 0, top = 40, + t0 = "2021-07-01", t1 = "2021-07-31"), + dx = 10, dy = 10, dt = "P1M", aggregation = a + )$aggregation + } + for (a in drift:::.cube_view_aggregations) expect_identical(honoured(a), a) + expect_identical(honoured("count"), "none") +}) + +test_that("dft_stac_cube and dft_stac_fetch refuse a bad aggregation before any network call (#92)", { + skip_if_not_installed("gdalcubes") + aoi <- sf::st_read( + system.file("extdata", "example_aoi.gpkg", package = "drift"), + quiet = TRUE + ) + cache <- withr::local_tempdir() + # reaching either query means validation ran too late + testthat::local_mocked_bindings( + stac_cube_items = function(...) stop("reached the network"), + stac_items_paged = function(...) stop("reached the network") + ) + for (b in list("count", "sum", NA_character_, c("median", "mean"))) { + expect_error(dft_stac_cube(aoi, aggregation = b, cache_dir = cache), + class = "drift_bad_aggregation") + expect_error(dft_stac_fetch(aoi, source = "io-lulc", years = 2017, + aggregation = b, cache_dir = cache), + class = "drift_bad_aggregation") + } + # nothing was written for a refused call + expect_length(list.files(cache, recursive = TRUE), 0L) +}) From ba0b2039924919ec80dac9e78a4b438e375d99a9 Mon Sep 17 00:00:00 2001 From: almac2022 Date: Wed, 30 Sep 2026 09:55:26 -0700 Subject: [PATCH 03/10] Count clear days with dft_stac_composite(aggregation = "count") (#92) A count reads the window one day per time step (dt = "P1D", day aggregation "first") and reduces with reduce_time("count()"), so overlapping MGRS tiles count one acquisition once and a masked tile does not blank a clear one. No scale, no offset, no offset split (terra::cover would drop the post side). A pixel with no clear day is NA: gdalcubes returns 0 or NaN depending on chunk layout, which follows `parallel`. The count keys under its own tag and writes count_.tif (INT2U, NEAREST overviews), so 0.19.x reflectance-as-count files are never read and every composite key is unchanged. stac_cube_assemble() refuses an unhonoured aggregation as a last guard. Relates to #92 Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01PRhUJsuKABLfpBGktPoiBN --- R/dft_stac_composite.R | 110 ++++++++-- R/dft_stac_cube.R | 3 + inst/notes/gdalcubes-pc-gotchas.md | 13 ++ planning/active/findings.md | 67 ++++++ planning/active/progress.md | 4 + planning/active/review-p1-round3.md | 56 +++++ planning/active/review-p2-round1.md | 22 ++ planning/active/review-p2-round2.md | 63 ++++++ planning/active/review-p2-round3.md | 50 +++++ planning/active/task_plan.md | 22 +- tests/testthat/test-dft_stac_composite.R | 267 +++++++++++++++++++++++ 11 files changed, 645 insertions(+), 32 deletions(-) create mode 100644 planning/active/review-p1-round3.md create mode 100644 planning/active/review-p2-round1.md create mode 100644 planning/active/review-p2-round2.md create mode 100644 planning/active/review-p2-round3.md diff --git a/R/dft_stac_composite.R b/R/dft_stac_composite.R index 8d897ca..e7e5fe9 100644 --- a/R/dft_stac_composite.R +++ b/R/dft_stac_composite.R @@ -150,7 +150,10 @@ dft_stac_composite <- function(aoi, "i" = "Use {.fn dft_stac_fetch} for categorical rasters." )) } - aggregation <- aggregation_check(aggregation) + aggregation <- aggregation_check(aggregation, + c(.cube_view_aggregations, "count")) + is_count <- identical(aggregation, "count") + family <- if (is_count) "count" else "composite" years <- composite_years_check(years) months <- composite_months_check(months) band_assets <- composite_band_assets(bands, cfg$roles) @@ -176,30 +179,39 @@ dft_stac_composite <- function(aoi, # Every band scaled to reflectance with the offset for the side of the # 2022-01-25 boundary the window's scenes fall on (a window straddling it is - # refused by composite_offset_check()). - pixel_fn <- function(cube, offset_use) { - exprs <- vapply(band_assets, scale_token, character(1), - scale = scale, offset = offset_use, USE.NAMES = FALSE) - gdalcubes::apply_pixel(cube, exprs, names = bands) + # refused by composite_offset_check()). A count applies neither: it counts + # clear days, and the offset only moves values, not whether they are masked. + pixel_fn <- if (is_count) { + function(cube, offset_use) composite_count_cube(cube, band_assets, bands) + } else { + function(cube, offset_use) { + exprs <- vapply(band_assets, scale_token, character(1), + scale = scale, offset = offset_use, USE.NAMES = FALSE) + gdalcubes::apply_pixel(cube, exprs, names = bands) + } } out <- lapply(years, function(year) { w <- composite_window(year, months) label <- composite_label(year, months) + # A count is read one day per time step (see composite_count_cube()), so + # that is the step it reads with and the step it keys on. + dt_read <- if (is_count) "P1D" else w$dt cache_key <- stac_composite_cache_key( - aoi_target, res, target_crs, w$datetime, w$dt, aggregation, resampling, + aoi_target, res, target_crs, w$datetime, dt_read, aggregation, resampling, cfg$stac_url, cfg$collection, band_assets, bands, cloud_cover_max, - mask_values, scale, offset, offset_before, months, clip, tile_size + mask_values, scale, offset, offset_before, months, clip, tile_size, + family = family ) cache_file <- file.path(cache_source_dir, - paste0("composite_", cache_key, ".tif")) + paste0(family, "_", cache_key, ".tif")) if (!force && file.exists(cache_file) && - cache_hit_ok(cache_file, "composite")) { + cache_hit_ok(cache_file, family)) { r <- stac_cube_cache_read(cache_file, cfg$collection, w$datetime, - label = "composite") + label = family) if (!is.null(r)) { - message(" composite ", label, ": cached") + message(" ", family, " ", label, ": cached") return(composite_finish(r, bands, w$t0)) } } @@ -207,15 +219,29 @@ dft_stac_composite <- function(aoi, build <- function() { fetched <- stac_cube_items(cfg, aoi_wgs84, w$query, cloud_cover_max, months, sign_fn) - composite_offset_check(fetched$is_pre, label, cfg$offset_boundary) + if (is_count) { + # No offset split. stac_cube_assemble() coalesces the two sides with + # terra::cover(pre, post), which keeps the pre-side count wherever it is + # not NA and so would drop every post-side day. The offset does not + # change which pixels are masked, so one pass counts the window whole. + fetched$is_pre[] <- FALSE + } else { + composite_offset_check(fetched$is_pre, label, cfg$offset_boundary) + } stk <- stac_cube_assemble( fetched, cfg, aoi_target, target_crs, t0 = w$t0, t1 = w$t1, res = res, - dt = w$dt, aggregation = aggregation, resampling = resampling, + dt = dt_read, + # "count" is not a cube_view aggregation: gdalcubes would read it as + # "none" and return reflectance (#92). Within a day, "first" skips + # masked (NaN) items, so a clear tile is not blanked by a cloudy one. + aggregation = if (is_count) "first" else aggregation, + resampling = resampling, band_assets = band_assets, mask_values = mask_values, offset = offset, offset_before = offset_before, pixel_fn = pixel_fn, tile_size = tile_size ) stk <- composite_layers_order(stk, bands, label, w) + if (is_count) stk <- count_zero_na(stk) if (isTRUE(clip)) stk <- stac_cube_clip(stk, aoi_target) cube_check_nonempty(stk, cfg$collection, w$datetime, cached = FALSE) # Names, and NO time, before the write. The stack arrives carrying a time @@ -225,11 +251,15 @@ dft_stac_composite <- function(aoi, # below, as it is on a cache hit. names(stk) <- bands terra::time(stk) <- NULL + # A count is a whole number of days, stored as one; its overviews take + # the nearest cell rather than an average, which would not be a count. cache_write_atomic(cache_file, function(path) { terra::writeRaster( - stk, path, filetype = "COG", datatype = "FLT4S", overwrite = TRUE, + stk, path, filetype = "COG", + datatype = if (is_count) "INT2U" else "FLT4S", overwrite = TRUE, gdal = c("COMPRESS=DEFLATE", "PREDICTOR=YES", "BLOCKSIZE=512", - "OVERVIEW_RESAMPLING=AVERAGE") + paste0("OVERVIEW_RESAMPLING=", + if (is_count) "NEAREST" else "AVERAGE")) ) }) # Return what was cached, not the in-memory double stack: the file is @@ -258,6 +288,38 @@ dft_stac_composite <- function(aoi, } +#' Count the clear days per pixel in a daily masked cube, one layer per band +#' +#' `cube` is the masked `raster_cube()` over a window at `dt = "P1D"`, so each +#' time step is one day and same-day items (overlapping MGRS tiles) have already +#' been reduced to one value. `count()` counts the non-NaN, i.e. unmasked, days. +#' The built-in string reducer runs in C++, so the reduce_time() closure trap in +#' inst/notes/gdalcubes-pc-gotchas.md does not apply. `names = bands` gives the +#' layers their role names, which composite_layers_order() requires. The window +#' always spans more than one day, which matters: reduce_time() passes a +#' single-step cube through unchanged. +#' @noRd +composite_count_cube <- function(cube, band_assets, bands) { + gdalcubes::reduce_time(cube, paste0("count(", band_assets, ")"), + names = bands) +} + + +#' A pixel with no clear day is NA, never 0 +#' +#' gdalcubes returns such a pixel as 0 when its chunk holds a clear pixel +#' somewhere and as NaN when the whole chunk is empty, and the chunk size follows +#' `parallel`. Measured on a 128 x 128 fixture: 12,544 zeros at 256 px chunks, +#' 4,352 zeros and 8,192 NaN at 64 px. Mapping 0 to NA makes the output the same +#' at every chunking, so `parallel` stays a cost-only setting. It is NA rather +#' than 0 because a chunk that failed to read is NaN as well, and a failed read +#' must not be published as a count of zero clear days. +#' @noRd +count_zero_na <- function(stk) { + terra::classify(stk, cbind(0, NA)) +} + + #' Warn that a year was dropped, and return NULL for it #' @noRd composite_skip <- function(label, why) { @@ -419,20 +481,26 @@ composite_label <- function(year, months) { #' Cache key for one composite #' #' Its own function rather than a reuse of stac_cube_cache_key(), whose legacy -#' shape is frozen. The leading `"composite"` tag keeps the two families apart -#' even if their parameter lists ever coincide; the filename prefix already -#' does, and this makes it true of the hash as well. Bands and their assets are +#' shape is frozen. The leading tag keeps the families apart even if their +#' parameter lists ever coincide; the filename prefix already does, and this +#' makes it true of the hash as well. Bands and their assets are #' order-sensitive (they are the layer order); months and mask values are not. +#' +#' `family = "count"` is the clear-observation count (#92). It keys apart from +#' every composite, and in particular from the 0.19.x `composite_.tif` +#' files written under `aggregation = "count"`, which hold reflectance: the +#' tag, and the `P1D` read step the caller passes, both change the hash. The +#' default keeps every existing composite key unchanged. #' @noRd stac_composite_cache_key <- function(aoi_target, res, target_crs, datetime, dt, aggregation, resampling, stac_url, collection, band_assets, bands, cloud_cover_max, mask_values, scale, offset, offset_before, months, clip, - tile_size = NULL) { + tile_size = NULL, family = "composite") { geom_wkb <- sf::st_as_binary(sf::st_geometry(aoi_target), endian = "little") parts <- list( - "composite", geom_wkb, as.numeric(res), target_crs, datetime, dt, + family, geom_wkb, as.numeric(res), target_crs, datetime, dt, aggregation, resampling, stac_url, collection, band_assets, bands, as.numeric(cloud_cover_max), sort(as.numeric(mask_values)), as.numeric(scale), as.numeric(offset), as.numeric(offset_before), diff --git a/R/dft_stac_cube.R b/R/dft_stac_cube.R index f952839..a122c5c 100644 --- a/R/dft_stac_cube.R +++ b/R/dft_stac_cube.R @@ -567,6 +567,9 @@ stac_cube_assemble <- function(fetched, cfg, aoi_target, target_crs, t0, t1, res, dt, aggregation, resampling, band_assets, mask_values, offset, offset_before, pixel_fn, tile_size = NULL) { + # the last point before cube_view(): an aggregation it does not honour is + # read as "none" and returns plausible pixels (#92) + aggregation_check(aggregation) mask_asset <- cfg$roles$mask features <- fetched$features is_pre <- fetched$is_pre diff --git a/inst/notes/gdalcubes-pc-gotchas.md b/inst/notes/gdalcubes-pc-gotchas.md index 9f686c7..c4e0422 100644 --- a/inst/notes/gdalcubes-pc-gotchas.md +++ b/inst/notes/gdalcubes-pc-gotchas.md @@ -211,3 +211,16 @@ Measured on gdalcubes 0.7.5, rstac 1.0.1 and terra 1.9.50, with the BULK floodpl - Writing that to a COG makes terra emit a `.aux.json` sidecar. So do units, varnames, longnames, metags and scoff; names alone do not. - **A date-only STAC end bound is read as 00:00Z**, so it drops that day's scenes, which in BC land around 19:00Z. A window ending on a scene day returned 22 of 23 items; the same window with `T23:59:59Z` returned 23. The composite uses explicit times. The cube still passes the user's string (#83). - **CMR-STAC (NASA LPCLOUD, for HLS) differs from PC.** It returns 500 on `intersects` (use `bbox`) and ignores the CQL2 `eo:cloud_cover` filter (filter client-side). Its assets need a real Earthdata Login; `earthdatalogin`'s defaults return 401. See #82. + +## Silent fallbacks and counting clear observations (#92, 2026-09-30) + +- **`cube_view()` does not refuse an aggregation it does not know. It reads it as `"none"`.** The R wrapper checks only that the value is one string. The C++ side lower-cases it, and maps anything unrecognised to `AGG_NONE`, which copies every image in. That copy includes NaNs, so a masked item can blank a clear one. This is how `aggregation = "count"` came back as plausible reflectance. + - Measured by round-tripping values through `cube_view()$aggregation`. `min`, `max`, `mean`, `median`, `first` and `last` survive, and so do the undocumented `count_values` and `count_images`. `count`, `sum` and `""` become `none`. + - drift passes only the first six (`aggregation_check()`). `count_*` count **items**, so overlapping MGRS tiles double-count one acquisition, and every drift caller would scale the result as reflectance. + - `resampling` has the same fallback, to `near` (`bilinaer -> near`), in #96. +- **Count clear days, not items: `dt = "P1D"`, then `reduce_time(cube, "count(B04)", names = ...)`.** Masked pixels are NaN before aggregation, and every aggregation tried (first, max, median) skips NaN within a day. So two same-day tiles give one clear day, and a cloudy tile does not blank a clear one. The string reducer is C++, so the R-callback closure trap above does not apply. + - `reduce_time()` passes a single-time-step cube through unchanged. A count over a one-day window is therefore the input, not a count. +- **A zero-clear pixel is 0 or NaN depending on chunk layout.** An all-NaN chunk stays empty (NaN), while a pixel in a chunk that holds a clear value somewhere reads 0. Measured on a 128 x 128 fixture: 12,544 zeros at 256 px chunks, against 4,352 zeros plus 8,192 NaN at 64 px. Chunk size follows `gdalcubes_options(parallel =)`, so the raw output depends on a setting drift documents as cost-only. + - `count_zero_na()` maps 0 to NA, after which the chunkings agree. + - NA rather than 0, because a failed chunk read is NaN too. +- **`create_image_collection(one_band_per_file = FALSE)` on two-band GeoTIFFs segfaulted a worker** (0.7.5). Separate files per band with a format JSON work. The #92 test fixture uses them. diff --git a/planning/active/findings.md b/planning/active/findings.md index e38e5e0..66641f0 100644 --- a/planning/active/findings.md +++ b/planning/active/findings.md @@ -72,8 +72,75 @@ Probe: 6 local items on a 4x4 grid, EPSG:32609, July 2021. The inputs: - Not acted on: the misleading "offset split" message from `stac_cube_items()` under count (cosmetic). - `create_image_collection(one_band_per_file = FALSE)` on 2-band files segfaulted a gdalcubes worker (0.7.5). Separate files per band plus a format JSON work, so the tests use those. +## Live check (2026-09-30, Phase 2 snapshot, `data`: scratchpad `count_live*.R`) + +A 2 km square at the packaged AOI's centroid, July, `bands = "red"`, `res = 20`, EPSG:32609, 10,000 cells: + +| run | items | distinct dates | count min / median / max | NA | wall | peak RSS | +|---|---|---|---|---|---|---| +| 2021, `cloud_cover_max = 20` | 4 | 4 | 4 / 4 / 4 | 0 | 9.6 s | 0.30 GiB (both years) | +| 2023, `cloud_cover_max = 20` | 6 | 4 | 4 / 4 / 4 | 0 | 10.9 s | | +| 2021, `cloud_cover_max = 100` | 12 | 12 | 4 / 6 / 9 | 0 | 19.7 s | 0.36 GiB (both years) | +| 2023, `cloud_cover_max = 100` | 20 | 12 | 5 / 8 / 9 | 0 | 31.9 s | | + +What the runs show: + +- **Overlapping MGRS tiles collapse to one day.** 2023 returned 6 items on 4 dates and counts 4 everywhere, not 6. +- **The mask is applied per pixel.** With cloudy scenes allowed, counts vary spatially and stay at or below the distinct-date count. +- **Every value is whole and every file is INT2U**, with the scale not applied. + +The issue's own call (2021 July, NECR) returned continuous values from 0.0077 to 0.18. + +## Scale: BULK count chips (2026-09-30) + +`data-raw/benchmark_composite_bulk.R chips-count` covers the first 20 of #79's seed-79 chips: 300 m buffers, 2023 Jul-Aug, `bands = "red"`, `aggregation = "count"`, run under `/usr/bin/time -l`. + +- **Wall clock:** 577.5 s real (575.1 s in the loop), per chip median 24.7 s, max 59.6 s. +- **Peak RSS:** 489,013,248 B (0.46 GiB). +- **Output:** 3,660 cells a chip, per-chip max count 5-9, 0 NA cells. +- **Against the median chips:** the same 20 chips as median true-colour composites (#79's `chips_per_chip.csv` rows 1-20) took 1,051.6 s, median 44.4 s. #79 recorded a 0.48 GiB peak over all 100 chips. +- **Why the count is faster:** it streams one band against three, so daily time steps did not cost more at chip scale. + +The floodplain-wide read is out of scope while #88 is open. + +## P2 code-check termination: enumeration of the count build's wiring + +Round 3 named the mechanism: a mock that receives an argument and returns a fixed result hides how the call site wires that argument. Round 1's fix (capturing `pixel_fn`) was an instance of it. Below is every argument the count build passes to its collaborators, with whether a test now observes it. + +**`stac_cube_items()`** + +| argument | observed? | +|---|---| +| `datetime` (`w$query`) | ✓, and a mutation fails | +| `cloud_cover_max` | ✓, and a mutation fails | +| `months` | ✓, and a mutation fails | +| `cfg`, `aoi_wgs84`, `sign_fn` | not observed | + +**`stac_cube_assemble()`** + +| argument | observed? | +|---|---| +| `fetched$is_pre`, `t0`, `t1`, `dt`, `aggregation`, `band_assets`, `pixel_fn`, `resampling`, `mask_values` | ✓, and each mutation fails | +| `offset`, `offset_before` | ignored by the count's `pixel_fn` by design | +| `cfg`, `aoi_target`, `target_crs`, `res`, `tile_size` | not observed | + +**Key call site** + +| argument | observed? | +|---|---| +| `family`, `dt_read` | not observed at the call site | + +- **The unobserved `stac_cube_items()` and `stac_cube_assemble()` arguments** are shared, unchanged lines of the median path, not introduced by #92. They are accepted as pre-existing. +- **The key call site's `family` and `dt_read`:** dropping either is not a defect. The asserted `count_` prefix already separates the families, so the tag and `P1D` in the hash are defence in depth. +- **The overview resampling was the one count-specific setting left unobserved** (round 3: a fixture smaller than one COG block builds no overviews). It is now tested with a 1100 px checkerboard read at `OVERVIEW_LEVEL=0`: NEAREST gives {1, 5} and AVERAGE gives {3}. + +Six round-3 mutations each turn a test red: overview AVERAGE, `cfg$mask_values`, `months = NULL`, `w$datetime` for `w$query`, constant resampling, constant cloud cover. + ## Errors Encountered | Error | Resolution | |-------|------------| +| P2 round 3: count overviews untested (fixture below one 512 px block); `mask_values`/`months`/query/`resampling` reached the mocks unobserved | 1100 px checkerboard overview test; count test captures and asserts the four with non-default values; each of 6 mutations red. Enumeration above ends the loop | +| P2 round 1: no test drove `dft_stac_composite()`'s count `pixel_fn` (a reflectance closure there stayed green) | The mocked assemble now captures `pixel_fn` and runs it on a fixture cube, expecting 2 clear days named `red`. The mutation turns 2 tests red | +| Round 3: the error message claimed gdalcubes "would not honour" values it does honour (none, count_*) | Reworded to state the fallback and drift's policy. Enumerated all 8 gdalcubes claims in the Phase 1 diff, all measured or test-asserted. Sibling defect on `resampling` filed as #96 | | Round 1: `.cube_view_aggregations` pinned to `?cube_view`, which refuses `"last"` (honoured) and case variants | Probed `cube_view()$aggregation`: honours min/max/mean/median/first/last, count_values, count_images; lower-cases; else `"none"`. Set is now the six, case-insensitive; test pins behaviour via round trip, not the Rd | diff --git a/planning/active/progress.md b/planning/active/progress.md index 93be213..c2844db 100644 --- a/planning/active/progress.md +++ b/planning/active/progress.md @@ -6,3 +6,7 @@ - Created branch `92-dft-stac-composite-aggregation-count-si` off main - Scaffolded PWF baseline from issue #92 with approved phases - Next: start Phase 1 +- Plan review (Plan agent) folded in: `"count"` never reaches cube_view, and the 0/NaN chunk dependence was measured, giving the rule 0 -> NA. Phase 2 was rewritten accordingly. +- Phase 1 committed (727c5d9) after 3 code-check rounds. Round 1 found `last` refused and was fixed by a behaviour round trip. Round 2 was clean. Round 3's message claim was fixed and then enumerated (8 claims). +- Filed #96 (resampling has the same silent fallback) +- Phase 2 implemented. 7 mutations each redden the composite tests. diff --git a/planning/active/review-p1-round3.md b/planning/active/review-p1-round3.md new file mode 100644 index 0000000..4ec12a3 --- /dev/null +++ b/planning/active/review-p1-round3.md @@ -0,0 +1,56 @@ +# Review — #92 Phase 1, round 3 + +Snapshot: `scratchpad/p1snap`, gdalcubes 0.7.5. Brief: find every place in the diff or in the +three functions where a membership list is asserted rather than measured against the consumer. + +## Measured first + +The round-trip test checks the gdalcubes parser (`cube_view()$aggregation`), which is a stand-in +for the reduction that actually runs. So I ran the real thing: three single-band scenes (values +10, 30, 20 on 07-05, 07-15, 07-25), each value put through `raster_cube()` and then `as_array()`. + +| aggregation | parsed as | reduced value | honoured? | +|---|---|---|---| +| min / max / mean / median | same | 10 / 30 / 20 / 20 | yes | +| first / last | same | 10 / 20 | yes (last = the latest date) | +| MEDIAN / Last | median / last | 20 / 20 | yes, ignoring case | +| count, sum | none | 20 | silent fallback | +| none | none | 20 | — | +| count_values / count_images | same | 3 / 3 | yes | + +So the allowed six are correct as behaviour, not only as parser output. The code comment's +claims hold: it lower-cases, `count` falls back to `none` and comes back as a value, and +`count_*` are honoured. Probe script: `scratchpad/aggprobe/probe.R`. Both test files pass +(141 and 61 passing). The new tests do not skip; the only skips are pre-existing +(missing-package path, network). + +## Findings + +- **[fragile]** R/dft_stac_cube.R:322 (`aggregation_check()`, the "x" bullet). The refusal says + `Got "", which gdalcubes would not honour.` That is false for three values the diff + refuses on purpose: `"none"`, `"count_values"` and `"count_images"`. I measured all three as + honoured, and the comment eight lines above says so for the `count_*` pair. It is the same + mechanism as round 1: a claim about the consumer that was not measured. A user who passes + `"count_values"` is told gdalcubes cannot do it, when drift is refusing it by policy (items + are double-counted over MGRS overlap). The message will also be wrong in Phase 2 for any value + the composite's extended `allowed` refuses. Suggested fix: say it is not supported by drift, + not that gdalcubes would not honour it. The two test titles that make the same claim + (`test-dft_stac_cube.R:758`, `test-dft_stac_composite.R:299`) are only labels, not defects. + +- **[bug, pre-existing, not introduced by this diff]** R/dft_stac_cube.R:601, + R/dft_stac_fetch.R:616, and the composite path through `stac_cube_assemble()`. `resampling` + reaches `gdalcubes::cube_view()` in the same three functions with no validation, and it + falls back the same silent way, to `"near"`: + `"bilinaer" -> near`, `"foo" -> near`, `"Bilinear" -> bilinear`, `"nearest" -> near`. + So a typo in `resampling` on the continuous or composite path returns a + nearest-neighbour cube, cached under the typo's key, with no error. This is the #92 defect + on the sibling argument. Out of this phase's scope; file it rather than widen the PR. + +No other asserted list is reached: +- The roxygen `@param aggregation` text in all three functions names only the defaults. It + lists no set and makes no claim about gdalcubes. +- The test's refused list is drift policy, and it agrees with the measurement above. +- The allowed list in the test (`c("min", ... "last")`) duplicates `.cube_view_aggregations`, + but the round-trip test iterates the constant itself, so the constant cannot drift away from + gdalcubes. +- No data-raw script, vignette or test passes an aggregation outside the six. diff --git a/planning/active/review-p2-round1.md b/planning/active/review-p2-round1.md new file mode 100644 index 0000000..2a3fb90 --- /dev/null +++ b/planning/active/review-p2-round1.md @@ -0,0 +1,22 @@ +# Code-check review: #92 Phase 2, round 1 + +Snapshot: `scratchpad/p2snap`, gdalcubes 0.7.5. The composite test file passes (89 pass, 2 network skips). + +## Findings + +- **[fragile]** tests/testthat/test-dft_stac_composite.R, test "dft_stac_composite(aggregation = 'count') reads daily, unsplit, and caches integers" (the `stac_cube_assemble` mock). Nothing tests that `dft_stac_composite()` connects the count path to `composite_count_cube()`. The mock takes `pixel_fn` through `...` and never calls it. The two fixture tests call `drift:::composite_count_cube()` directly and do not go through `dft_stac_composite()`. **Mutation run on a copy (`scratchpad/p2mut`):** I changed the count branch of `pixel_fn` (R/dft_stac_composite.R:186) to the reflectance `apply_pixel(scale_token(...))` closure. The suite stayed green, `[ FAIL 0 | SKIP 2 | PASS 89 ]`. So if `pixel_fn` were inverted or dropped, the count path would return scaled reflectance read at P1D with no `reduce_time`, which is several layers. The tests would stay green for that part. The `composite_layers_order()` layer-count abort would still fire at run time, but only against a live STAC. This is the one wiring defect of the #92 class that the 7 accepted mutations do not cover. A fix that can fail: in the mock, capture `pixel_fn`, and assert that `seen$pixel_fn` applied to a `count_fixture()` cube returns the day counts. Or, more cheaply, assert that it returns a `reduce_time` cube. For example, `inherits(seen$pixel_fn(cube, 0), "reduce_time_cube")`. + +## Probed and clean + +- **NA round trip through INT2U COG:** a 600x600x3 stack of {0, NaN, NA, 1, 7, 40} went through `count_zero_na()`, then `writeRaster(COG, INT2U, NEAREST)`. The NoData value is 65535, and every 0, NaN and NA read back as NA. There was no 0 and no 65535, and the NA mask matched the pre-write mask exactly. `cache_invalid_reason()` gives NA, and `stac_cube_cache_read()` returns the raster and not NULL. `gdalinfo` shows UInt16 and NEAREST overviews. +- **Count across time chunks:** I made the clear days fall only in the later of the 16-day time chunks, with the whole grid cloudy in the first. At chunkings `c(16,256,256)`, `c(16,64,64)`, `c(1,64,64)`, `c(4,128,128)` and `c(31,64,64)`, the counts were identical and correct (1 for 12,544 pixels, 2 for 3,840). +- **Order within a day:** I swapped which same-day item is cloudy, so the cloudy one sorts first. `first`, `last`, `max` and `median` all still count the day. The test fixture's ordering is therefore not load-bearing. +- **Window edges:** items on 07-01 and 07-31 are both counted with `t1 = "2021-07-31"` and `dt = "P1D"`. +- **Multi-band naming:** `composite_count_cube(cube, c("B04","SCL"), c("red","blue"))`, written to NetCDF, is read by terra as `blue, red` (alphabetical). `composite_layers_order()` restores `red, blue` by name. +- **Can "count" reach cube_view:** no. `dft_stac_cube()` and `dft_stac_fetch()` keep the default allowed set. The composite passes `"first"`, and `stac_cube_assemble()` rejects `"count"` with the default set. +- **clip / tile_size ordering:** `count_zero_na()` runs after `mosaic_stacks()`, which receives the whole assembled stack, and before `stac_cube_clip()`. Tiles do not overlap, and both orderings give the same NA set. `tile_size` is in the count key. +- **Cache family:** no other code matches on the `composite_` or `cube_` prefix, so `count_` breaks no cache consumer. The default family keeps the literal `"composite"` first element. + +## Note (not a finding) + +`count_zero_na()` adds one full-grid `terra::classify()` intermediate, in memory, next to the existing in-memory `stac_cube_clip()`. The composite is sized for AOIs and chips, so this is only relevant if a count is run at BULK scale. diff --git a/planning/active/review-p2-round2.md b/planning/active/review-p2-round2.md new file mode 100644 index 0000000..9a809aa --- /dev/null +++ b/planning/active/review-p2-round2.md @@ -0,0 +1,63 @@ +# Code-check review: #92 Phase 2, round 2 + +Snapshot: `scratchpad/p2snap`, gdalcubes 0.7.5. Work copy: `scratchpad/p2r2/mut`, restored and +checked with `cmp` afterwards. On the unmutated copy the composite test file passes: +`[ FAIL 0 | SKIP 2 | PASS 92 ]`. + +## Clean + +No bug, security or fragility finding. The round-1 fix discriminates, and the rest of the diff +held up under probing. + +## The fix: does the pixel_fn assertion discriminate? + +I swapped the count branch of `pixel_fn` (R/dft_stac_composite.R:186) for each variant below and +ran the test file: + +| pixel_fn mutation | result | +|---|---| +| identity (`cube`, a 31-step passthrough) | red (2 fail) | +| `reduce_time(sum(B04))`, which is reflectance-scale | red | +| `count(...)` with `names = band_assets` | red | +| scale-then-count (`apply_pixel(B04*1e-4-0.1)` then `count`) | green, but equivalent: scaling does not change which cells are NaN, so the count is the same | +| `count(SCL)` | green, but equivalent: `image_mask` NaN-masks every band, the mask band included, so it gives the same count | + +Both survivors produce the same output as the real code, so they are not defects the test ought +to catch. Round 1 already accepted the reflectance `apply_pixel` closure as a red mutation. + +- **Counting items rather than days** cannot happen inside `pixel_fn`. Same-day items are + collapsed by the `cube_view` (P1D, `"first"`) before `pixel_fn` sees the cube. So that property + belongs to the view, and it is pinned by `seen$dt == "P1D"` and `seen$aggregation == "first"`, + plus the separate fixture test with duplicate same-day items. The pixel_fn fixture having one + item per day does not weaken it. +- **Single-step passthrough:** the fixture view runs 07-01..07-31 at P1D, which is 31 steps, so + it cannot reach it. I measured the passthrough claim in the docstring and the notes: a one-day + view (`t0 = t1`) returns `500` (the input), and a two-day view returns `1`. A composite window + is always at least one whole month, so production never reaches it. + +## Probed, not covered by round 1 + +- **Swath-edge nodata as a "clear day".** SCL 0 (NO_DATA) is not in `mask_values`, and + `stac_image_collection()` sets band `nodata = ""`. So the only guard is the COG's own + NoData. Fixture probe: when the B04 file declares `NoData=0`, the 0 pixels are not counted + (1 vs 2). When the file declares no NoData, they are counted, which would inflate the count. + Checked live: a PC `sentinel-2-l2a` item (S2A_MSIL2A_20210710T192911_R142_T09UXA) has + `NoData Value=0` on both B04 and SCL. So production is safe. The test fixture declares + `"nodata": 0` in its format JSON instead of in the file, and the effect is the same. +- **Case:** `aggregation_check()` returns lower case, so `"COUNT"` → `is_count` TRUE. +- **Claims in the notes:** #96 exists and is the resampling fallback. +- **Cache-hit path:** a count file is not passed through `count_zero_na()` again, and does not + need to be (NA is stored as NoData 65535). + +## Notes (not findings) + +- The docstring of `count_zero_na()` justifies NA-not-0 by saying that "a failed chunk read is + NaN". That is true when every time chunk of a spatial block fails. A **partial** failure + (#87: an expired token or an HTTP error on some days' chunks, reported only on stderr) is + reduced as missing days and published as a plausible undercount (for example 5 instead of 7). + The median composite has the same exposure (a plausible wrong median), so this is the + existing #87 class and not new with this diff. The Phase 3 docs could state that the count is + a lower bound under read failure. +- `stac_cube_items()` still prints "offset split at …: N pre / M post" on a straddling count + window before `is_pre` is reset, so the message describes a split that does not happen. This + is cosmetic. diff --git a/planning/active/review-p2-round3.md b/planning/active/review-p2-round3.md new file mode 100644 index 0000000..4820e56 --- /dev/null +++ b/planning/active/review-p2-round3.md @@ -0,0 +1,50 @@ +# Code-check review: #92 Phase 2, round 3 + +Snapshot: `scratchpad/p2snap`, copied to `scratchpad/p2r3/base` and mutated only in `scratchpad/p2r3/mut_*`. Baseline for `test-dft_stac_composite.R`: `[ FAIL 0 | SKIP 2 | PASS 92 ]`. No other test file calls `dft_stac_composite()` or `stac_composite_cache_key()`, so this file is the whole guard. Harness: `scratchpad/p2r3/run_mut.sh`. + +## Mechanism + +Round 1's defect was not specific to `...`. **A mock that receives an argument and returns a fixed result makes the call site's wiring of that argument unobservable.** Passing the argument through `...` and naming it as a formal the mock never reads are the same hole. The only wiring a test asserts is what the mock copies into `seen`, so the assertions cover exactly the arguments the author already suspected. Round 1 was that selection missing `pixel_fn`. + +The same mechanism shows up at a second seam. `stac_composite_cache_key()` is not mocked in the count tests, but its call-site output is checked only by the filename regex `^count_[0-9a-f]{16}\.tif$`. So every argument the call site feeds it is unobserved too, the same way. + +A third, adjacent shape is a fixture that cannot reach the failure mode. The mocked stack is 63 x 65 px, and a COG with `BLOCKSIZE=512` builds no overviews at that size. So the overview-resampling argument is written but never exercised. + +## Where it reaches in this diff's tests + +| seam | args the mock ignores | wiring asserted elsewhere? | +|---|---|---| +| `stac_cube_items` (`function(...)` in every mock: l.169, 272, 302, 429, 490, 504) | cfg, aoi_wgs84, **w$query**, cloud_cover_max, **months**, sign_fn | no. The only check is `composite_window()$query` in isolation, and the network tests are skipped by default | +| `stac_cube_assemble`, count mock (l.432) | named but not read: cfg, aoi_target, target_crs, res, **resampling**, **mask_values**, offset_before. Through `...`: **tile_size** | no. The composite build mock (l.273) and the empty-count mock (l.505) are `function(...)` and read nothing | +| `stac_cube_assemble`, count mock | captured: is_pre, dt, aggregation, t0, t1, band_assets, pixel_fn, offset | yes. Every count-specific value is asserted (offset only feeds pixel_fn, which ignores it for a count) | +| `stac_composite_cache_key` call site (l.200-205) | all args, including **dt_read** and **family** | only the `count_` prefix and 16-hex shape. The `ckey()` tests call the function directly, not through the call site | +| COG write (l.258-263) | `OVERVIEW_RESAMPLING` | no. The fixture is too small for overviews | + +## Mutations run (each on a fresh copy, each green = `[ FAIL 0 | SKIP 2 | PASS 92 ]`) + +| mutation (R/dft_stac_composite.R) | result | real defect? | +|---|---|---| +| l.262 count overviews `NEAREST` -> `AVERAGE` | green | yes (Finding 1) | +| l.220 `w$query` -> `w$datetime` | green | yes, pre-existing (Finding 2) | +| l.221 `months` -> `NULL` | green | yes, pre-existing (Finding 2) | +| l.239 `mask_values = mask_values` -> `cfg$mask_values` | green | yes, pre-existing (Finding 2) | +| l.238 `resampling = resampling` -> `"near"` | green | yes, pre-existing (Finding 2) | +| l.241 drop `tile_size = tile_size` | green | cost only: an untiled read with identical values | +| l.204 drop `family = family` from the key call | green | no, see Note | +| l.201 `dt_read` -> `w$dt` in the key call | green | no, see Note | + +## Findings + +- **[fragile]** R/dft_stac_composite.R:262. The count's `OVERVIEW_RESAMPLING=NEAREST` is new in this diff and nothing asserts it. Reverting it to `AVERAGE` stays green, because the mocked 63 x 65 stack is below the 512 px block size and GDAL builds no overviews. This is a fixture that cannot reach the failure mode, not a mock swallowing an argument. Any AOI wider than 512 px at `res` would get averaged INT2U overviews: titiler and QGIS show a rounded mean at low zoom, which is not a count, and the code comment on l.254-255 says exactly that. For the guard to fire, the mocked stack has to exceed 512 px on one side, for example 1024 x 16 at `resolution = 10` over a matching extent. Then read an overview back (`gdalinfo` "Overviews:" plus a value check, or `terra::rast(f, lyrs=…)` of the overview). A cheaper option is to assert the `OVERVIEW_RESAMPLING` metadata that `gdalinfo` reports in the COG's IMAGE_STRUCTURE. + +- **[fragile, pre-existing, not introduced by this diff]** R/dft_stac_composite.R:220-221 and 238-239. Four user-facing arguments reach the two stubbed seams and are never observed there. Each wiring defect below stays green, and each changes pixels, not just cost: + - `w$query` -> `w$datetime` brings back #83 (a date-only end bound drops the last day's scenes). A count silently loses one clear day. + - `months` -> `NULL` removes the month filter (`stac_cube_items()` filters only when `!is.null(months)`). `months = c(6, 8)` then counts or medians July scenes as well. + - `mask_values` -> `cfg$mask_values` ignores the user's mask. For a count this argument is the definition of "clear". The cache key still hashes the user's value, so the file is labelled with a mask it was not built with. + - `resampling` -> a constant has the same shape: it is keyed on the user's value and read with another. + + The count mock already names `mask_values` and `resampling` as formals and a `...` that holds `tile_size`. Copying them into `seen`, and giving `stac_cube_items` named formals that capture `datetime` and `months`, would make all four red. Before this diff the composite build test mock was `function(...)`, so this gap is older than #92. It is listed because the new count test repeats the same hole on the path whose meaning depends most on `mask_values`. + +## Note (not a finding): correction to the accepted list + +"Key without family" is red only when the mutation is inside `stac_composite_cache_key()`, where `ckey(family = "count") != ckey()` catches it. At the call site (l.204), dropping `family = family` stays green, and so does keying on `w$dt` instead of `dt_read` (l.201). Neither is a real defect. The asserted `count_` file prefix already separates the count files from every `composite_*.tif`, including the 0.19.x reflectance files written under `aggregation = "count"`. That separation holds even with both mutations applied together, which reproduces the stale hash `08e0c5510e8ae297` but under a `count_` name. The tag and `P1D` in the hash are defence in depth, not the separator. The docstring on l.490-494 ("the tag, and the `P1D` read step the caller passes, both change the hash") describes the function correctly, but nothing pins the caller to it. diff --git a/planning/active/task_plan.md b/planning/active/task_plan.md index 17749a5..84f6f82 100644 --- a/planning/active/task_plan.md +++ b/planning/active/task_plan.md @@ -35,42 +35,42 @@ Decisions taken at the gate: Revised after the plan review (see findings.md, "Plan review" and "Count semantics"). -- [ ] The count path calls `stac_cube_assemble()` with `dt = "P1D"` and a day +- [x] The count path calls `stac_cube_assemble()` with `dt = "P1D"` and a day aggregation of `"first"`, never `"count"`. `"count"` must never reach `cube_view`, because it maps to `AGG_NONE` (the #92 bug). As a last guard, `aggregation_check()` also runs inside `stac_cube_assemble()`. -- [ ] Add a helper, `composite_count_cube(cube, band_assets, bands)`, that returns +- [x] Add a helper, `composite_count_cube(cube, band_assets, bands)`, that returns `reduce_time(cube, paste0("count(", band_assets, ")"), names = bands)`. It applies no scale or offset, and `names = bands` satisfies `composite_layers_order()`. -- [ ] **A pixel with zero clear days is NA, always** (`count_zero_na()`, `0 -> NA`). +- [x] **A pixel with zero clear days is NA, always** (`count_zero_na()`, `0 -> NA`). gdalcubes returns 0 or NaN for it depending on chunk layout, and chunk size follows `parallel`. This rule is therefore what keeps `parallel` a cost-only knob. It fails toward NA because a failed chunk read cannot be told apart from "all cloudy". An all-NA result is caught by `cube_check_nonempty()` and the year is dropped (`drift_empty_cube`). -- [ ] **No offset split for count.** Set `is_pre` to all FALSE and skip +- [x] **No offset split for count.** Set `is_pre` to all FALSE and skip `composite_offset_check()`, because `terra::cover(pre, post)` would drop the post side. A window that straddles 2022-01-25 is valid for a count. -- [ ] Give the count its own cache family: the `"count"` tag in +- [x] Give the count its own cache family: the `"count"` tag in `stac_composite_cache_key()` (a family argument whose default is `"composite"`, so the median key is byte-identical, pinned to `03ee8ecc66b832a8`), the key hashing `dt = "P1D"`, and the file prefix `count_`. The count key must not equal the stale 0.19 key `08e0c5510e8ae297`. Write the file as INT2U with `OVERVIEW_RESAMPLING=NEAREST`, and use `"count"` as the label in the cache-read and hit messages. -- [ ] Offline tests: - - [ ] Test `composite_count_cube()` on a local `create_image_collection()` +- [x] Offline tests: + - [x] Test `composite_count_cube()` on a local `create_image_collection()` fixture of separate B04/SCL files. Cover: a clear same-day item not blanked by a masked one, a count of distinct days rather than items, integer output with no scale, and names equal to the roles. - - [ ] Test that the result is chunking-invariant after `count_zero_na()` + - [x] Test that the result is chunking-invariant after `count_zero_na()` (`chunking = c(16, 64, 64)` against `c(16, 256, 256)`). - - [ ] Build `dft_stac_composite(aggregation = "count")` offline with a mocked + - [x] Build `dft_stac_composite(aggregation = "count")` offline with a mocked `stac_cube_assemble` that captures its arguments. Assert `"first"` and `"P1D"`, that `is_pre` is all FALSE with no refusal for a straddling window, the `count_.tif` name, INT2U, and that zeros become NA. - - [ ] Restore each defect and confirm the matching test goes red. -- [ ] Add the `AGG_NONE` root cause and the chunk-dependent 0/NaN behaviour to + - [x] Restore each defect and confirm the matching test goes red. +- [x] Add the `AGG_NONE` root cause and the chunk-dependent 0/NaN behaviour to `inst/notes/gdalcubes-pc-gotchas.md`. ## Phase 3: Docs diff --git a/tests/testthat/test-dft_stac_composite.R b/tests/testthat/test-dft_stac_composite.R index 157541a..16c71ad 100644 --- a/tests/testthat/test-dft_stac_composite.R +++ b/tests/testthat/test-dft_stac_composite.R @@ -309,3 +309,270 @@ test_that("dft_stac_composite refuses an aggregation gdalcubes would not honour } expect_length(list.files(cache, recursive = TRUE), 0L) }) + +# A local, dated Sentinel-2-shaped collection: one B04 and one SCL GeoTIFF per +# item, on an n x n grid at 10 m in EPSG:32609. `items` is a list of +# list(id, date, scl), `scl` a function of cell x/y returning SCL classes. Separate +# files per band: two-band files with one_band_per_file = FALSE segfaulted a +# gdalcubes 0.7.5 worker (#92). +count_fixture <- function(items, n = 4) { + d <- withr::local_tempdir(.local_envir = parent.frame()) + for (it in items) { + r <- terra::rast(nrows = n, ncols = n, xmin = 0, xmax = n * 10, ymin = 0, + ymax = n * 10, crs = "EPSG:32609") + xy <- terra::xyFromCell(r, seq_len(terra::ncell(r))) + b <- r + terra::values(b) <- 500 + s <- r + terra::values(s) <- it$scl(xy[, 1], xy[, 2]) + stem <- file.path(d, paste0(format(as.Date(it$date), "%Y%m%d"), "_", it$id)) + terra::writeRaster(b, paste0(stem, "_B04.tif"), datatype = "INT2U") + terra::writeRaster(s, paste0(stem, "_SCL.tif"), datatype = "INT1U") + } + # gdalcubes collection format, written literally (jsonlite is not a dependency) + fmt_file <- file.path(d, "format.json") + writeLines(c( + '{"description": "drift #92 fixture", "tags": ["test"],', + ' "pattern": ".*\\\\.tif",', + ' "images": {"pattern": ".*/([0-9]{8}_[a-z]+)_.*"},', + ' "datetime": {"pattern": ".*/([0-9]{8})_.*", "format": "%Y%m%d"},', + ' "bands": {"B04": {"pattern": ".*_B04\\\\.tif", "nodata": 0},', + ' "SCL": {"pattern": ".*_SCL\\\\.tif"}}}' + ), fmt_file) + col <- gdalcubes::create_image_collection( + list.files(d, pattern = "\\.tif$", full.names = TRUE), format = fmt_file, + out_file = file.path(d, "col.db") + ) + list(col = col, n = n) +} + +# Run the count path's reduction over a fixture: the same view settings and the +# same helpers dft_stac_composite() uses for aggregation = "count". +count_run <- function(fx, chunking = NULL, t0 = "2021-07-01", + t1 = "2021-07-31", aggregation = "first") { + v <- gdalcubes::cube_view( + srs = "EPSG:32609", + extent = list(left = 0, right = fx$n * 10, bottom = 0, top = fx$n * 10, + t0 = t0, t1 = t1), + dx = 10, dy = 10, dt = "P1D", aggregation = aggregation, + resampling = "near" + ) + m <- gdalcubes::image_mask("SCL", values = c(3, 8, 9, 10, 11)) + cube <- if (is.null(chunking)) { + gdalcubes::raster_cube(fx$col, v, mask = m) + } else { + gdalcubes::raster_cube(fx$col, v, mask = m, chunking = chunking) + } + tmp <- withr::local_tempfile(fileext = ".nc", .local_envir = parent.frame()) + gdalcubes::write_ncdf(drift:::composite_count_cube(cube, "B04", "red"), tmp) + drift:::count_zero_na(terra::rast(tmp)) +} + +test_that("the count path counts clear DAYS, not items or reflectance (#92)", { + skip_if_not_installed("gdalcubes") + clear <- function(x, y) rep(4, length(x)) + cloud <- function(x, y) rep(8, length(x)) + left_clear <- function(x, y) ifelse(x < 20, 4, 9) + fx <- count_fixture(list( + # same day, two tiles, both clear: one acquisition, counted once + list(id = "a", date = "2021-07-03", scl = clear), + list(id = "b", date = "2021-07-03", scl = clear), + list(id = "a", date = "2021-07-08", scl = cloud), + list(id = "a", date = "2021-07-13", scl = left_clear), + # same day: tile a cloudy, tile b clear. Under the #92 "none" aggregation + # the masked item overwrote the clear one; the count must keep it. + list(id = "a", date = "2021-07-18", scl = cloud), + list(id = "b", date = "2021-07-18", scl = clear) + )) + r <- count_run(fx) + expect_equal(names(r), "red") + v <- terra::values(r)[, 1] + x <- terra::xyFromCell(r, seq_len(terra::ncell(r)))[, 1] + # left half clear on 07-03, 07-13, 07-18; right half on 07-03, 07-18 + expect_equal(unique(v[x < 20]), 3) + expect_equal(unique(v[x > 20]), 2) + # whole numbers of days, not reflectance and not a scaled count + expect_true(all(v == round(v))) + # the day-level aggregation does not change the count: every one skips NaN + for (agg in c("max", "median")) { + expect_equal(terra::values(count_run(fx, aggregation = agg))[, 1], v) + } +}) + +test_that("a pixel with no clear day is NA at every chunking (#92)", { + skip_if_not_installed("gdalcubes") + # 128 x 128: clear only in a left strip, and in a corner on a second day. + # Raw gdalcubes gives a zero-clear pixel 0 or NaN depending on whether its + # chunk holds a clear pixel; count_zero_na() must make the two agree. + fx <- count_fixture(list( + list(id = "a", date = "2021-07-05", scl = function(x, y) ifelse(x < 300, 4, 8)), + list(id = "a", date = "2021-07-10", + scl = function(x, y) ifelse(x < 300 & y < 300, 4, 8)), + list(id = "a", date = "2021-07-15", scl = function(x, y) rep(8, length(x))) + ), n = 128) + coarse <- count_run(fx, chunking = c(16, 256, 256)) + fine <- count_run(fx, chunking = c(16, 64, 64)) + expect_equal(terra::values(fine), terra::values(coarse)) + v <- terra::values(coarse)[, 1] + expect_false(any(v == 0, na.rm = TRUE)) + expect_equal(sum(is.na(v)), 128^2 - 30 * 128) + expect_equal(sort(unique(stats::na.omit(v))), c(1, 2)) +}) + +test_that("dft_stac_composite(aggregation = 'count') reads daily, unsplit, and caches integers (#92)", { + skip_if_not_installed("gdalcubes") + aoi <- aoi_pkg() + aoi_t <- sf::st_transform(aoi, 32609) + seen <- NULL + queried <- NULL + testthat::local_mocked_bindings( + # a straddling window: the composite would refuse it; a count must not + stac_cube_items = function(cfg, aoi_wgs84, datetime, cloud_cover_max, + months, sign_fn) { + queried <<- list(datetime = datetime, cloud_cover_max = cloud_cover_max, + months = months) + list(features = list(), is_pre = c(TRUE, FALSE, FALSE)) + }, + stac_cube_assemble = function(fetched, cfg, aoi_target, target_crs, t0, t1, + res, dt, aggregation, resampling, band_assets, + mask_values, offset, offset_before, pixel_fn, + ...) { + seen <<- list(is_pre = fetched$is_pre, dt = dt, aggregation = aggregation, + t0 = t0, t1 = t1, band_assets = band_assets, + pixel_fn = pixel_fn, offset = offset, + resampling = resampling, mask_values = mask_values) + r <- terra::rast(terra::ext(aoi_t), resolution = 50, crs = "EPSG:32609") + terra::values(r) <- rep_len(c(0, 1, 2, 5), terra::ncell(r)) + names(r) <- "red" + terra::time(r) <- as.Date("1970-01-01") + r + } + ) + cache <- withr::local_tempdir() + # non-default mask, cloud cover and resampling: for a count the mask defines + # "clear", so what the caller passes must be what the read uses + out <- suppressMessages(dft_stac_composite( + aoi, years = 2022, months = 1, bands = "red", aggregation = "count", + mask_values = c(3L, 8L, 9L), cloud_cover_max = 55, resampling = "near", + cache_dir = cache + )) + expect_equal(queried$datetime, "2022-01-01T00:00:00Z/2022-01-31T23:59:59Z") + expect_equal(queried$months, 1L) + expect_equal(queried$cloud_cover_max, 55) + expect_equal(seen$mask_values, c(3L, 8L, 9L)) + expect_equal(seen$resampling, "near") + expect_equal(seen$aggregation, "first") + expect_equal(seen$dt, "P1D") + expect_false(any(seen$is_pre)) + expect_equal(c(seen$t0, seen$t1), c("2022-01-01", "2022-01-31")) + # the pixel function the build hands to gdalcubes is the day count, not the + # reflectance expression: run it on a fixture cube and read the result + expect_equal(seen$band_assets, "B04") + fx <- count_fixture(list( + list(id = "a", date = "2021-07-03", scl = function(x, y) rep(4, length(x))), + list(id = "a", date = "2021-07-09", scl = function(x, y) rep(4, length(x))), + list(id = "a", date = "2021-07-15", scl = function(x, y) rep(9, length(x))) + )) + v <- gdalcubes::cube_view( + srs = "EPSG:32609", + extent = list(left = 0, right = 40, bottom = 0, top = 40, + t0 = "2021-07-01", t1 = "2021-07-31"), + dx = 10, dy = 10, dt = "P1D", aggregation = "first", resampling = "near" + ) + cube <- gdalcubes::raster_cube( + fx$col, v, mask = gdalcubes::image_mask("SCL", values = c(3, 8, 9, 10, 11)) + ) + nc <- withr::local_tempfile(fileext = ".nc") + gdalcubes::write_ncdf(seen$pixel_fn(cube, seen$offset), nc) + pix <- terra::rast(nc) + expect_equal(names(pix), "red") + expect_equal(unique(terra::values(pix)[, 1]), 2) + dir <- drift:::cache_scheme_dir(cache, "sentinel-2-l2a") + files <- list.files(dir, all.files = TRUE, no.. = TRUE) + expect_length(files, 1L) + expect_match(files, "^count_[0-9a-f]{16}\\.tif$") + r <- out[[1]] + expect_equal(names(r), "red") + expect_equal(terra::datatype(terra::rast(file.path(dir, files))), "INT2U") + v <- terra::values(r)[, 1] + # zeros are NA; counts are the whole numbers written + expect_false(any(v == 0, na.rm = TRUE)) + expect_setequal(stats::na.omit(v), c(1, 2, 5)) + # served from the count cache on a second call, not rebuilt + testthat::local_mocked_bindings( + stac_cube_items = function(...) stop("rebuilt instead of served") + ) + again <- suppressMessages(dft_stac_composite( + aoi, years = 2022, months = 1, bands = "red", aggregation = "count", + mask_values = c(3L, 8L, 9L), cloud_cover_max = 55, resampling = "near", + cache_dir = cache + )) + expect_equal(terra::values(again[[1]]), terra::values(r)) +}) + +test_that("a count COG takes its overviews by nearest cell, so they are still counts (#92)", { + skip_if_not_installed("gdalcubes") + # Wider than one 512 px block, so the COG writer builds overviews. A 1/5 + # checkerboard averages to 3, which no cell holds; nearest keeps 1 or 5. + testthat::local_mocked_bindings( + stac_cube_items = function(...) list(features = list(), is_pre = FALSE), + stac_cube_assemble = function(...) { + r <- terra::rast(nrows = 1100, ncols = 1100, xmin = 0, xmax = 11000, + ymin = 0, ymax = 11000, crs = "EPSG:32609") + cells <- seq_len(terra::ncell(r)) + terra::values(r) <- ifelse( + (terra::rowFromCell(r, cells) + terra::colFromCell(r, cells)) %% 2 == 0, + 1, 5 + ) + names(r) <- "red" + r + } + ) + cache <- withr::local_tempdir() + suppressMessages(dft_stac_composite(aoi_pkg(), years = 2021, bands = "red", + aggregation = "count", cache_dir = cache)) + f <- list.files(drift:::cache_scheme_dir(cache, "sentinel-2-l2a"), + full.names = TRUE) + ov <- terra::rast(f, opts = "OVERVIEW_LEVEL=0") + expect_lt(terra::ncol(ov), 1100) + expect_true(all(terra::values(ov)[, 1] %in% c(1, 5))) +}) + +test_that("a count with no clear day anywhere drops the year, never caches zeros (#92)", { + skip_if_not_installed("gdalcubes") + aoi <- aoi_pkg() + aoi_t <- sf::st_transform(aoi, 32609) + testthat::local_mocked_bindings( + stac_cube_items = function(...) list(features = list(), is_pre = FALSE), + stac_cube_assemble = function(...) { + r <- terra::rast(terra::ext(aoi_t), resolution = 50, crs = "EPSG:32609", + vals = 0) + names(r) <- "red" + r + } + ) + cache <- withr::local_tempdir() + expect_error( + expect_warning( + suppressMessages(dft_stac_composite(aoi, years = 2021, bands = "red", + aggregation = "count", + cache_dir = cache)), + "no clear pixels" + ), + "No year produced" + ) + expect_length(list.files(cache, pattern = "\\.tif$", recursive = TRUE), 0L) +}) + +test_that("the count keys apart from every composite, and composite keys are unchanged (#92)", { + # pinned before #92: the family argument must not move an existing key + expect_identical(ckey(), "03ee8ecc66b832a8") + count <- ckey(aggregation = "count", dt = "P1D", family = "count") + # 0.19.x wrote reflectance under this key for aggregation = "count" + expect_false(identical(count, "08e0c5510e8ae297")) + expect_false(identical(count, ckey(aggregation = "count"))) + # the tag alone separates the families + expect_false(identical(ckey(family = "count"), ckey())) + expect_identical(count, ckey(aggregation = "count", dt = "P1D", + family = "count")) +}) From 18acaf6ef85c68dac8b1c787995febd95d1060ef Mon Sep 17 00:00:00 2001 From: almac2022 Date: Wed, 30 Sep 2026 10:16:32 -0700 Subject: [PATCH 04/10] Document the count and the aggregation set; hash aggregation as given (#92) Roxygen for dft_stac_composite() gains a section on counting clear days (distinct days, snow in the default mask, NA for no clear day, a lower bound where reads fail), the count_.tif cache name, the valid aggregation values, and a count example. dft_stac_cube() and dft_stac_fetch() list their valid values. aggregation_check() now returns the value as given rather than lower-cased: every caller hashes it into a cache key, and lower-casing moved the key of any mixed-case caller (v0.19.2 "Median" keyed 150c8ca5003bbe78), silently re-streaming cached cubes. Only the new count family is normalised. A caller-level test drives fetch, cube and composite with "Median" and records what reaches each key function. Adds a DRIFT_TEST_NETWORK test of a live count and a COG-overview test. Relates to #92 Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01PRhUJsuKABLfpBGktPoiBN --- R/dft_stac_composite.R | 50 +++++++-- R/dft_stac_cube.R | 16 +-- R/dft_stac_fetch.R | 5 +- inst/notes/gdalcubes-pc-gotchas.md | 2 +- man/dft_stac_composite.Rd | 45 +++++++- man/dft_stac_cube.Rd | 5 +- man/dft_stac_fetch.Rd | 5 +- planning/active/findings.md | 3 + planning/active/progress.md | 6 ++ planning/active/review-p345-round1.md | 104 +++++++++++++++++++ planning/active/review-p345-round2.md | 126 +++++++++++++++++++++++ planning/active/review-p345-round3.md | 91 ++++++++++++++++ planning/active/task_plan.md | 10 +- tests/testthat/test-dft_stac_composite.R | 75 +++++++++++++- tests/testthat/test-dft_stac_cube.R | 51 ++++++++- 15 files changed, 561 insertions(+), 33 deletions(-) create mode 100644 planning/active/review-p345-round1.md create mode 100644 planning/active/review-p345-round2.md create mode 100644 planning/active/review-p345-round3.md diff --git a/R/dft_stac_composite.R b/R/dft_stac_composite.R index e7e5fe9..5217e16 100644 --- a/R/dft_stac_composite.R +++ b/R/dft_stac_composite.R @@ -21,9 +21,28 @@ #' shared stretch to every composite it is given, so a brightness difference #' between years is a real one. #' +#' @section Counting clear observations: +#' `aggregation = "count"` returns, instead of reflectance, the number of +#' **distinct clear days** per pixel in each window. That is the number a caller +#' needs to choose composite windows. It is read one day per time step, so +#' where adjacent MGRS tiles overlap, their two items of one acquisition count +#' once, and a masked tile does not hide a clear one from the same day. "Clear" +#' means not in `mask_values` (by default cloud, cloud shadow, cirrus **and +#' snow**, so a spring or autumn count excludes snow cover as well as cloud), and +#' only scenes passing `cloud_cover_max` are counted. The counts are whole +#' numbers with no scale or offset applied, stored as integers, one layer per +#' band. The mask is shared, so bands count alike and `bands = "red"` is enough. +#' A pixel with no clear day is `NA`, never 0, because gdalcubes cannot tell a +#' chunk it failed to read from one that was all cloud. For the same reason a +#' day whose read failed for part of a chunk goes uncounted there, so where +#' reads fail a count is a lower bound; gdalcubes reports such failures only on +#' the console, as a composite's are. A count is not refused +#' across the Sentinel-2 2022-01-25 offset change, which does not affect it. +#' #' @section Caching: #' Each year's composite is written once under [dft_cache_path()] as -#' `/composite_.tif`, a Cloud Optimized GeoTIFF. The key hashes the +#' `/composite_.tif`, a Cloud Optimized GeoTIFF; a count is written +#' as `/count_.tif` and keyed apart from every composite. The key hashes the #' AOI geometry and every parameter that changes the pixels: bands (in order), #' months, the year's window, resolution, CRS, aggregation, resampling, cloud #' cover, mask values, reflectance scale and offset, `clip` and `tile_size`. @@ -68,7 +87,12 @@ #' @param crs Character. Target CRS as an EPSG string. When `NULL`, #' auto-detected from the AOI centroid's UTM zone. #' @param aggregation Character. How scenes within a window are reduced -#' (default `"median"`). +#' (default `"median"`): one of `"median"`, `"mean"`, `"min"`, `"max"`, +#' `"first"` or `"last"`, or `"count"` for the number of clear days per pixel +#' instead of reflectance (see the section on counting). Anything else is +#' refused. gdalcubes reads a value it does not know as no aggregation at all +#' and returns reflectance with no error, which is how `"count"` behaved in +#' drift 0.18.0 to 0.19.2, so drift passes it only values measured to work. #' @param resampling Character. Spatial resampling (default `"bilinear"`). #' @param clip Logical. Clip the output to the AOI polygon (default `FALSE`). #' Reference imagery is read around a place, not only inside it, so the @@ -78,7 +102,9 @@ #' (default 20, stricter than the cube's 60 because a composite of few clear #' scenes looks better than one of many cloudy ones). #' @param mask_values Integer vector of mask-band classes to exclude. When -#' `NULL`, taken from [dft_stac_config()]. +#' `NULL`, taken from [dft_stac_config()]; for Sentinel-2 that is the SCL +#' cloud, cloud-shadow, cirrus and **snow** classes, so snow cover is masked +#' as well as cloud. #' @param tile_size Numeric or `NULL` (default). Read-tiling edge length in CRS #' units; only tiles intersecting the AOI are streamed. A memory knob, not a #' speed one; see [dft_stac_cube()]. @@ -94,8 +120,10 @@ #' @return A named list of [terra::SpatRaster]s, one per year, each with one #' layer per band. Names label the window, e.g. `"2017 Jun–Jul"`, and #' become layer labels in [dft_map_interactive()]. Each raster's -#' [terra::time()] is the window start. A year with no usable scenes is -#' dropped with a warning; if every year is empty, the call aborts. +#' [terra::time()] is the window start. With `aggregation = "count"` each layer +#' holds integer clear-day counts rather than reflectance. A year with no +#' usable scenes (or, for a count, no clear day anywhere) is dropped with a +#' warning; if every year is empty, the call aborts. #' #' @seealso [dft_map_interactive()] (`rgb =`) to display them, #' [dft_stac_cube()] for index time series. @@ -121,6 +149,11 @@ #' }) #' #' dft_map_interactive(rgb = c(tc, fc[1]), aoi = aoi) +#' +#' # How many clear days each July had, to choose a composite window +#' n <- dft_stac_composite(aoi, years = 2017:2023, months = 7, bands = "red", +#' aggregation = "count") +#' sapply(n, function(r) terra::global(r, "max", na.rm = TRUE)[[1]]) #' } #' #' @export @@ -152,7 +185,10 @@ dft_stac_composite <- function(aoi, } aggregation <- aggregation_check(aggregation, c(.cube_view_aggregations, "count")) - is_count <- identical(aggregation, "count") + # matched without case, like every aggregation; the count family is new, so + # normalising it moves no existing key + is_count <- identical(tolower(aggregation), "count") + if (is_count) aggregation <- "count" family <- if (is_count) "count" else "composite" years <- composite_years_check(years) months <- composite_months_check(months) @@ -487,7 +523,7 @@ composite_label <- function(year, months) { #' order-sensitive (they are the layer order); months and mask values are not. #' #' `family = "count"` is the clear-observation count (#92). It keys apart from -#' every composite, and in particular from the 0.19.x `composite_.tif` +#' every composite, and in particular from the 0.18.0-0.19.2 `composite_.tif` #' files written under `aggregation = "count"`, which hold reflectance: the #' tag, and the `P1D` read step the caller passes, both change the hash. The #' default keeps every existing composite key unchanged. diff --git a/R/dft_stac_cube.R b/R/dft_stac_cube.R index a122c5c..6ce3b6f 100644 --- a/R/dft_stac_cube.R +++ b/R/dft_stac_cube.R @@ -36,7 +36,10 @@ #' (default `"P1M"`, monthly). The cadence [dft_rast_break()]'s `frequency` #' must agree with. #' @param aggregation Character. Temporal aggregation for multiple scenes in one -#' `dt` window (default `"median"`). +#' `dt` window (default `"median"`): one of `"median"`, `"mean"`, `"min"`, +#' `"max"`, `"first"` or `"last"`. Anything else is refused: gdalcubes reads a +#' value it does not know as no aggregation, with no error, so drift passes it +#' only values measured to work. #' @param resampling Character. Spatial resampling (default `"bilinear"`). #' @param clip Logical. When `TRUE` (default), clip the returned stack to the AOI #' polygon with `terra::mask()`, so @@ -302,8 +305,8 @@ dft_stac_cube <- function(aoi, # round-tripping each value through cube_view()$aggregation: the six below # survive. It also honours "count_values" and "count_images", which are held back # on purpose: they count items, so overlapping MGRS tiles double-count one -# acquisition, and every caller here would then apply a reflectance scale or an -# index expression to a count. "none" is the silent fallback itself. +# acquisition, and every caller here would then read the count as reflectance, +# an index or a class code. "none" is the silent fallback itself. .cube_view_aggregations <- c("min", "max", "mean", "median", "first", "last") #' Refuse an `aggregation` outside `allowed`, before any network call @@ -311,8 +314,9 @@ dft_stac_cube <- function(aoi, #' Shared by [dft_stac_fetch()], [dft_stac_cube()] and [dft_stac_composite()], #' the three callers of `gdalcubes::cube_view()`. An unsupported value is not #' refused downstream, so it has to be refused here (#92). Matching ignores case, -#' as gdalcubes does, and the value is returned lower-cased so `"Median"` and -#' `"median"` share a cache key. +#' as gdalcubes does, and the value is returned exactly as given: every caller +#' hashes it into a cache key, and lower-casing it would move the key of any +#' mixed-case caller, silently re-streaming a cube that is already cached. #' @noRd aggregation_check <- function(aggregation, allowed = .cube_view_aggregations) { if (!is.character(aggregation) || length(aggregation) != 1L || @@ -328,7 +332,7 @@ aggregation_check <- function(aggregation, allowed = .cube_view_aggregations) { } ), class = "drift_bad_aggregation") } - tolower(aggregation) + aggregation } diff --git a/R/dft_stac_fetch.R b/R/dft_stac_fetch.R index 5a0b3ea..9a6db78 100644 --- a/R/dft_stac_fetch.R +++ b/R/dft_stac_fetch.R @@ -48,7 +48,10 @@ #' @param dt Character. ISO 8601 duration for the temporal aggregation window #' (default `"P1Y"`). #' @param aggregation Character. Temporal aggregation method (default -#' `"first"`). Use `"median"` for multi-scene composites. +#' `"first"`). Use `"median"` for multi-scene composites. One of `"first"`, +#' `"last"`, `"median"`, `"mean"`, `"min"` or `"max"`; anything else is +#' refused: gdalcubes reads a value it does not know as no aggregation, with +#' no error, so drift passes it only values measured to work. #' @param resampling Character. Spatial resampling method (default `"near"` #' for categorical data). #' @param tile_size Numeric or `NULL` (default). Edge length, in CRS units diff --git a/inst/notes/gdalcubes-pc-gotchas.md b/inst/notes/gdalcubes-pc-gotchas.md index c4e0422..90672a4 100644 --- a/inst/notes/gdalcubes-pc-gotchas.md +++ b/inst/notes/gdalcubes-pc-gotchas.md @@ -216,7 +216,7 @@ Measured on gdalcubes 0.7.5, rstac 1.0.1 and terra 1.9.50, with the BULK floodpl - **`cube_view()` does not refuse an aggregation it does not know. It reads it as `"none"`.** The R wrapper checks only that the value is one string. The C++ side lower-cases it, and maps anything unrecognised to `AGG_NONE`, which copies every image in. That copy includes NaNs, so a masked item can blank a clear one. This is how `aggregation = "count"` came back as plausible reflectance. - Measured by round-tripping values through `cube_view()$aggregation`. `min`, `max`, `mean`, `median`, `first` and `last` survive, and so do the undocumented `count_values` and `count_images`. `count`, `sum` and `""` become `none`. - - drift passes only the first six (`aggregation_check()`). `count_*` count **items**, so overlapping MGRS tiles double-count one acquisition, and every drift caller would scale the result as reflectance. + - drift passes only the first six (`aggregation_check()`). `count_*` count **items**, so overlapping MGRS tiles double-count one acquisition, and every drift caller would read the result as reflectance, an index or a class code. - `resampling` has the same fallback, to `near` (`bilinaer -> near`), in #96. - **Count clear days, not items: `dt = "P1D"`, then `reduce_time(cube, "count(B04)", names = ...)`.** Masked pixels are NaN before aggregation, and every aggregation tried (first, max, median) skips NaN within a day. So two same-day tiles give one clear day, and a cloudy tile does not blank a clear one. The string reducer is C++, so the R-callback closure trap above does not apply. - `reduce_time()` passes a single-time-step cube through unchanged. A count over a one-day window is therefore the input, not a count. diff --git a/man/dft_stac_composite.Rd b/man/dft_stac_composite.Rd index 5892198..90b8607 100644 --- a/man/dft_stac_composite.Rd +++ b/man/dft_stac_composite.Rd @@ -49,7 +49,12 @@ colour. Layers are named by role, in the order given.} auto-detected from the AOI centroid's UTM zone.} \item{aggregation}{Character. How scenes within a window are reduced -(default \code{"median"}).} +(default \code{"median"}): one of \code{"median"}, \code{"mean"}, \code{"min"}, \code{"max"}, +\code{"first"} or \code{"last"}, or \code{"count"} for the number of clear days per pixel +instead of reflectance (see the section on counting). Anything else is +refused. gdalcubes reads a value it does not know as no aggregation at all +and returns reflectance with no error, which is how \code{"count"} behaved in +drift 0.18.0 to 0.19.2, so drift passes it only values measured to work.} \item{resampling}{Character. Spatial resampling (default \code{"bilinear"}).} @@ -63,7 +68,9 @@ default keeps the whole AOI bounding box. \code{TRUE} uses the same rule as scenes looks better than one of many cloudy ones).} \item{mask_values}{Integer vector of mask-band classes to exclude. When -\code{NULL}, taken from \code{\link[=dft_stac_config]{dft_stac_config()}}.} +\code{NULL}, taken from \code{\link[=dft_stac_config]{dft_stac_config()}}; for Sentinel-2 that is the SCL +cloud, cloud-shadow, cirrus and \strong{snow} classes, so snow cover is masked +as well as cloud.} \item{tile_size}{Numeric or \code{NULL} (default). Read-tiling edge length in CRS units; only tiles intersecting the AOI are streamed. A memory knob, not a @@ -85,8 +92,10 @@ replacement is atomic.} A named list of \link[terra:SpatRaster-class]{terra::SpatRaster}s, one per year, each with one layer per band. Names label the window, e.g. \code{"2017 Jun–Jul"}, and become layer labels in \code{\link[=dft_map_interactive]{dft_map_interactive()}}. Each raster's -\code{\link[terra:time]{terra::time()}} is the window start. A year with no usable scenes is -dropped with a warning; if every year is empty, the call aborts. +\code{\link[terra:time]{terra::time()}} is the window start. With \code{aggregation = "count"} each layer +holds integer clear-day counts rather than reflectance. A year with no +usable scenes (or, for a count, no clear day anywhere) is dropped with a +warning; if every year is empty, the call aborts. } \description{ Reference imagery for windows you choose. For each year, streams every scene @@ -111,10 +120,31 @@ contrast stretching to display time. \code{\link[=dft_map_interactive]{dft_map_i shared stretch to every composite it is given, so a brightness difference between years is a real one. } +\section{Counting clear observations}{ + +\code{aggregation = "count"} returns, instead of reflectance, the number of +\strong{distinct clear days} per pixel in each window. That is the number a caller +needs to choose composite windows. It is read one day per time step, so +where adjacent MGRS tiles overlap, their two items of one acquisition count +once, and a masked tile does not hide a clear one from the same day. "Clear" +means not in \code{mask_values} (by default cloud, cloud shadow, cirrus \strong{and +snow}, so a spring or autumn count excludes snow cover as well as cloud), and +only scenes passing \code{cloud_cover_max} are counted. The counts are whole +numbers with no scale or offset applied, stored as integers, one layer per +band. The mask is shared, so bands count alike and \code{bands = "red"} is enough. +A pixel with no clear day is \code{NA}, never 0, because gdalcubes cannot tell a +chunk it failed to read from one that was all cloud. For the same reason a +day whose read failed for part of a chunk goes uncounted there, so where +reads fail a count is a lower bound; gdalcubes reports such failures only on +the console, as a composite's are. A count is not refused +across the Sentinel-2 2022-01-25 offset change, which does not affect it. +} + \section{Caching}{ Each year's composite is written once under \code{\link[=dft_cache_path]{dft_cache_path()}} as -\verb{/composite_.tif}, a Cloud Optimized GeoTIFF. The key hashes the +\verb{/composite_.tif}, a Cloud Optimized GeoTIFF; a count is written +as \verb{/count_.tif} and keyed apart from every composite. The key hashes the AOI geometry and every parameter that changes the pixels: bands (in order), months, the year's window, resolution, CRS, aggregation, resampling, cloud cover, mask values, reflectance scale and offset, \code{clip} and \code{tile_size}. @@ -167,6 +197,11 @@ chips <- lapply(seq_len(nrow(buf)), function(i) { }) dft_map_interactive(rgb = c(tc, fc[1]), aoi = aoi) + +# How many clear days each July had, to choose a composite window +n <- dft_stac_composite(aoi, years = 2017:2023, months = 7, bands = "red", + aggregation = "count") +sapply(n, function(r) terra::global(r, "max", na.rm = TRUE)[[1]]) } } diff --git a/man/dft_stac_cube.Rd b/man/dft_stac_cube.Rd index 4d13fad..325b671 100644 --- a/man/dft_stac_cube.Rd +++ b/man/dft_stac_cube.Rd @@ -47,7 +47,10 @@ from the AOI centroid's UTM zone.} must agree with.} \item{aggregation}{Character. Temporal aggregation for multiple scenes in one -\code{dt} window (default \code{"median"}).} +\code{dt} window (default \code{"median"}): one of \code{"median"}, \code{"mean"}, \code{"min"}, +\code{"max"}, \code{"first"} or \code{"last"}. Anything else is refused: gdalcubes reads a +value it does not know as no aggregation, with no error, so drift passes it +only values measured to work.} \item{resampling}{Character. Spatial resampling (default \code{"bilinear"}).} diff --git a/man/dft_stac_fetch.Rd b/man/dft_stac_fetch.Rd index 9f7406b..e5aad4f 100644 --- a/man/dft_stac_fetch.Rd +++ b/man/dft_stac_fetch.Rd @@ -46,7 +46,10 @@ When \code{NULL}, auto-detected from the AOI centroid's UTM zone.} (default \code{"P1Y"}).} \item{aggregation}{Character. Temporal aggregation method (default -\code{"first"}). Use \code{"median"} for multi-scene composites.} +\code{"first"}). Use \code{"median"} for multi-scene composites. One of \code{"first"}, +\code{"last"}, \code{"median"}, \code{"mean"}, \code{"min"} or \code{"max"}; anything else is +refused: gdalcubes reads a value it does not know as no aggregation, with +no error, so drift passes it only values measured to work.} \item{resampling}{Character. Spatial resampling method (default \code{"near"} for categorical data).} diff --git a/planning/active/findings.md b/planning/active/findings.md index 66641f0..d0de5ed 100644 --- a/planning/active/findings.md +++ b/planning/active/findings.md @@ -140,6 +140,9 @@ Six round-3 mutations each turn a test red: overview AVERAGE, `cfg$mask_values`, | Error | Resolution | |-------|------------| +| P3-5 round 3: the `"Median"` pin called the key function directly, so lowering inside any caller stayed green; the "3 mutations red" claim overstated it | Caller-level test drives fetch/cube/composite with `"Median"` and records what reaches each key function. Mutations lowering in each caller, and in the helper, all turn red. Round 3 computed every non-count key against v0.19.2 (aggregation in any case, months, mask_values, tile_size, clip, bands, dt, res, resampling, cloud_cover_max): none moved. With every key flow tested, the enumeration ends the loop. Cost: 1 plan review + 9 code-check rounds | +| P3-5 round 2: the P1 round-1 fix lower-cased `aggregation`, which moved the cache key of every mixed-case caller (v0.19.2 `"Median"` keyed `150c8ca5003bbe78`, and the fix made it `03ee8ecc66b832a8`). For a cube that means a silent multi-hour re-stream plus an orphan `superseded` cannot reclaim. Also, the "every caller would scale" reason is false for `dft_stac_fetch()` | `aggregation_check()` returns the value as given, and only the new count family is normalised. Enumeration: 3 callers x (cache key, `cube_view`) = 6 flows, all hashing or passing the value as given, the same as v0.19.2, and gdalcubes lower-cases on its own side. Pinned the `"Median"` key; 3 mutations turn red. Reason reworded in NEWS, the gotchas note and the code comment | +| P3-5 round 1, all in the prose. 1: the docs and NEWS said drift refuses what "gdalcubes would not honour", which is false for `count_*` (the P1 round 3 claim, back again). 2: the bug and stale files dated "0.19.x" when it was 0.18.0. 3: the reclaim remedy would delete cubes. 4: "every chip 5-9" (only the per-chip max was recorded). 5: the RSS evidence was only in a gitignored log. 6: an unmeasured causal claim | Reworded as drift policy, with `count_*` named. Changed to 0.18.0-0.19.2 (verified with `git show v0.18.0`). Remedy now names `composite_*.tif` under `v2/sentinel-2-l2a` and spares cubes. Changed to "each chip's highest count". Committed `time_chips_count.txt`. Causal claim softened, with the uncontrolled comparison stated | | P2 round 3: count overviews untested (fixture below one 512 px block); `mask_values`/`months`/query/`resampling` reached the mocks unobserved | 1100 px checkerboard overview test; count test captures and asserts the four with non-default values; each of 6 mutations red. Enumeration above ends the loop | | P2 round 1: no test drove `dft_stac_composite()`'s count `pixel_fn` (a reflectance closure there stayed green) | The mocked assemble now captures `pixel_fn` and runs it on a fixture cube, expecting 2 clear days named `red`. The mutation turns 2 tests red | | Round 3: the error message claimed gdalcubes "would not honour" values it does honour (none, count_*) | Reworded to state the fallback and drift's policy. Enumerated all 8 gdalcubes claims in the Phase 1 diff, all measured or test-asserted. Sibling defect on `resampling` filed as #96 | diff --git a/planning/active/progress.md b/planning/active/progress.md index c2844db..ec94be9 100644 --- a/planning/active/progress.md +++ b/planning/active/progress.md @@ -10,3 +10,9 @@ - Phase 1 committed (727c5d9) after 3 code-check rounds. Round 1 found `last` refused and was fixed by a behaviour round trip. Round 2 was clean. Round 3's message claim was fixed and then enumerated (8 claims). - Filed #96 (resampling has the same silent fallback) - Phase 2 implemented. 7 mutations each redden the composite tests. +- Phase 2 committed (ba0b203) after 3 code-check rounds and an enumeration of the count build's wiring +- Phases 3-5 reviewed over 3 rounds: + - round 1: six prose claims fixed; + - round 2: the Phase 1 lower-casing had moved mixed-case cache keys, fixed by returning the value as given; + - round 3: a caller-level key test was added; the enumeration against v0.19.2 keys shows none moved. +- Live checks: counts 4-9 on a 2 km square. BULK: 20 count chips in 9.6 min, 0.46 GiB peak. Full suite: 1431 pass, 0 fail. diff --git a/planning/active/review-p345-round1.md b/planning/active/review-p345-round1.md new file mode 100644 index 0000000..de122f8 --- /dev/null +++ b/planning/active/review-p345-round1.md @@ -0,0 +1,104 @@ +# Review p345 round 1: #92 docs, network test, BULK count benchmark, NEWS + +Snapshot: scratchpad/p345snap (staged index). Offline suite on the snapshot: +`test-dft_stac_composite.R` [ FAIL 0 | WARN 0 | SKIP 3 | PASS 99 ]. + +## Findings + +- **[bug: false claim] NEWS.md:5, R/dft_stac_composite.R:92-94, R/dft_stac_cube.R:39-41, R/dft_stac_fetch.R:51-53 (and the three man/*.Rd).** + NEWS says drift now refuses any `aggregation` "outside the six gdalcubes actually honours". The three + `@param aggregation` docs say anything else is refused "because gdalcubes would silently read it as no + aggregation". Both are false for `count_values` and `count_images`. gdalcubes 0.7.5 honours both + (re-probed here: `cube_view(..., aggregation = "count_values")$aggregation` returns `"count_values"`, + and `"count_images"` returns `"count_images"`). drift refuses them **by policy**, as its own comment + at R/dft_stac_cube.R:305-308 says. Phase 1 round 3 reworded the error message to fix exactly this + claim (findings.md, Errors Encountered, "claimed gdalcubes would not honour values it does honour"), + and the docs and NEWS now bring it back. Suggested fix: "the six drift passes to gdalcubes (it also + honours `count_values`/`count_images`, which count items and are held back)", and in the params: + "anything else is refused; an unknown value would be read by gdalcubes as no aggregation". + +- **[bug: false claim, understates scope] NEWS.md:12 ("Cached 'count' files from 0.19.x"), R/dft_stac_composite.R:94 ("how `"count"` behaved in drift 0.19.x"), tests/testthat/test-dft_stac_composite.R:~573 comment ("0.19.x wrote reflectance under this key").** + The defect and its stale cache files date from **0.18.0**, not 0.19.x. `dft_stac_composite()` was + introduced in 0.18.0 (#79, NEWS.md:41). In `git show v0.18.0:R/dft_stac_composite.R`, `aggregation` + already goes straight to `stac_cube_assemble()` (line 212), and the file is already written as + `composite_.tif` (line 194). The file has no key or aggregation diff between v0.18.0 and + v0.19.2, and the cache scheme has been `v2` since 0.12.0. A reader on 0.18.0, which is where the + floodplains#93 work started, would read the entry as not applying to them, both for the wrong + values and for the cache files. It should say "0.18.0-0.19.2". + +- **[fragile: remedy over-deletes / is not actionable] NEWS.md:12, "clear the current scheme or delete them by hand".** + The `superseded` part is correct. `superseded` = `setdiff(list.dirs(base), /v2)` (R/dft_cache.R), + and these files are in `v2`, so it does not reclaim them. + - **"Clear the current scheme"** means `dft_cache_clear(scheme = "current")`. That deletes every + current entry for every source, including `dft_stac_cube()`'s `cube_*.tif`. The code calls those a + "multi-hour Sentinel-2 stream" (R/dft_stac_cube.R:291). Even with `source = "sentinel-2-l2a"`, + `unlink()` removes the whole directory, cubes included. + - **"Delete them by hand"** gives no way to find them, and the same sentence says they "cannot be + told from real composites". + - **A remedy that works and spares cubes:** delete `composite_*.tif` under + `/v2/sentinel-2-l2a/` (composites rebuild cheaply; `cube_*` and the new `count_*` are + untouched). + +- **[bug: false claim] NEWS.md:13, "Every chip counted 5-9 clear days, with no `NA` cells".** + `chips_count_per_chip.csv` records only `count_max` per chip (range 5-9, recomputed). It does not + record a per-pixel minimum. The true statement is "each chip's **maximum** count was 5-9". As + written, it reads as every pixel having at least 5 clear days, and that is exactly the number a + reader uses to choose a window. The live 2 km square had a minimum of 4-5, and nothing bounds the + BULK chips' minimum. The "no NA cells" half is correct: `count_na` is 0 in all 20 rows. + +- **[fragile: evidence not committed] NEWS.md:13, "0.46 GiB peak (`/usr/bin/time -l`) ... Evidence is in `data-raw/logs/benchmark_composite_bulk/*count*`".** + The committed `*count*` files are `chips_count_per_chip.csv` and `timings_chips-count.csv`, and + neither holds the peak. The 489,013,248 B maximum resident set size and the 577.48 s real are only + in `run_chips_count.log`. That file exists in the working tree and is gitignored + (`data-raw/logs/**/*.log`). So on any other checkout the glob does not contain the evidence for the + memory figure. #79's 0.48 GiB, by contrast, has `rss_chips.txt` committed. The gitignore rule is + accepted, so the claim is the issue: commit the `/usr/bin/time` block as, for example, + `time_chips-count.txt`, or drop "Evidence is in" for the peak. + +- **[low: unmeasured causal claim] NEWS.md:13, "The count reads daily steps but only one band, so it is the cheaper of the two."** + That the count was cheaper (17.5 vs 9.6 min on the same 20 chips) is measured. The reason given + is inferred. No one-band median control was run, the two runs were two days apart (2026-09-28 vs + 09-30), and the count run was network-bound (44.5 s user of 577 s real). The 0.46 vs 0.48 GiB + comparison also sets a true peak (`time -l`) against a 2 s sampler, which reads low (CLAUDE.md). + NEWS does not claim the count uses less memory, so this is only a caution. Either soften the + wording to "likely because" or drop the clause. + +## Verified correct (recomputed from the artifacts) + +- **Count timings** (`chips_count_per_chip.csv` / `timings_chips-count.csv`): + - 9.6 min: loop 575.06 s = 9.58 min, or `time -l` real 577.48 s = 9.62 min. + - Median 24.7 s (24.696), max 59.6 s. + - Count max 5-9, 0 NA cells, 3,660 cells a chip. + - 0.46 GiB = 489,013,248 B / 2^30 = 0.455, from the local log. +- **The #79 median chips** (`chips_per_chip.csv`): + - Rows 1-20 sum to 1,051.62 s = 17.5 min, median 44.37 s. + - The 100-chip peak is 503,296 KiB = 0.48 GiB (`rss_chips.txt` max). + - The same seed-79 draw, so chip i is the same geometry, cells 3,660. +- **Live numbers match findings.md:** 20 items / 12 dates / 5-9, and 6 items / 4 dates / 4 everywhere. + The issue's reflectance was 0.0077-0.1835, which NEWS gives as 0.008-0.18. +- **`?cube_view` (0.7.5) omits `last`:** Rd lists "min", "max", "mean", "median", or "first". `last` + round-trips. +- **Case:** `Median` becomes `median`. +- **`resampling`:** "foo" becomes `near`, and #96's title matches. +- **SCL default `c(3, 8, 9, 10, 11)`:** 11 is snow/ice (R/dft_stac_config.R:67-71). +- **The code matches every count-behaviour claim in the docs:** + - INT2U, and the `count_` prefix. + - 0 becomes NA (`count_zero_na`) before `cube_check_nonempty`. + - No offset split (`is_pre[] <- FALSE`). + - A P1D read with `first`, and NEAREST overviews. + - An empty year is dropped with "no clear pixels". + - The composite key is pinned. +- **Network test** (tests/testthat/test-dft_stac_composite.R:580-617): it would fail on reflectance. + - Values of about 0.03 fail the whole-number check, `min >= 1` and `median >= 1`. + - Reflectance truncated to 0 by INT2U fails `min >= 1`. + - An all-NA raster fails `expect_gte(median(numeric(0)), 1)`. Checked: an NA comparison is an + expectation failure in testthat 3.3.2. + - The bound `max <= distinct dates` is sound, because drift reads only items from the same + intersects/cloud-cover/datetime query. + - `rstac` is in Imports and `withr` in Suggests. +- **Benchmark stage `chips-count`** (data-raw/benchmark_composite_bulk.R:82-114) runs as written: + - It uses the same `set.seed(79)` / `st_sample` draw, then `[1:20, ]`. + - It writes both CSVs, and the shared tail writes `timings_.csv` and "ALL STAGES DONE". + - The local log confirms a completed run. + - Cosmetic only: `benchmark_composite_bulk-run.sh`'s usage line and summary grep do not know + `chips-count`, but the header routes that stage through `/usr/bin/time -l` instead. diff --git a/planning/active/review-p345-round2.md b/planning/active/review-p345-round2.md new file mode 100644 index 0000000..def982c --- /dev/null +++ b/planning/active/review-p345-round2.md @@ -0,0 +1,126 @@ +# Review p345 round 2: #92 docs, network test, BULK count benchmark, NEWS + +Snapshot: scratchpad/p345snap (staged index). Read-only review. Nothing under the repo was edited +except this file. + +## Findings + +- **[low: false claim] NEWS.md:5, "drift refuses both, because ... every caller would scale the result as reflectance".** + This is false for `dft_stac_fetch()`, one of the three callers the same bullet names. Fetch + applies no scale: `grep scale R/dft_stac_fetch.R` finds only prose, and the sources it serves + are categorical. There, a `count_values` result would be read as **class codes** by + `dft_rast_classify()`, not scaled as reflectance. The policy is still right, but the reason given + does not hold for fetch. The same sentence also appears in committed text outside this diff: + `inst/notes/gdalcubes-pc-gotchas.md:219` ("every drift caller would scale the result as + reflectance") and the comment at `R/dft_stac_cube.R:306-308` ("a reflectance scale or an index + expression"). Suggested wording: "every caller would read the result as reflectance, an index + or a class code". Fix it in all three places, not only NEWS. This is round 1's mechanism: a + reason recalled from the composite/cube path and applied to every caller without reading fetch. + +- **[low: false claim] NEWS.md:12, "every genuine composite key is unchanged (pinned in the tests)".** + This holds only when `aggregation` is lower case. `aggregation_check()` now returns + `tolower(aggregation)`, and that lowered value is the one hashed. I computed this against the + v0.19.2 tree (`git archive v0.19.2`, `pkgload::load_all`, and the test file's own `key_args()`): + - Default: `03ee8ecc66b832a8`, which matches the pin. + - `aggregation = "count"`: `08e0c5510e8ae297`, which matches the test. + - `aggregation = "Median"`: **`150c8ca5003bbe78`**. The same call now keys as + `03ee8ecc66b832a8`. + + So a composite cached under a mixed-case value moves key and is rebuilt. The old file becomes an + orphan inside `v2`, where `dft_cache_clear(scheme = "superseded")` cannot reclaim it. + `stac_cube_cache_key()` (R/dft_stac_cube.R:224-225) and `stac_cache_key()` (fetch) also hash + `aggregation`, so a mixed-case `dft_stac_cube()` caller silently re-streams what the code calls a + "multi-hour Sentinel-2 stream". The test pins only the lower-case default, so it cannot see this. + Mixed-case callers are rare, and sharing a key afterwards is the right behaviour. The fault is + the unqualified claim. Suggested wording: "every composite key for a lower-case `aggregation` is + unchanged; a mixed-case value (`"Median"`) now keys as its lower-case form, so a cube or + composite cached under one is rebuilt once". Add this near "Case is ignored" on line 5. + +## Round-1 fixes: verified + +- **The remedy path is exact.** `cache_scheme_dir(cache_dir, source)` is + `file.path(dft_cache_path(cache_dir), "v2", source)` (R/dft_cache.R:98-106). `dft_stac_composite()` + writes to `cache_scheme_dir(cache_dir, source)` (line 210), and the only `cube = TRUE` source in + `dft_stac_config()` is `sentinel-2-l2a`, so every composite is under + `/v2/sentinel-2-l2a/`. The scheme was already `v2` at v0.18.0 + (`git show v0.18.0:R/dft_cache.R`:106). +- **Deleting `composite_*.tif` spares `cube_*.tif` and `count_*.tif`.** The cube writes + `paste0("cube_", key, ".tif")` (R/dft_stac_cube.R:230). The count writes + `paste0(family, "_", key, ".tif")` with family `"count"`. v0.18.0 already stripped `time` before + the write (line 226), so no `.aux.json` sidecars are left behind the glob. The local cache holds + only `composite_*.tif`. +- **"0.18.0-0.19.2" is correct.** v0.17.0 has no `R/dft_stac_composite.R`. v0.18.0, v0.19.0 and + v0.19.2 all pass `aggregation` straight through (line 212) and write `composite_.tif` + (line 194). The pinned stale key `08e0c5510e8ae297` is what v0.19.2 computes for + `aggregation = "count"` (computed above). +- **"superseded does not reclaim them" / "current would delete cubes too" is correct.** + `superseded` = `setdiff(list.dirs(base), /v2)`, and `current` unlinks all of `/v2`. +- **"refused by policy" is correct.** The NEWS wording no longer says gdalcubes cannot honour + `count_*`. The error message and the three `@param` texts say "reads a value it does not know as + no aggregation", which is true. +- **The per-chip max, the committed RSS evidence and the softened causal claim are all fixed.** + +## Claim enumeration: NEWS entry + +| # | Claim | Source | OK | +|---|---|---|---| +| 1 | `aggregation` went straight into `cube_view()` | v0.18.0-v0.19.2 line 212 | yes | +| 2 | Unknown value read as no aggregation, no error, copies every scene; masked can overwrite clear | gotchas.md:217-218 (round trip measured). "Copies incl. NaN" is from the C++ reading, hedged "could" | yes | +| 3 | Looked like red reflectance 0.008-0.18 | findings.md:21 (0.0077-0.1835) | yes | +| 4 | Three functions refuse all but the six | aggregation_check() at fetch:105, cube:169, composite:186 | yes | +| 5 | Each of the six round-trips; `?cube_view` omits `last` | test-dft_stac_cube.R:778-; findings Errors row | yes | +| 6 | gdalcubes honours `count_values`/`count_images` | gotchas.md:218, round 1 re-probe | yes | +| 7 | They count items; **every caller would scale as reflectance** | fetch applies no scale | **no (finding 1)** | +| 8 | Case ignored, as gdalcubes ignores it | gotchas.md:217; aggregation_check | yes | +| 9 | Check runs before any network call | first statements after check_gdalcubes/config in all three | yes | +| 10 | `resampling` falls back to `near`, #96 | round 1 probe, #96 title | yes | +| 11 | count = distinct clear days per pixel | P1D + first + reduce_time count; offline probe findings:59 | yes | +| 12 | Overlapping tiles count once | findings:88, fixture test 371- | yes | +| 13 | 20 items/12 dates -> 5-9; 6 items/4 dates -> 4 everywhere | findings live table (cc_max 100 / 20 rows) | yes (the conditions that separate the two runs are not stated, but nothing in the sentence is false) | +| 14 | Integer, no scale/offset, INT2U COG `count_.tif` | composite writeRaster datatype, cache_file | yes | +| 15 | NA never 0; 0/NA follows chunking/`parallel`; failed chunk NA | count_zero_na doc (measured 12,544 / 4,352+8,192); gotchas.md:206 | yes | +| 16 | Year with no clear day dropped with warning | count_zero_na before cube_check_nonempty -> drift_empty_cube -> composite_skip; test at ~555 | yes | +| 17 | Count spans the offset change without refusal | `is_pre[] <- FALSE`, no composite_offset_check | yes | +| 18 | Default mask includes snow + cloud/shadow/cirrus | dft_stac_config.R:71 `c(3,8,9,10,11)` | yes | +| 19 | Since 0.18.0 cached reflectance as `composite_.tif` | v0.18.0 tree | yes | +| 20 | Count keys under own tag + prefix; nothing reads old files | family tag in key, `count_` prefix, is_count forces family | yes | +| 21 | **Every genuine composite key unchanged** | lower-casing moves mixed-case keys | **no (finding 2)** | +| 22 | In `v2`, superseded does not reclaim | dft_cache.R | yes | +| 23 | Remedy path, spares cubes; `current` deletes cubes | above | yes | +| 24 | 20 chips, 300 m, Jul-Aug 2023, red | benchmark script chips-count stage | yes | +| 25 | 9.6 min, median 24.7 s, 0.46 GiB | recomputed: 575.06 s loop / 577.48 s real; 24.696; 489,013,248 B = 0.455 GiB | yes | +| 26 | Each chip's highest count 5-9, no NA | CSV count_max 5..9, count_na all 0 | yes | +| 27 | Same 20 as #79: 17.5 min, median 44.4 s | chips_per_chip.csv rows 1-20: 1,051.6 s, 44.37 | yes | +| 28 | #79 0.48 GiB over 100 chips, sampled every 2 s | rss_chips.txt max 503,296 KiB; run.sh `sleep 2` | yes | +| 29 | Count reads one band vs three | chips stage default bands = RGB; count bands = "red" | yes | +| 30 | Runs two days apart | wallclock_chips.txt 2026-09-28, time_chips_count.txt 2026-09-30 | yes | +| 31 | Evidence in `*count*` | three committed files hold claims 24-26. Claims 27-28 are in chips_per_chip.csv / rss_chips.txt (committed, not matched by the glob) | acceptable | + +## Claim enumeration: new roxygen (and man/, which matches it) + +| Claim | Source | OK | +|---|---|---| +| One day per time step; overlapping tiles once; masked tile does not hide a clear same-day one | findings:59 probe (first/max/median skip NaN); fixture test | yes | +| Snow in default mask; `cloud_cover_max` applies | config:71; stac_cube_items filter | yes | +| Whole numbers, no scale/offset, integers, one layer per band | pixel_fn, INT2U | yes | +| Mask shared, so bands count alike | image_mask applies to every band. Band-level nodata could differ at edges, but that is marginal and `bands = "red"` advice is sound | yes | +| NA never 0, because failed chunk = NaN | count_zero_na doc; gotchas:206 | yes | +| Partial read failure -> uncounted -> lower bound; stderr only | inference from gotchas:206 (failed chunks written NA, reported on stderr). The NA day is skipped by count() | yes (inferred, not measured, and stated as a consequence) | +| Not refused across 2022-01-25 | code | yes | +| `count_.tif`, keyed apart | code | yes | +| @param: the six or count; else refused; 0.18.0-0.19.2 behaviour | code, tags | yes | +| @param mask_values: SCL cloud, shadow, cirrus, snow | config:71 | yes | +| @return: integer counts; no clear day anywhere -> dropped | code | yes | +| Example: `aoi` defined earlier in the same `\dontrun` block; valid call | examples | yes | +| cube/fetch @param: six, else refused | aggregation_check | yes | + +## Notes (not findings) + +- **The test title at `tests/testthat/test-dft_stac_composite.R:299` still reads "refuses an aggregation gdalcubes would not honour".** + This is the round-1 wording, left in committed code outside this diff. It is not user-facing, + but the same `grep` that found round 1 finds it. +- **"which rebuild on their next call" is true, but says nothing of the cost.** + A floodplain-wide composite is hours, and with #88 open a BULK-size one does not complete. In + practice nobody is likely to hold such a cache, so this is not a false claim. +- **The network test mirrors drift's own query:** the same `stac_search(intersects = sfg)`, + `ext_filter`, `limit = 500` and `items_fetch()`. `aoi_pkg()` exists at test file line 3. diff --git a/planning/active/review-p345-round3.md b/planning/active/review-p345-round3.md new file mode 100644 index 0000000..218bf42 --- /dev/null +++ b/planning/active/review-p345-round3.md @@ -0,0 +1,91 @@ +# Review p345 round 3: #92 key/filename flows, the case-preserving fix, anything missed + +Snapshot: scratchpad/p345snap (staged index). Read-only review; mutations run in copies +(scratchpad/mut, scratchpad/mut2), nothing under the repo edited except this file. + +## Findings + +- **[fragile] tests/testthat/test-dft_stac_composite.R:568-571 — the `"Median"` pin cannot see the defect its comment says it guards, and nothing pins the flow at any of the three callers.** + The comment reads "neither may the aggregation check, which keeps the caller's case", but + `ckey()` calls `stac_composite_cache_key()` directly and never passes through + `aggregation_check()` or `dft_stac_composite()`. Measured by mutation (tabulated with + `reporter = "silent"`, since the summary reporter's failure lines do not match a + `Failure|Error` grep): + - `aggregation_check()` reverted to `tolower(aggregation)` (the round-2 bug itself): + test-dft_stac_composite.R **0 failed / 103 passed**, so the pin stays green. Only the unit + test in test-dft_stac_cube.R ("aggregation_check refuses anything...", 2 expectations) + goes red. + - `aggregation <- tolower(aggregation)` inserted right after the check in + `dft_stac_composite()`, `dft_stac_cube()` and `dft_stac_fetch()` (the same key move, + one line later): composite **0 failed / 103**, cube **0 / 141**, fetch **0 / 154**. + All green. + + So the round-2 mechanism, a normalisation that leaks into the cache key, is guarded only + by the helper's return value. It is not guarded at any of the six flows the findings.md + row enumerates, and the row's "Pinned the `"Median"` key; 3 mutations turn red" overstates + what the pin covers. The shipped code is correct (see below). This is about the guard. + **Fix:** drive the exported function with a mixed-case value and assert the key it produces. + For example, mock `stac_composite_cache_key` (and `stac_cube_cache_key` / + `stac_cache_key`) with a wrapper that records `aggregation` and then stops. Assert that + `"Median"` arrives as given, or assert the key against the v0.19.2 value. The probe below + already does this and could be lifted into the tests. Also correct the comment. + +## The mechanism check: every existing key, v0.19.2 vs now (all computed) + +Computed with `scratchpad/keyprobe.R`, run under `pkgload::load_all()` of each tree +(`git archive v0.19.2`, and the snapshot). The probe calls the **exported** functions on the +packaged AOI. It wraps the real key function, records the key, then stops before the network. +Fetch additionally mocks `stac_items_paged` and `gdalcubes::stac_image_collection`. + +| call | v0.19.2 | now | | +|---|---|---|---| +| composite default | 48e93230456245ee | 48e93230456245ee | same | +| composite `"Median"` | 15e150a8c375b3a3 | 15e150a8c375b3a3 | same | +| composite `"MEAN"` / `first` / `last` | b26b.. / 7e6c.. / 169e.. | identical | same | +| composite months `c(7,6)` / `7L` / `c(6,7,8)` | 4893.. / 0d20.. / d99f.. | identical | same | +| composite mask `c(11,3)` / `c(3L,8L)` | 8ae3.. / c88b.. | identical | same | +| composite tile 5000 / 5003 (snapped) | 77fb.. / 77fb.. | identical | same | +| composite clip `TRUE` / `"TRUE"` | d59d.. / d59d.. | identical | same | +| composite bands nir,red,green / `"red"`; res 20; resampling near; cc 100 | 72c7.. / 71a8.. / fda1.. / c102.. / 8fc3.. | identical | same | +| cube default / `"Median"` / `"LAST"` / mean | ccae.. / 0ea4.. / 6fc5.. / 6ccd.. | identical | same | +| cube dt P1Y / months 6:9 / mask / tile 5003 / clip 1 / ndvi | d5aa.. / 1cab.. / c07b.. / a61e.. / ccae.. / cc47.. | identical | same | +| fetch default / `"First"` / median / `"MEDIAN"` / tile / dt P1M | 3115.. / 0e86.. / 3631.. / 888b.. / a243.. / cc49.. | identical | same | +| composite `"count"` | 5dacaac91ec8ecad (`composite_`) | 4ed9b46904b621ef (`count_`) | intended | +| composite `"Count"` | cdcb020e1afcbf67 (`composite_`) | 4ed9b46904b621ef (`count_`) | intended (case-folded, new family) | +| cube `"count"`, fetch `"none"` | keyed | refused | intended, in NEWS | + +- **No existing non-count key moved**, for `aggregation` in any case, `bands`, `months`, + `mask_values`, `dt`, `tile_size`, `clip`, `res`, `resampling` or `cloud_cover_max`. +- **Branch changes by parameter.** The branch changes no normalisation other than + `aggregation` (`git diff --stat v0.19.2 HEAD -- R/` plus the staged diff touch only the three + `stac` files). `clip`, `tile_size`, `months` and `mask_values` normalise exactly as in + v0.19.2. +- **The `"COUNT"` to `"count"` fold.** It reaches only the new `count_` family. `family`, + `dt_read = "P1D"` and the `count_` filename prefix all derive from `is_count`. +- **gdalcubes sees the same values as in v0.19.2.** A non-count value goes to `cube_view()` as + given, the same as before; gdalcubes lower-cases on its own side. A count sends `"first"`. +- **`stac_cube_assemble()`'s new `aggregation_check()`** is a pure check. Its return value is + discarded, so it moves nothing. + +## The code fix itself: correct + +- `is_count <- identical(tolower(aggregation), "count")` is safe: + `aggregation_check()` guarantees a length-1, non-NA character. +- Only the count family is rewritten. +- The new "matched without case" test turns red if either the `tolower` or the + `aggregation <- "count"` line is dropped: two builds, two files, and `seen` of length 2 or + `"COUNT"`. +- Offline suites on the snapshot: + - composite: 0 failed / 3 skipped / 103 passed + - cube: 0 / 3 / 141 + - fetch: 0 / 4 / 154 +- `devtools::document()` in a copy reproduces `man/` and `NAMESPACE` byte for byte. + +## Also checked, no finding + +- NEWS "no existing cache key moves": true, per the table above. +- NEWS "the check runs before any network call": fetch and cube check first. Composite + checks right after `dft_stac_config()`, which is local. +- NEWS "0.46 GiB peak": 489,013,248 B = 0.455 GiB. "9.6 min": 577.48 s. +- `planning/` and `data-raw/` are both in `.Rbuildignore`, so the review files and the + benchmark logs do not ship. diff --git a/planning/active/task_plan.md b/planning/active/task_plan.md index 84f6f82..25666bf 100644 --- a/planning/active/task_plan.md +++ b/planning/active/task_plan.md @@ -74,19 +74,19 @@ Revised after the plan review (see findings.md, "Plan review" and "Count semanti `inst/notes/gdalcubes-pc-gotchas.md`. ## Phase 3: Docs -- [ ] Description, `@return` and the Caching section carry the count exception +- [x] Description, `@return` and the Caching section carry the count exception (integer counts, not reflectance; `count_.tif`). -- [ ] `@param aggregation` in `dft_stac_composite()` lists the valid values and describes +- [x] `@param aggregation` in `dft_stac_composite()` lists the valid values and describes `"count"`: distinct clear days per pixel, per band, integer, no scale; the `cloud_cover_max` pre-filter applies, so the count is of scenes that pass it. `dft_stac_cube()` and `dft_stac_fetch()` list their valid values. -- [ ] `@param mask_values`: the default SCL classes include **snow**, so a spring or autumn +- [x] `@param mask_values`: the default SCL classes include **snow**, so a spring or autumn "clear" count excludes snow as well as cloud (ask 4). -- [ ] Add a count example to `@examples` (inside the existing `\dontrun{}`, since it needs +- [x] Add a count example to `@examples` (inside the existing `\dontrun{}`, since it needs the network). Run `devtools::document()` and `pkgdown::check_pkgdown()`. ## Phase 4: Live verification and scale -- [ ] Network test (`DRIFT_TEST_NETWORK`): on the packaged AOI, the count is an integer, +- [x] Network test (`DRIFT_TEST_NETWORK`): on the packaged AOI, the count is an integer, ≥ 0, and ≤ the number of distinct item dates in the window. It must not look like reflectance, so its maximum must be ≥ 1 and must be a whole number. - [ ] Reproduce the issue's check on a small square: the old call gave continuous values diff --git a/tests/testthat/test-dft_stac_composite.R b/tests/testthat/test-dft_stac_composite.R index 16c71ad..5295292 100644 --- a/tests/testthat/test-dft_stac_composite.R +++ b/tests/testthat/test-dft_stac_composite.R @@ -296,7 +296,7 @@ test_that("a built composite leaves only its COG in the cache, no sidecar", { expect_equal(terra::time(out[[1]])[1], as.Date("2023-06-01")) }) -test_that("dft_stac_composite refuses an aggregation gdalcubes would not honour (#92)", { +test_that("dft_stac_composite refuses an aggregation outside the set drift passes (#92)", { skip_if_not_installed("gdalcubes") testthat::local_mocked_bindings( stac_cube_items = function(...) stop("reached the network") @@ -565,10 +565,13 @@ test_that("a count with no clear day anywhere drops the year, never caches zeros }) test_that("the count keys apart from every composite, and composite keys are unchanged (#92)", { - # pinned before #92: the family argument must not move an existing key + # pinned against v0.19.2: the family argument must not move an existing key. + # (That callers hand the key function the caller's case is tested in + # test-dft_stac_cube.R, through the exported functions.) expect_identical(ckey(), "03ee8ecc66b832a8") + expect_identical(ckey(aggregation = "Median"), "150c8ca5003bbe78") count <- ckey(aggregation = "count", dt = "P1D", family = "count") - # 0.19.x wrote reflectance under this key for aggregation = "count" + # 0.18.0-0.19.2 wrote reflectance under this key for aggregation = "count" expect_false(identical(count, "08e0c5510e8ae297")) expect_false(identical(count, ckey(aggregation = "count"))) # the tag alone separates the families @@ -576,3 +579,69 @@ test_that("the count keys apart from every composite, and composite keys are unc expect_identical(count, ckey(aggregation = "count", dt = "P1D", family = "count")) }) + +test_that("a live count is whole clear days, at most the distinct dates, and not reflectance (#92)", { + skip_if(Sys.getenv("DRIFT_TEST_NETWORK") != "true", + "network test; set DRIFT_TEST_NETWORK=true to run") + skip_if_not_installed("gdalcubes") + # A 2 km square at the packaged AOI's centroid, July 2023, cloudy scenes + # admitted so the mask has work to do. Measured 2026-09-30: 20 items on 12 + # dates (overlapping MGRS tiles), counts 5-9. + ctr <- sf::st_centroid(sf::st_union(sf::st_transform(aoi_pkg(), 32609))) + sq <- sf::st_as_sf(sf::st_buffer(ctr, 1000, endCapStyle = "SQUARE")) + cache <- withr::local_tempdir() + n <- suppressMessages(dft_stac_composite( + sq, years = 2023, months = 7, bands = "red", aggregation = "count", + res = 20, cloud_cover_max = 100, cache_dir = cache + ))[[1]] + v <- terra::values(n)[, 1] + # the distinct acquisition dates the same query returns + items <- rstac::stac("https://planetarycomputer.microsoft.com/api/stac/v1") |> + rstac::stac_search( + collections = "sentinel-2-l2a", + intersects = sf::st_geometry(sf::st_transform(sq, 4326))[[1]], + datetime = "2023-07-01T00:00:00Z/2023-07-31T23:59:59Z", limit = 500 + ) |> + rstac::ext_filter(`eo:cloud_cover` <= 100) |> + rstac::post_request() |> + rstac::items_fetch() + dates <- unique(substr(vapply(items$features, + function(f) f$properties$datetime, ""), 1, 10)) + expect_gt(length(dates), 1) + expect_true(all(v == round(v), na.rm = TRUE)) + expect_gte(min(v, na.rm = TRUE), 1) + expect_lte(max(v, na.rm = TRUE), length(dates)) + # the #92 value was reflectance, median ~0.03-0.04; a count is at least 1 + expect_gte(stats::median(v, na.rm = TRUE), 1) + f <- list.files(drift:::cache_scheme_dir(cache, "sentinel-2-l2a"), + full.names = TRUE) + expect_match(basename(f), "^count_[0-9a-f]{16}\\.tif$") + expect_equal(terra::datatype(terra::rast(f)), "INT2U") +}) + +test_that("a count is matched without case and keys as 'count' (#92)", { + skip_if_not_installed("gdalcubes") + aoi_t <- sf::st_transform(aoi_pkg(), 32609) + seen <- NULL + testthat::local_mocked_bindings( + stac_cube_items = function(...) list(features = list(), is_pre = FALSE), + stac_cube_assemble = function(fetched, cfg, aoi_target, target_crs, t0, t1, + res, dt, aggregation, ...) { + seen <<- c(seen, aggregation) + r <- terra::rast(terra::ext(aoi_t), resolution = 50, crs = "EPSG:32609", + vals = 3) + names(r) <- "red" + r + } + ) + cache <- withr::local_tempdir() + for (a in c("COUNT", "count")) { + suppressMessages(dft_stac_composite(aoi_pkg(), years = 2021, bands = "red", + aggregation = a, cache_dir = cache)) + } + # the upper-case call took the count path; the lower-case one hit its cache + expect_equal(seen, "first") + files <- list.files(drift:::cache_scheme_dir(cache, "sentinel-2-l2a")) + expect_length(files, 1L) + expect_match(files, "^count_") +}) diff --git a/tests/testthat/test-dft_stac_cube.R b/tests/testthat/test-dft_stac_cube.R index 37a65d2..5f96d7f 100644 --- a/tests/testthat/test-dft_stac_cube.R +++ b/tests/testthat/test-dft_stac_cube.R @@ -759,9 +759,10 @@ test_that("aggregation_check refuses anything outside the set drift passes to cu for (ok in c("min", "max", "mean", "median", "first", "last")) { expect_identical(drift:::aggregation_check(ok), ok) } - # case-insensitive, as gdalcubes is, and returned lower-cased for the cache key - expect_identical(drift:::aggregation_check("Median"), "median") - expect_identical(drift:::aggregation_check("LAST"), "last") + # case-insensitive, as gdalcubes is, and returned as given: callers hash it + # into cache keys, so lower-casing it would move a mixed-case caller's key + expect_identical(drift:::aggregation_check("Median"), "Median") + expect_identical(drift:::aggregation_check("LAST"), "LAST") bad <- list("count", "sum", "none", "count_values", "count_images", "", NA_character_, NA, c("median", "mean"), character(0), 1, NULL) for (b in bad) { @@ -814,3 +815,47 @@ test_that("dft_stac_cube and dft_stac_fetch refuse a bad aggregation before any # nothing was written for a refused call expect_length(list.files(cache, recursive = TRUE), 0L) }) + +test_that("fetch, cube and composite hash the caller's aggregation as given (#92)", { + # v0.19.2 hashed the value as passed, so "Median" and "median" had different + # keys. Lowering it anywhere between the check and the key moves every + # mixed-case caller's key and silently re-streams a cached cube. Drive the + # exported functions and record what reaches each key function. + skip_if_not_installed("gdalcubes") + aoi <- sf::st_read( + system.file("extdata", "example_aoi.gpkg", package = "drift"), + quiet = TRUE + ) + cache <- withr::local_tempdir() + hashed <- list() + record <- function(name) { + function(aoi_target, res, target_crs, a, b, ...) { + # the argument order differs: fetch/cube pass dt then aggregation, the + # composite passes datetime, dt, then aggregation + args <- list(...) + hashed[[name]] <<- if (name == "composite") args[[1]] else b + stop("recorded") + } + } + testthat::local_mocked_bindings( + stac_cache_key = record("fetch"), + stac_cube_cache_key = record("cube"), + stac_composite_cache_key = record("composite"), + # fetch queries and builds its collection before it keys + stac_items_paged = function(...) list(features = list(list(id = "x"))), + stac_cube_items = function(...) stop("reached the network") + ) + testthat::local_mocked_bindings( + stac_image_collection = function(...) NULL, .package = "gdalcubes" + ) + expect_error(suppressMessages( + dft_stac_fetch(aoi, source = "io-lulc", years = 2017, + aggregation = "Median", cache_dir = cache) + ), "recorded") + expect_error(dft_stac_cube(aoi, aggregation = "Median", cache_dir = cache), + "recorded") + expect_error(dft_stac_composite(aoi, years = 2021, aggregation = "Median", + cache_dir = cache), "recorded") + expect_identical(hashed, list(fetch = "Median", cube = "Median", + composite = "Median")) +}) From f3ef444d012788959a6e983c418c978cba3d9427 Mon Sep 17 00:00:00 2001 From: almac2022 Date: Wed, 30 Sep 2026 10:16:43 -0700 Subject: [PATCH 05/10] Benchmark clear-day count chips on BULK (#92) New chips-count stage of data-raw/benchmark_composite_bulk.R: the first 20 of #79's seed-79 chips (300 m buffers, 2023 Jul-Aug) as clear-day counts, bands = "red". 577.5 s real, median 24.7 s a chip, 0.46 GiB peak RSS (/usr/bin/time -l, committed as time_chips_count.txt); every chip's highest count 5-9, no NA cells. The same 20 chips as median true-colour composites took 1,051.6 s (#79). Relates to #92 Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01PRhUJsuKABLfpBGktPoiBN --- data-raw/benchmark_composite_bulk.R | 41 ++++++++++++++++++- .../chips_count_per_chip.csv | 21 ++++++++++ .../time_chips_count.txt | 4 ++ .../timings_chips-count.csv | 2 + planning/active/task_plan.md | 4 +- 5 files changed, 69 insertions(+), 3 deletions(-) create mode 100644 data-raw/logs/benchmark_composite_bulk/chips_count_per_chip.csv create mode 100644 data-raw/logs/benchmark_composite_bulk/time_chips_count.txt create mode 100644 data-raw/logs/benchmark_composite_bulk/timings_chips-count.csv diff --git a/data-raw/benchmark_composite_bulk.R b/data-raw/benchmark_composite_bulk.R index 1109581..a79fd75 100644 --- a/data-raw/benchmark_composite_bulk.R +++ b/data-raw/benchmark_composite_bulk.R @@ -9,9 +9,14 @@ # composited on its own: the floodplains#93 review path. # floodplain one floodplain-wide true-colour composite, tiled: the COG path # for titiler browsing. +# chips-count the first 20 of the same chips, as clear-day counts +# (aggregation = "count", bands = "red", #92). A count reads one +# time step per day rather than one per window, so it is the path +# whose memory and time could differ from a median. # # Rscript data-raw/benchmark_composite_bulk.R chips # Rscript data-raw/benchmark_composite_bulk.R floodplain +# /usr/bin/time -l Rscript data-raw/benchmark_composite_bulk.R chips-count # # Run through data-raw/benchmark_composite_bulk-run.sh, which samples RSS every # 2 s from outside the process and records the wall clock. @@ -20,7 +25,7 @@ pkgload::load_all(quiet = TRUE) stage <- commandArgs(trailingOnly = TRUE)[1] -stopifnot(stage %in% c("chips", "floodplain")) +stopifnot(stage %in% c("chips", "floodplain", "chips-count")) out_dir <- file.path("data-raw", "logs", "benchmark_composite_bulk") dir.create(out_dir, recursive = TRUE, showWarnings = FALSE) @@ -74,6 +79,40 @@ if (stage == "chips") { as.integer(stats::median(cells)), cache_mb)) } +if (stage == "chips-count") { + # the same draw as `chips`, so chip i here is chip i there + set.seed(79) + pts <- sf::st_as_sf(sf::st_sample(sf::st_union(aoi), 100)) + buf <- sf::st_buffer(pts, 300)[1:20, ] + per_chip <- numeric(nrow(buf)) + cells <- numeric(nrow(buf)) + count_max <- numeric(nrow(buf)) + count_na <- numeric(nrow(buf)) + tick("wall", { + for (i in seq_len(nrow(buf))) { + t0 <- Sys.time() + chip <- suppressMessages( + dft_stac_composite(buf[i, ], years = 2023, months = 7:8, bands = "red", + aggregation = "count", cache_dir = cache) + ) + per_chip[i] <- as.numeric(difftime(Sys.time(), t0, units = "secs")) + v <- terra::values(chip[[1]])[, 1] + cells[i] <- length(v) + count_max[i] <- max(v, na.rm = TRUE) + count_na[i] <- sum(is.na(v)) + } + }) + utils::write.csv(data.frame(chip = seq_along(per_chip), secs = per_chip, + cells = cells, count_max = count_max, + count_na = count_na), + file.path(out_dir, "chips_count_per_chip.csv"), + row.names = FALSE) + message(sprintf("count chips: %d, median %.1f s, max %.1f s, cells/chip %d, count max %d-%d, NA cells %d", + length(per_chip), stats::median(per_chip), max(per_chip), + as.integer(stats::median(cells)), as.integer(min(count_max)), + as.integer(max(count_max)), as.integer(sum(count_na)))) +} + if (stage == "floodplain") { # On failure, inventory R's temp files BEFORE the session deletes them: run 4 # (2026-09-28) read all 30 tiles, then terra::mask() could not read the merged diff --git a/data-raw/logs/benchmark_composite_bulk/chips_count_per_chip.csv b/data-raw/logs/benchmark_composite_bulk/chips_count_per_chip.csv new file mode 100644 index 0000000..0158ebd --- /dev/null +++ b/data-raw/logs/benchmark_composite_bulk/chips_count_per_chip.csv @@ -0,0 +1,21 @@ +"chip","secs","cells","count_max","count_na" +1,20.7404742240906,3660,7,0 +2,24.7813889980316,3660,9,0 +3,24.0691511631012,3660,9,0 +4,24.6106159687042,3660,8,0 +5,19.6276619434357,3660,7,0 +6,26.125659942627,3660,5,0 +7,26.3550560474396,3660,9,0 +8,58.7533140182495,3660,9,0 +9,25.2467560768127,3660,8,0 +10,26.88019490242,3660,5,0 +11,23.7180621623993,3660,9,0 +12,23.7300660610199,3660,8,0 +13,24.9521441459656,3660,6,0 +14,59.606479883194,3660,9,0 +15,37.0566339492798,3660,9,0 +16,36.0856280326843,3660,9,0 +17,22.4143421649933,3660,8,0 +18,23.1965670585632,3660,9,0 +19,23.1127099990845,3660,8,0 +20,23.9492878913879,3660,9,0 diff --git a/data-raw/logs/benchmark_composite_bulk/time_chips_count.txt b/data-raw/logs/benchmark_composite_bulk/time_chips_count.txt new file mode 100644 index 0000000..899d04a --- /dev/null +++ b/data-raw/logs/benchmark_composite_bulk/time_chips_count.txt @@ -0,0 +1,4 @@ +# /usr/bin/time -l Rscript data-raw/benchmark_composite_bulk.R chips-count (2026-09-30, macOS, 64 GB) + 577.48 real 44.54 user 8.04 sys + 489013248 maximum resident set size + 319064960 peak memory footprint diff --git a/data-raw/logs/benchmark_composite_bulk/timings_chips-count.csv b/data-raw/logs/benchmark_composite_bulk/timings_chips-count.csv new file mode 100644 index 0000000..64bfc0f --- /dev/null +++ b/data-raw/logs/benchmark_composite_bulk/timings_chips-count.csv @@ -0,0 +1,2 @@ +"stage","what","secs" +"chips-count","wall",575.057301998138 diff --git a/planning/active/task_plan.md b/planning/active/task_plan.md index 25666bf..ac92313 100644 --- a/planning/active/task_plan.md +++ b/planning/active/task_plan.md @@ -89,10 +89,10 @@ Revised after the plan review (see findings.md, "Plan review" and "Count semanti - [x] Network test (`DRIFT_TEST_NETWORK`): on the packaged AOI, the count is an integer, ≥ 0, and ≤ the number of distinct item dates in the window. It must not look like reflectance, so its maximum must be ≥ 1 and must be a whole number. -- [ ] Reproduce the issue's check on a small square: the old call gave continuous values +- [x] Reproduce the issue's check on a small square: the old call gave continuous values around 0.03, and the new max must be ≤ the distinct dates (6 in the issue's July 2021 square). -- [ ] Scale: count chips on BULK (about 20 points, 300 m buffers, one month) under +- [x] Scale: count chips on BULK (about 20 points, 300 m buffers, one month) under `/usr/bin/time -l`, recording peak RSS and wall clock, beside #79's 100-chip median run (0.48 GiB, 84.1 min). The floodplain-wide read is out of scope while #88 is open. Numbers go in the PR body and the archive README. From 72d78fb74f86cc97881f9d0fb03ac0f4fba8156d Mon Sep 17 00:00:00 2001 From: almac2022 Date: Wed, 30 Sep 2026 10:16:43 -0700 Subject: [PATCH 06/10] NEWS for #92 Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01PRhUJsuKABLfpBGktPoiBN --- NEWS.md | 14 ++++++++++++++ planning/active/task_plan.md | 8 ++++---- 2 files changed, 18 insertions(+), 4 deletions(-) diff --git a/NEWS.md b/NEWS.md index a5eada9..027bfb4 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,17 @@ +# drift 0.20.0 + +- **`dft_stac_composite(aggregation = "count")` returned reflectance, not a count, and now returns clear-day counts (#92).** + - **The cause:** `aggregation` went straight into `gdalcubes::cube_view()`. That function reads any string it does not know as no aggregation at all, with no error, and then copies every scene in, so a masked scene could overwrite a clear one. What came back looked like July red reflectance, 0.008-0.18. + - **Refused up front:** `dft_stac_composite()`, `dft_stac_cube()` and `dft_stac_fetch()` now refuse any `aggregation` other than `median`, `mean`, `min`, `max`, `first` and `last`. Each of the six was measured to survive a round trip through `cube_view()`; the set was not taken from `?cube_view`, which omits `last`. gdalcubes also honours `count_values` and `count_images`. drift refuses both, because they count items rather than days and every caller would read the result as reflectance, an index or a class code. Case is ignored, as gdalcubes ignores it, and the value is hashed as given, so no existing cache key moves. The check runs before any network call. `resampling` has the same silent fallback, to `near`, and is #96. +- **Clear-day counts.** `aggregation = "count"` is new in `dft_stac_composite()`. It returns the number of **distinct clear days** per pixel in each window, the number needed to choose composite windows. + - **What counts as one day:** each day is read as one time step, so overlapping MGRS tiles count one acquisition once. Live, on a 2 km square in July 2023, 20 items on 12 dates gave counts of 5-9, and 6 items on 4 clear dates gave 4 everywhere. + - **Format:** counts are whole numbers with no scale or offset, stored as an INT2U COG, `count_.tif`. + - **No clear day is `NA`, never 0.** gdalcubes returns such a pixel as 0 or `NA` depending on its chunk layout, which follows `parallel`, and a chunk that failed to read is `NA` too. A year with no clear day anywhere is dropped with a warning, like an empty composite. + - **The offset change:** a count spans the Sentinel-2 2022-01-25 offset change without refusal, since the offset does not change which pixels are masked. + - **"Clear" excludes snow too:** it means outside `mask_values`, whose Sentinel-2 default includes the SCL **snow** class as well as cloud, shadow and cirrus. A spring or autumn count therefore excludes snow cover. The docs now say so. +- **Cached "count" files from 0.18.0-0.19.2 are never read.** Since `dft_stac_composite()` arrived in 0.18.0, it cached what it returned for `aggregation = "count"`, reflectance, as `composite_.tif`. A count now keys under its own tag and file prefix, so nothing reads those files again, and every genuine composite key is unchanged (pinned in the tests). They sit in the current `v2` scheme and cannot be told apart from real composites, so `dft_cache_clear(scheme = "superseded")` does not reclaim them. To reclaim the space, delete `composite_*.tif` under `/v2/sentinel-2-l2a/`. That removes genuine composites too, which rebuild on their next call, but it leaves the `cube_*.tif` index caches alone. `dft_cache_clear(scheme = "current")` would delete those as well. +- **Floodplain scale.** 20 count chips on BULK (300 m buffers, 2023 Jul-Aug, `bands = "red"`): 9.6 min in total, median 24.7 s a chip, 0.46 GiB peak (`/usr/bin/time -l`). Each chip's highest count was 5-9 clear days, with no `NA` cells. The same 20 chips as #79's median true-colour composites took 17.5 min, median 44.4 s a chip. #79 recorded a 0.48 GiB peak over its full 100 chips, sampled every 2 s, so that figure reads low against the kernel peak above. The count reads daily steps but one band where the median reads three, which is likely why it is the cheaper of the two. The runs were two days apart, so this is not a controlled comparison. Evidence is in `data-raw/logs/benchmark_composite_bulk/*count*`. + # drift 0.19.2 - **`dft_rast_classify()` now classifies a raster that is already a factor (#91).** Before, it returned zero levels and no colour table, with no error or warning. This affected the published floodplain rasters, which carry their own attribute table, and a raster re-classified after an earlier `dft_rast_classify()`. The cause was `terra::unique()`, which returns a factor's labels rather than its codes, so no code in `class_table` matched. A factor's levels are now stripped on a copy before its codes are read. The output takes its levels and colours from `class_table`, and the caller's raster keeps its own. A `remap =` that matched a class already worked, because `terra::classify()` reads raw codes. A `remap =` that matched nothing had the bug, and now works. diff --git a/planning/active/task_plan.md b/planning/active/task_plan.md index ac92313..e3840f6 100644 --- a/planning/active/task_plan.md +++ b/planning/active/task_plan.md @@ -98,7 +98,7 @@ Revised after the plan review (see findings.md, "Plan review" and "Count semanti Numbers go in the PR body and the archive README. ## Phase 5: Release notes -- [ ] NEWS entry covering: +- [x] NEWS entry covering: - the silent reflectance, and which callers now refuse what - the new count and its semantics - the orphaned 0.19.x `composite_*.tif` files written under `"count"`, which are never @@ -108,7 +108,7 @@ Revised after the plan review (see findings.md, "Plan review" and "Count semanti ## Validation -- [ ] Tests pass (`devtools::test()`), `lintr::lint_package()` clean -- [ ] `/code-check` clean on each commit -- [ ] PWF checkboxes match landed work +- [x] Tests pass (`devtools::test()`), `lintr::lint_package()` clean +- [x] `/code-check` clean on each commit +- [x] PWF checkboxes match landed work - [ ] `/planning-archive` on completion, then `/gh-pr-push` (PR body: `Relates to NewGraphEnvironment/sred-2025-2026#16`) From 80617663bed19fd911230f10c6d3cceed48aea2f Mon Sep 17 00:00:00 2001 From: almac2022 Date: Wed, 30 Sep 2026 10:17:32 -0700 Subject: [PATCH 07/10] Archive planning files for issue #92 Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01PRhUJsuKABLfpBGktPoiBN --- .../README.md | 59 +++++++++++++++++++ .../findings.md | 0 .../progress.md | 0 .../review-p1-round1.md | 0 .../review-p1-round2.md | 0 .../review-p1-round3.md | 0 .../review-p2-round1.md | 0 .../review-p2-round2.md | 0 .../review-p2-round3.md | 0 .../review-p345-round1.md | 0 .../review-p345-round2.md | 0 .../review-p345-round3.md | 0 .../task_plan.md | 2 +- 13 files changed, 60 insertions(+), 1 deletion(-) create mode 100644 planning/archive/2026-09-issue-92-aggregation-count/README.md rename planning/{active => archive/2026-09-issue-92-aggregation-count}/findings.md (100%) rename planning/{active => archive/2026-09-issue-92-aggregation-count}/progress.md (100%) rename planning/{active => archive/2026-09-issue-92-aggregation-count}/review-p1-round1.md (100%) rename planning/{active => archive/2026-09-issue-92-aggregation-count}/review-p1-round2.md (100%) rename planning/{active => archive/2026-09-issue-92-aggregation-count}/review-p1-round3.md (100%) rename planning/{active => archive/2026-09-issue-92-aggregation-count}/review-p2-round1.md (100%) rename planning/{active => archive/2026-09-issue-92-aggregation-count}/review-p2-round2.md (100%) rename planning/{active => archive/2026-09-issue-92-aggregation-count}/review-p2-round3.md (100%) rename planning/{active => archive/2026-09-issue-92-aggregation-count}/review-p345-round1.md (100%) rename planning/{active => archive/2026-09-issue-92-aggregation-count}/review-p345-round2.md (100%) rename planning/{active => archive/2026-09-issue-92-aggregation-count}/review-p345-round3.md (100%) rename planning/{active => archive/2026-09-issue-92-aggregation-count}/task_plan.md (99%) diff --git a/planning/archive/2026-09-issue-92-aggregation-count/README.md b/planning/archive/2026-09-issue-92-aggregation-count/README.md new file mode 100644 index 0000000..9328993 --- /dev/null +++ b/planning/archive/2026-09-issue-92-aggregation-count/README.md @@ -0,0 +1,59 @@ +# #92: `dft_stac_composite(aggregation = "count")` returned reflectance + +## Outcome + +`aggregation` passed straight into `gdalcubes::cube_view()`. That function reads any value it does not know as `"none"` and raises no error. So `"count"` came back as plausible surface reflectance. + +**What was fixed.** + +- `dft_stac_fetch()`, `dft_stac_cube()` and `dft_stac_composite()` now refuse any aggregation outside six values: `min`, `max`, `mean`, `median`, `first` and `last`. Each of the six was measured to survive a `cube_view()` round trip. +- The check is case-insensitive, but it hashes the value exactly as given, so no existing cache key moves. +- `dft_stac_composite(aggregation = "count")` now returns distinct clear days per pixel. It reads with `dt = "P1D"` and `"first"`, then applies `reduce_time("count()")`. It applies no scale, offset or offset split, and turns 0 into NA. +- The count has its own cache family, written as `count_.tif` in INT2U with NEAREST overviews. + +**What was learned.** + +- `?cube_view` omits `last`, and gdalcubes also honours `count_values` and `count_images`. drift refuses those two by policy, because they count items rather than days. +- A pixel with no clear day comes back as 0 or NaN depending on chunk layout, and chunk layout follows `parallel`. +- `resampling` has the same silent fallback, to `near`. That is filed as #96. +- **Mistakes along the way.** Two were mine. The first all-masked probe used `t0 == t1`, and `reduce_time` passes a single time step through unchanged, so that probe proved nothing and was superseded. And the Phase 1 fix lower-cased `aggregation`, which would have silently re-streamed cached cubes for mixed-case callers. A later review round caught that, and it was reverted before merge. + +## Measurement + +**Chunk dependence (128 x 128 fixture).** At 256 px chunks the fixture gave 12,544 zeros. At 64 px chunks it gave 4,352 zeros and 8,192 NaN. After mapping 0 to NA, both chunkings are identical. This is why a pixel with no clear day is NA, never 0. + +**Live, on a 2 km square in July** (`res = 20`, 10,000 cells): + +| year | `cloud_cover_max` | items | distinct dates | counts | +|---|---|---|---|---| +| 2023 | 20 | 6 | 4 | 4 everywhere | +| 2023 | 100 | 20 | 12 | 5-9 | +| 2021 | 100 | 12 | 12 | 4-9 | + +- No NA cells in any run. +- Peak RSS 0.30-0.36 GiB. +- Overlapping MGRS tiles de-duplicate on real data: the 6 items above count as 4 days. + +**BULK, 20 count chips** (300 m buffers, 2023 Jul-Aug, `bands = "red"`): + +| run | total | median per chip | peak | +|---|---|---|---| +| count | 577.5 s | 24.7 s | 0.46 GiB (kernel peak) | +| same 20 chips as median composites (#79) | 1,051.6 s | 44.4 s | not measured for these 20 | + +- Each count chip's maximum was 5-9, with no NA cells. +- The count reads one band where the median reads three. That is the likely reason it is cheaper, but the two runs were not controlled against each other. + +**Review cost and yield.** One plan review and nine `/code-check` rounds. Every round that found something produced a mutation-verified fix. Rounds that found defects inside the previous fix were ended by enumeration, not by a reviewer calling the work finished: + +- the count build's collaborator arguments; +- every key flow, checked against v0.19.2. + +## Evidence + +- `data-raw/logs/benchmark_composite_bulk/*count*`: the BULK count chips. `time_chips_count.txt` holds the `/usr/bin/time -l` block. +- `data-raw/benchmark_composite_bulk.R`, stage `chips-count`. +- `inst/notes/gdalcubes-pc-gotchas.md`, section "Silent fallbacks and counting clear observations (#92)". +- `review-*.md` in this directory: each round's findings. + +Closed by: PR (link added on open) diff --git a/planning/active/findings.md b/planning/archive/2026-09-issue-92-aggregation-count/findings.md similarity index 100% rename from planning/active/findings.md rename to planning/archive/2026-09-issue-92-aggregation-count/findings.md diff --git a/planning/active/progress.md b/planning/archive/2026-09-issue-92-aggregation-count/progress.md similarity index 100% rename from planning/active/progress.md rename to planning/archive/2026-09-issue-92-aggregation-count/progress.md diff --git a/planning/active/review-p1-round1.md b/planning/archive/2026-09-issue-92-aggregation-count/review-p1-round1.md similarity index 100% rename from planning/active/review-p1-round1.md rename to planning/archive/2026-09-issue-92-aggregation-count/review-p1-round1.md diff --git a/planning/active/review-p1-round2.md b/planning/archive/2026-09-issue-92-aggregation-count/review-p1-round2.md similarity index 100% rename from planning/active/review-p1-round2.md rename to planning/archive/2026-09-issue-92-aggregation-count/review-p1-round2.md diff --git a/planning/active/review-p1-round3.md b/planning/archive/2026-09-issue-92-aggregation-count/review-p1-round3.md similarity index 100% rename from planning/active/review-p1-round3.md rename to planning/archive/2026-09-issue-92-aggregation-count/review-p1-round3.md diff --git a/planning/active/review-p2-round1.md b/planning/archive/2026-09-issue-92-aggregation-count/review-p2-round1.md similarity index 100% rename from planning/active/review-p2-round1.md rename to planning/archive/2026-09-issue-92-aggregation-count/review-p2-round1.md diff --git a/planning/active/review-p2-round2.md b/planning/archive/2026-09-issue-92-aggregation-count/review-p2-round2.md similarity index 100% rename from planning/active/review-p2-round2.md rename to planning/archive/2026-09-issue-92-aggregation-count/review-p2-round2.md diff --git a/planning/active/review-p2-round3.md b/planning/archive/2026-09-issue-92-aggregation-count/review-p2-round3.md similarity index 100% rename from planning/active/review-p2-round3.md rename to planning/archive/2026-09-issue-92-aggregation-count/review-p2-round3.md diff --git a/planning/active/review-p345-round1.md b/planning/archive/2026-09-issue-92-aggregation-count/review-p345-round1.md similarity index 100% rename from planning/active/review-p345-round1.md rename to planning/archive/2026-09-issue-92-aggregation-count/review-p345-round1.md diff --git a/planning/active/review-p345-round2.md b/planning/archive/2026-09-issue-92-aggregation-count/review-p345-round2.md similarity index 100% rename from planning/active/review-p345-round2.md rename to planning/archive/2026-09-issue-92-aggregation-count/review-p345-round2.md diff --git a/planning/active/review-p345-round3.md b/planning/archive/2026-09-issue-92-aggregation-count/review-p345-round3.md similarity index 100% rename from planning/active/review-p345-round3.md rename to planning/archive/2026-09-issue-92-aggregation-count/review-p345-round3.md diff --git a/planning/active/task_plan.md b/planning/archive/2026-09-issue-92-aggregation-count/task_plan.md similarity index 99% rename from planning/active/task_plan.md rename to planning/archive/2026-09-issue-92-aggregation-count/task_plan.md index e3840f6..99bcda1 100644 --- a/planning/active/task_plan.md +++ b/planning/archive/2026-09-issue-92-aggregation-count/task_plan.md @@ -111,4 +111,4 @@ Revised after the plan review (see findings.md, "Plan review" and "Count semanti - [x] Tests pass (`devtools::test()`), `lintr::lint_package()` clean - [x] `/code-check` clean on each commit - [x] PWF checkboxes match landed work -- [ ] `/planning-archive` on completion, then `/gh-pr-push` (PR body: `Relates to NewGraphEnvironment/sred-2025-2026#16`) +- [x] `/planning-archive` on completion, then `/gh-pr-push` (PR body: `Relates to NewGraphEnvironment/sred-2025-2026#16`) From 6f56c25acc54982c39e93c352295c06d08084b05 Mon Sep 17 00:00:00 2001 From: almac2022 Date: Wed, 30 Sep 2026 10:17:45 -0700 Subject: [PATCH 08/10] CLAUDE.md: the clear-day count call (#92) Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01PRhUJsuKABLfpBGktPoiBN --- CLAUDE.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CLAUDE.md b/CLAUDE.md index 98154bf..baf2226 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -83,6 +83,8 @@ dft_map_interactive(classified, aoi = aoi, rgb = tc) # or rgb = # sample points: ONE CALL PER BUFFERED POINT (small, separately cached chips). # All points in one aoi spans the floodplain bbox under one cache key. chips <- lapply(seq_len(nrow(buf)), \(i) dft_stac_composite(buf[i, ], years = 2023)) +# clear days per pixel, to choose a window (#92): integers, NA = no clear day, snow masked too +n <- dft_stac_composite(aoi, years = 2017:2023, months = 7, bands = "red", aggregation = "count") ``` Accuracy and error-adjusted area (#81; Olofsson et al. 2014, Stehman 2014). Mapped area is From 2e2597f9af8dcf2f10e5afd1c76454e5e2cf58c9 Mon Sep 17 00:00:00 2001 From: almac2022 Date: Wed, 30 Sep 2026 10:17:45 -0700 Subject: [PATCH 09/10] Release v0.20.0 (#92) Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01PRhUJsuKABLfpBGktPoiBN --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 95a5f44..adeff0b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: drift Title: Detecting Riparian and Inland Floodplain Transitions -Version: 0.19.2 +Version: 0.20.0 Date: 2026-09-29 Authors@R: c( person("Allan", "Irvine", , "al@newgraphenvironment.com", role = c("aut", "cre"), From 713979e44f933e7e73731fdcc585499c40410145 Mon Sep 17 00:00:00 2001 From: almac2022 Date: Wed, 30 Sep 2026 10:18:30 -0700 Subject: [PATCH 10/10] Archive README: link PR #97 Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01PRhUJsuKABLfpBGktPoiBN --- planning/archive/2026-09-issue-92-aggregation-count/README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/planning/archive/2026-09-issue-92-aggregation-count/README.md b/planning/archive/2026-09-issue-92-aggregation-count/README.md index 9328993..e5d36d8 100644 --- a/planning/archive/2026-09-issue-92-aggregation-count/README.md +++ b/planning/archive/2026-09-issue-92-aggregation-count/README.md @@ -56,4 +56,4 @@ - `inst/notes/gdalcubes-pc-gotchas.md`, section "Silent fallbacks and counting clear observations (#92)". - `review-*.md` in this directory: each round's findings. -Closed by: PR (link added on open) +Closed by: PR #97 (https://github.com/NewGraphEnvironment/drift/pull/97)