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 documents 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.
Corrected during the fix (PR #97). ?cube_view is incomplete. Round-tripping each value through cube_view()$aggregation shows gdalcubes also honours "last", "count_values" and "count_images", lower-cases its input, and maps anything else to "none". "count" is one of those unknown values, which is why it comes back as reflectance. The bug dates from 0.18.0, when dft_stac_composite() was introduced, not 0.19.0.
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 <- 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
- 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.
- Provide clear-observation counts per window. The obvious way is
dt = "P1D" in the view, then gdalcubes::reduce_time(cube, "count(<band>)") 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_<key>.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.
Resolution (PR #97, v0.20.0)
- Validation.
dft_stac_fetch(), dft_stac_cube() and dft_stac_composite() refuse anything outside min, max, mean, median, first and last, which is the measured set. The check ignores case and hashes the value exactly as given, so no existing cache key moves.
count_values and count_images are refused by policy. They count items, not days, and every caller would read the result as reflectance, an index or a class code.
- Counts: an
aggregation = "count" path, as proposed here.
- It reads with
dt = "P1D" and a day aggregation of "first", then applies reduce_time("count(<band>)"). There is no scale, no offset and no offset split.
- A pixel with no clear day is
NA, never 0. gdalcubes returns 0 or NaN for such a pixel depending on chunk layout, which follows parallel, and a failed read is NaN too.
- Cache. Counts key under their own tag and write
count_<key>.tif, so the 0.18.0-0.19.2 files are never read. Every existing composite key is unchanged.
- Acquisitions, not items. Same-day items collapse to one day. Live, 6 items on 4 dates counted 4. The docs now say that the default mask excludes snow.
resampling has the same silent fallback (to near), filed as #96.
Problem
dft_stac_composite(aggregation = "count")returns reflectance, not a count. It raises no error and no warning.aggregationgoes straight intogdalcubes::cube_view(aggregation = )(R/dft_stac_cube.R,stac_cube_assemble), which documents 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: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
aggregationagainst whatcube_viewsupports, indft_stac_composite()anddft_stac_cube(), and refuse anything else. A wrong value should never return plausible reflectance.dt = "P1D"in the view, thengdalcubes::reduce_time(cube, "count(<band>)")over the window, with no scale or offset applied to the result. Either a newaggregation = "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 everycomposite_<key>.tifwritten by 0.19.0 withaggregation = "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 defaultmask_values, so a spring or autumn "clear" count excludes snow as well as cloud. The docs should say so.Resolution (PR #97, v0.20.0)
dft_stac_fetch(),dft_stac_cube()anddft_stac_composite()refuse anything outsidemin,max,mean,median,firstandlast, which is the measured set. The check ignores case and hashes the value exactly as given, so no existing cache key moves.count_valuesandcount_imagesare refused by policy. They count items, not days, and every caller would read the result as reflectance, an index or a class code.aggregation = "count"path, as proposed here.dt = "P1D"and a day aggregation of"first", then appliesreduce_time("count(<band>)"). There is no scale, no offset and no offset split.NA, never 0. gdalcubes returns 0 or NaN for such a pixel depending on chunk layout, which followsparallel, and a failed read is NaN too.count_<key>.tif, so the 0.18.0-0.19.2 files are never read. Every existing composite key is unchanged.resamplinghas the same silent fallback (tonear), filed as #96.