filter_geom evaluated and rejected on measurement; ship gdalcubes parallel instead (was: adopt the fixed filter_geom) #47

Outcome: filter_geom was adopted-ready and rejected on measurement

This issue's original premise is refuted. The body below has been rewritten to
say what is true, rather than leaving a superseded proposal reading as a spec.
The original argument is preserved in the history and summarised under What the
issue originally claimed.

The fork works. Adoption never depended on appelmar/gdalcubes#111 merging.
gdalcubes::filter_geom() is simply not worth using — measured, not assumed.

The measurement

data-raw/benchmark_filter_geom.R, five arms over the packaged Neexdzii AOI
(4-month monthly kNDVI, Sentinel-2 L2A), one subprocess per arm with
CPL_CURL_VERBOSE=YES and stderr to its own file, so cost is an observed count
of HTTP range requests rather than an inference from chunk arithmetic. Raw
evidence: data-raw/logs/benchmark_filter_geom/summary.csv, equivalence.csv.

Why the ~10× was never reachable

The issue reasoned from AOI/bbox ≈ 0.102 — "~90% of the streamed bounding box is
outside the polygon". That is the wrong bound. filter_geom skips whole
chunks, not pixels: filter_geom_cube::read_chunk() returns an empty chunk
without calling _in_cube->read_chunk(id). And gdalcubes:::.default_chunk_size()
targets 2 * parallel spatial chunks with the edge clamped to [64, 1024] px.

At parallel = 1 that is a 2 × 2 grid on a 3.3 km reach — a corridor
intersects all four, so the skip is zero. Predicted 462 requests; observed
462. Exactly neutral.

Forcing chunking to the 64 px floor does skip 26.7% of the ground, and still
loses: Sentinel-2 L2A COGs are Block=512x512, so a 64 px chunk sits inside
one source block and the same bytes are refetched per chunk. Chunk arithmetic
predicted 339 requests. The wire said 693.

A structural property is not a performance measurement.

The trap that made it worse than neutral

terra::mask() defaults to touches = TRUE; filter_geom rasterizes at
cell centre. So the issue's central proposal — "drop the now-redundant
output-side terra::mask()" — is not a redundancy removal. It silently
shrinks the analysed footprint:

mask() default   : 150 cells      <- touches = TRUE
mask(touches=F)  : 122 cells
true polygon area: 123.4 cells

Measured end-to-end: 49,244 → 41,608 non-NA cells, −15.5%. Values agree
exactly where both have data (correlation 1.000, max abs diff 0) — it is the
footprint that moves. floodplains judges parity at ±1 ha on 943 ha (0.1%); this
is 30–150× that tolerance, and it would have shipped looking like an optimisation.

Any future adoption must buffer the polygon by >= res·√2/2 and keep the
terra::mask(), so filter_geom is a read-bound only and the footprint cannot move.

Two more things that were not true

• Remotes: would not have pinned anything. It is a remotes/pak field,
  ignored by install.packages(), and — with gdalcubes in Suggests: and no
  version constraint — satisfied by any already-installed copy. It reads as a
  guarantee and is a hint.
• A capability probe was designed and then not needed. Worth recording that a
  version check could never have worked: the fork and the broken CRAN build both
  report 0.7.4 and both export filter_geom. Presence is not provenance. Only a
  behavioural probe discriminates, and it would have had to run in a subprocess,
  because the failure is bimodal — a hard segfault tryCatch cannot catch, or a
  silent all-NA cube.

What shipped instead

The benchmark surfaced a much larger, much cheaper win, and two live defects.

• dft_stac_cube() gains parallel, ~2× faster by default. drift never
  called gdalcubes::gdalcubes_options() anywhere in R/, so every cube read has
  run single-threaded since the function existed — not a considered choice, just
  the consequence of never setting it. Output is byte-identical at
  every setting, so it stays out of the cache key and no cached cube is orphaned.

  Corrected after review: an earlier version of this body said the old default
  "cost twice over", because gdalcubes derives its chunk size from parallel.
  The branch's own arms refute that — B_fg_default and C_fg_128 differ only
  in chunk size, both at parallel = 1, and measured 236.9 s / 462 requests
  against 343.7 s / 693. Finer chunking is a cost here. The speedup is
  concurrency alone.

• The end-to-end cube test passed on an all-NA cube. Its coverage assertion
  was mean(rowSums(!is.na(values)) == 0) > 0.5, which an empty cube satisfies
  with 1.0; the layer-count, time and cache-file assertions passed too. Every
  assertion in drift's only network cube test was satisfied by a cube containing
  no data — the exact failure mode dft_stac_cube: restore AOI-polygon clip (filter_geom segfault workaround) #32 exists to prevent. Replaced, and
  dft_stac_cube() now aborts before caching an empty cube.

• stac_cube_clip() documented the wrong clip rule (cell-centre, actually
  touches = TRUE). Corrected and pinned by a test using an irregular polygon —
  the existing test could not catch it, being axis-aligned on a cell boundary
  where both rules agree.

• tile_size is documented with its real cost. At 640 m it is the slowest
  option measured (1263.6 s / 3213 requests), because every tile rebuilds the
  image collection and reopens the COGs. It also lands sub-pixel-offset
  (correlation 0.996, recall 95.8%). Prefer parallel; reach for tile_size
  only when peak memory rather than wall clock binds.

What the issue originally claimed

That NewGraphEnvironment/gdalcubes@newgraph unblocked #38, that calling
filter_geom() would push the AOI into the read for roughly a 10× saving on the
0.102 area ratio, that the output-side terra::mask() could then be dropped as
redundant, and that the posture question was capability-probe vs Remotes: fork
pin. The fork claim was correct; the rest did not survive measurement.

Related

• Bound dft_stac_cube() streaming to the AOI (fetch-time read still bbox-bound after #32) #38 — parent, closed having named filter_geom as the way out. It is not.
• dft_stac_cube: restore AOI-polygon clip (filter_geom segfault workaround) #32 — the terra::mask() clip, which stays, and whose touches = TRUE
  semantics are now documented correctly
• Tile dft_stac_fetch to bound download over sparse-floodplain bounding boxes #36 / Integrate drift#36 tiled STAC fetch (tile_size) + benchmark speedup, parity, accuracy floodplains#8 — client-side tiling, independently
  measured and rejected there, and now measured here as the slowest arm
• inst/notes/gdalcubes-pc-gotchas.md — the durable record, so this is not
  re-litigated from the area ratio again
• Cache key moved on an rlang upgrade, silently orphaning every cached cube and fetch #48 — found while running this benchmark: both frozen cache-key goldens fail on
  main, meaning the key moved and silently orphaned every cached cube and
  fetch. Unrelated to this issue, but it is why the suite is FAIL 2 either side.
• Upstream appelmar/gdalcubes#110 / Add R CMD check workflow #111 — the segfault and its fix. Worth
  reporting the .default_chunk_size() × filter_geom interaction there: at
  parallel = 1 the default chunking makes filter_geom a guaranteed no-op.

Resolved by PR #49

Merged with filter_geom not adopted, dft_stac_cube(parallel =) shipped in
its place (drift 0.9.0), and two independent defects fixed — the end-to-end cube
test that passed on an all-NA cube, and stac_cube_clip()'s wrong clip rule.
Three adversarial review rounds found seven further defects in that work, four of
them introduced by it; all are recorded in
planning/archive/2026-09-issue-47-filter-geom-rejected/.