diff --git a/.claude/sweep-performance-state.csv b/.claude/sweep-performance-state.csv index 3b09bd31c..77efd2a00 100644 --- a/.claude/sweep-performance-state.csv +++ b/.claude/sweep-performance-state.csv @@ -18,7 +18,7 @@ fire,2026-03-31T18:00:00Z,SAFE,compute-bound,0,, flood,2026-06-25,RISKY,compute-bound,1,3503,"Cat1 HIGH: _validate_mannings_n_dataarray used .values, eagerly materialized dask/cupy roughness raster (OOM); fixed lazy-safe. Core elementwise ops SAFE. LOW: dask+cupy host round-trip (map_blocks b.get) shared with cost_distance/surface_distance, documented not fixed." focal,2026-05-29,SAFE,compute-bound,1,2734,"HIGH: _hotspots_dask_cupy chunk fn round-tripped each chunk host<->GPU (cupy.asnumpy classify cupy.asarray); fixed PR 2739 to reuse _run_gpu_hotspots on device. LOW (not fixed): _apply_numpy/_hotspots_cupy use zeros_like where empty would suffice. CUDA kernels regs<=62, no register-pressure issue." geodesic,2026-03-31T18:00:00Z,N/A,compute-bound,0,, -geotiff,2026-07-05,SAFE,IO-bound,0,,"Pass 18 (2026-07-05): re-audit, 0 CRITICAL/0 HIGH/0 MEDIUM. Only 2 commits since Pass 17: #3604 (f0275b61, O(1) ndim gate in _writers/eager.py before dispatch -- no perf surface) and #3605 (6b0887fb, tests only). Category scans across the whole subpackage came back matching Pass 17's known-clean state: Cat1 .compute()/.get() sites are the documented streaming writes (_writer.py row bands, _writers/gpu.py per-band stream #3166/#3117) and batched D2H; Cat2 50k-task cap intact at _backends/dask.py:580, #3597 StreamingStats single-pass intact, _dask_finite_stats fuses 5 reductions into one dask.compute; Cat3 all Device().synchronize() sites batch-end/error-recovery (#2212/#2107 fixes intact); Cat5 _overview_kernels.py @ngjit row-major float-stable. Probes on-device this pass (PYTHONPATH=worktree, verified xrspatial.__file__): dask read 2560x2560 deflate chunks=256 = 400 tasks/100 chunks (4/chunk); eager numpy parity 0.0; eager GPU read returns cupy, parity 0.0, 400ms incl warmup, no host round trip; dask+cupy lazy (401 tasks), corner chunk is cupy, parity 0.0. _gpu_decode.py untouched since Pass 17 so recorded register values stand (inflate 67 / LZW 29; introspection still unavailable via numba _Kernel API). Prior deferred LOWs unchanged (redundant .copy() in _writer predictor encode, identity pack arithmetic, twice-built IFD byte length, per-level device sync in _block_reduce_2d_gpu, _nvcomp_batch_compress adler32 staging, _CloudSource.read_range per-call open). No HIGH/MEDIUM => no /rockout, no new benchmark required. SAFE/IO-bound holds. | Pass 17 (2026-07-02): re-audit, 0 CRITICAL/0 HIGH/0 MEDIUM. Prior row logged as needing GPU re-validation; CUDA IS available on this host and all GPU decode/read paths were re-run on-device this pass (no cuda-unavailable token). Probes: dask read 2560x2560 chunks=256 = 400 tasks / 100 chunks (4/chunk, 50k cap intact at _backends/dask.py:580); eager numpy parity 0.0; eager GPU read of deflate-tiled file returns cupy, parity 0.0, 165ms, no host round trip; dask+cupy read lazy (401 tasks) producing cupy corner parity 0.0; GPU deflate+fp-predictor(predictor=3) decode parity 0.0 (routes through nvcomp GPU path on this host, numba _inflate_tiles fallback not specialized). Cat1: no .values/np.asarray/np.array on dask or cupy arrays in _backends/; the _writer.py:1221-1394 .compute() calls are the intentional streaming row-band writes (single-pass verified in #3117/#3235/#3597, not per-chunk re-execution). Cat3: all Device().synchronize() in _gpu_decode.py (1008/1029/1264/1758/2531/2824) are batch-end or error-recovery, not per-tile in loops -- #2212/#2107 fixes intact. Cat5: _overview_kernels.py are @ngjit CPU kernels (row-major inner loop, float-stable accumulators), not CUDA. Reviewed the 6 commits since last inspection (#3595/#3599 PAM overwrite cleanup, #3592/#3598 + #3593/#3596 docstrings, #3594 non-finite RAT guard, #3600 color_ramp streaming accumulation = Pass 16's own fix, #3588/#3589 isort): all correctness/docs, no new performance surface, no regressions. Register introspection unavailable in this numba's _Kernel API; kernels unchanged since Pass 7 (inflate 67 regs / LZW 29 regs). Prior deferred LOWs unchanged (redundant .copy() in _writer predictor encode, identity pack arithmetic, twice-built IFD byte length). No HIGH/MEDIUM => no /rockout, no new benchmark required. SAFE/IO-bound holds. | Pass 16 (2026-07-01): 1 MEDIUM found and fixed, 0 HIGH. to_geotiff(dask, color_ramp=...) executed the source graph twice: streaming write computed every chunk, then _write_sidecars -> _finite_stats ran dask.compute over the same source for the PAM/QML statistics (measured 32 chunk executions for a 16-chunk source; color_ramp_range escape hatch stays at 16). Filed #3597, fixed by threading a chunk_observer through _write_streaming's three materialisation sites (row band, wide-raster segment, strip band) into a new StreamingStats accumulator (_symbology.py; Chan mean/M2 combine, float64 accumulators, ddof=0, nodata/finite exclusion matches _finite_stats; all-valid buffers skip the boolean-index copy and reductions use dtype=float64 on the original buffer to avoid an astype copy). Post-fix executions 16/16; round-trip bench 8192x8192 f32 deflate write+stats 1.53s -> 1.44s and total source reads halve (the real win is expensive-to-recompute sources: HTTP COGs, cold storage, long pipelines). GPU writer and VRT paths keep the documented second pass (GPU writer fully materialises anyway); docstrings updated. 13 new tests in test_color_ramp_single_pass_3597.py incl. execution counters (row-band, strip, segmented wide path with full-width chunks), nodata/int/all-NaN/multiband gating, accumulator-vs-_finite_stats parity, and a dask+cupy gpu=False streaming leg (run on-device). Audited all 26 geotiff commits since 2026-06-11: pack range guards #3272/#3277 correctly defer per-chunk scans into the write's single compute (no #3235 regression); #3374 fixed Pass 14's deferred LOW (chunked GPU read now parses header/IFDs once); xarray engine #3375/#3377/#3380 is a thin wrapper (no eager compute); PAM sidecar read on open_geotiff is a local os.path.exists probe (no HTTP cost); _CloudSource cat_file change is neutral-to-better. Probes on-device this pass: dask read 400 tasks/100 chunks (4/chunk, cap intact), eager GPU read cupy parity 0.0 (222ms), dask+GPU lazy with cupy meta, symbology _finite_stats GPU/CPU parity 0.0. LOWs unchanged from prior passes. SAFE/IO-bound holds. | Pass 15 (2026-06-11): 1 MEDIUM found and fixed. _pack (_attrs.py:~1795) guarded the no-sentinel integer restore with an eager bool(out.isnull().any()), which executed the whole upstream dask graph at to_geotiff(pack=True) call time; the streaming writer then executed it again, so every source chunk computed twice (measured 32 decode-task executions for 16 chunks on a 512x512 int16 SCALE/OFFSET no-GDAL_NODATA source; 71->33 total task starts post-fix). Filed #3235, fixed by mapping a per-chunk NaN guard (_pack_guard_no_nan) into the graph for dask-backed data (raises from the write's single compute; numpy keeps the eager call-time check; meta= preserves cupy backing). 9 new tests in test_pack_lazy_nan_guard_3235.py incl. fusion-proof execution counter and cupy-chunk guard unit test (dask+cupy e2e still blocked upstream by #3112). Scrutinised all 16 commits since 2026-06-08 (pack/unpack series #3065/#3075/#3079/#3129/#3174/#3175, VRT placement #3135, compression_level gate #3176, streaming banding #3136, dask+cupy writer order fix #3171): no other regressions; #3171's get-then-asarray order is intentional D2H for gpu=False. GPU validated on-device this pass: eager GPU unpack returns cupy with exact parity (387ms incl warmup, only 0-d scalar .get()s -- no bulk host round trip), dask+GPU unpack lazy (112 tasks/16 chunks, cupy meta, compute returns cupy, parity 0.0), GDS fast path intact without unpack (4 tasks/chunk); unpack disqualifying GDS is documented intentional. Dask CPU probe 4 tasks/chunk, 50k-task cap intact. Note: #1714 (_write_vrt_tiled synchronous scheduler) is now FIXED+CLOSED (scheduler='threads' at _writers/eager.py:1517) -- drop from the open-issue list. LOW noted (no fix): _pack does identity (data-0.0)/1.0 arithmetic allocating two full-array temporaries when scale==1/offset==0 (masked_nodata-only pack); prior deferred LOWs unchanged. SAFE/IO-bound holds. | Pass 14 (2026-06-09): MEDIUM found and fixed -- _write_streaming ran one dask .compute() per 256-row tile-row/strip, so a source chunk taller than the band re-executed once per band it overlapped (measured 2x at chunks=512, 4x at chunks=1024, whole upstream graph re-runs for computed pipelines). Filed #3117, fixed via _stream_row_bands: consecutive tile-rows/strips group into row bands sized by the source chunk-row span (one-chunk halo, #3007 accounting) under streaming_buffer_bytes; each band computes once and tiles/strips are carved from the materialised band. Wide rasters needing column segmentation keep the per-tile-row path. Post-fix per-chunk executions == 1 on the default read->write round trip. 5 new tests (TestRowBandRecompute3117 + _stream_row_bands unit); write/integration/parity suites pass (2195). LOW deferred (no fix): _read_geotiff_gpu_chunked parses header+all IFDs twice at graph build (_backends/gpu.py ~1367-1419, cap check then GDS probe; build-time only). GPU paths validated on-device this pass: eager gpu read returns cupy with parity, dask+GPU chunked read lazy (17 tasks/4 chunks) with parity; GPU writer full materialisation is documented intentional (streaming_buffer_bytes no-op). Read path keeps 50k-task graph cap; dask read probe 4 tasks/chunk. SAFE/IO-bound holds. | Pass 13 (2026-05-20): 1 MEDIUM found and fixed. _nvjpeg_batch_encode (_gpu_decode.py:~L1560) and _nvjpeg2k_batch_encode (~L2958) called cupy.cuda.Device().synchronize() inside the per-tile encode loops, a whole-device fence that blocked every CUDA stream and serialised concurrent work (e.g. predictor encodes on other streams). The decode-side counterpart _try_nvjpeg_batch_decode already used cupy.cuda.Stream.null.synchronize() at L1442; the encoder side was inconsistent. Filed #2212 and fixed both encoders to use Stream.null.synchronize(), scoping the per-tile sync to the default stream the encode/retrieve calls were issued on. nvJPEG / nvJPEG2000 encoders maintain a single shared state per encoder so encodes within a batch are inherently serial; the fix removes the device-wide blocker without changing the API ordering contract. 5 new tests in test_nvjpeg_encode_stream_sync_2212.py (AST checks that neither encoder contains Device().synchronize() inside a for-loop, that both call Stream.null.synchronize() in the loop, and that the decoder reference pattern stays pinned). All 5 new tests + 19 existing related encode/decode tests pass. nvjpeg/nvjpeg2k shared libs not present on this host so end-to-end encode verification is gated; add cuda-unavailable-libs note to re-validate on a host with the RAPIDS conda env. SAFE/IO-bound verdict holds; no change in dask graph cost. Dask probe: 2560x2560 deflate-tiled file via read_geotiff_dask(chunks=256) yields 400 tasks for 100 chunks (4 tasks/chunk), well under the 50K cap. LOW deferred (no fix in this PR): _build_ifd called twice per IFD level in _assemble_standard_layout (_writer.py:1531+1543), _assemble_cog_layout (1582+1625), and the COG overview path (2519+2546+2740) -- the first call's bytes are discarded; only the overflow byte length is used to compute pixel_data_offset. Cost is bounded by IFD count (typically 1-5 overview levels) so absolute impact is minor. Pre-existing pattern. | Pass 12 (2026-05-18): 1 MEDIUM found and fixed. _try_nvjpeg2k_batch_decode at _gpu_decode.py:~L2725-2778 allocated per-tile per-component cupy.empty buffers (N*S round-trips through the cupy memory pool) and called cupy.cuda.Device().synchronize() once per tile, forcing default-stream serialisation that defeats nvJPEG2000's internal pipelining. Filed #2107 and fixed: pre-allocate a single d_comp_pool sized n_tiles*samples*tile_height*pitch under a _check_gpu_memory guard, derive per-tile/per-component views as slab offsets, and replace the per-tile sync with a single batch-end sync. Same pattern as #1659 (_try_nvcomp_from_device_bufs), #1688 (_try_kvikio_read_tiles), #1712 (_nvcomp_batch_compress). 7 new tests in test_nvjpeg2k_single_alloc_2107.py: AST-level structural assertions confirm no cupy.empty inside the for-loop and no Device().synchronize() inside the loop, plus pool/per_tile_comp_bytes presence and _check_gpu_memory guard checks; lib-absent short-circuit; unsupported-dtype cleanup contract; cupy-only pool slab-non-overlap test (gpu-marked). libnvjpeg2k.so not present on this host so the end-to-end nvJPEG2000 decode is gated -- note added to re-validate on a host with the RAPIDS conda env. All 30 jpeg2000/compression tests + 7 new tests pass. SAFE/IO-bound verdict holds (no change in dask graph cost). Dask probe: 4096x4096 deflate-tiled file via read_geotiff_dask(chunks=512) yields 256 tasks for 64 chunks (4 tasks/chunk), well under the 50K cap. | Pass 11 (2026-05-18): 1 MEDIUM found and fixed. _read_strips (_reader.py:~L1972) and _fetch_decode_cog_http_strips (_reader.py:~L2670) decoded strips sequentially in a Python for-loop while the tile counterparts (_read_tiles L2146, _fetch_decode_cog_http_tiles L2898) gated parallel decode on _PARALLEL_DECODE_PIXEL_THRESHOLD via ThreadPoolExecutor. Filed #2100 and fixed: both strip paths now collect jobs, parallel-decode when n_strips > 1 and strip_pixels >= 64K, then place sequentially. Measured (uint16, 4-core): 4096x4096 deflate 130ms->34ms (3.82x), 8192x8192 deflate 531ms->146ms (3.63x), 8192x8192 zstd 211ms->85ms (2.48x), uncompressed 25ms->22ms (1.14x). 5 new tests in test_parallel_strip_decode_2100.py (parallel/serial parity, pool-engaged on multi-strip, serial-path for single-strip, windowed cross-strip read, HTTP COG strip parity). 3998 tests pass; 8 pre-existing failures predating this change (predictor2 BE + size_param_validation_gpu_vrt reference now-private read_to_array attr). SAFE/IO-bound verdict holds. | Pass 10 (2026-05-15): 1 new MEDIUM found and fixed; 2 LOW noted. MEDIUM (_reader.py:2737): _fetch_decode_cog_http_tiles decoded tiles sequentially in a Python for-loop after the concurrent fetch landed (issue #1480). Local _read_tiles parallelises decode whenever tile_pixels >= 64K via ThreadPoolExecutor (_reader.py:2017); the HTTP path was structurally similar but never picked up the same gate, so wide windowed reads of multi-tile COGs left deflate/zstd decode single-threaded. Mirrored the local-path threshold + pool. 5 new tests in test_cog_http_parallel_decode_2026_05_15.py (parallel + serial round-trip correctness, pool-instantiation branch selection above the threshold, single-tile path skips the pool, structural _decode_strip_or_tile call count == n_tiles). All 262 COG/HTTP tests pass; 3162 of 3164 selected geotiff tests pass overall (2 pre-existing failures predating Pass 9 per prior notes -- test_predictor2_big_endian_gpu_1517 references the now-private read_to_array attr, and the test_size_param_validation_gpu_vrt_1776 tile_size=4 validator failure). LOW deferred (no fix in this PR): (1) _block_reduce_2d_gpu (_gpu_decode.py:3142/3163/3189) does bool(mask.any().item()) per overview level when nodata is set, paying one device sync per level; the alternative (unconditional cupy.putmask) always pays the work cost and the short-circuit is correct under the current API. (2) _nvcomp_batch_compress adler32 staging (_gpu_decode.py:2543-2546) issues n_tiles slice-assign kernels into a fresh contig buffer despite all callers passing slices of a single underlying d_tile_buf; an API refactor to accept the source buffer directly would skip the rebuild. SAFE/IO-bound verdict holds. Dask probe: 2560x2560 chunks=256 yields 400 tasks (4 per chunk), well under the 50000 cap. GPU probe: 1024x1024 float32 zstd read returns CuPy-backed in 236 ms with no host round-trip. | Rockout 2026-05-15: LOW filed #1934 -- _apply_nodata_mask_gpu used cupy.where (allocating); switched to cupy.putmask on the already-owned buffer (float path) and on the post-astype float64 buffer (int path). Saves one chunk-sized device allocation per call. 7 new tests in test_apply_nodata_mask_gpu_inplace_1934.py; 52 related nodata tests pass. | Pass 8 (2026-05-12): 1 new MEDIUM found and fixed. _assemble_standard_layout/_assemble_cog_layout returned bytes(bytearray), doubling peak memory transiently during eager writes. Filed #1756, fixed by returning the bytearray directly. Measured: 95 MB uint8 raster peak drops 202 MB -> 107 MB. _write_bytes / parse_header already accepted the buffer protocol so the change is transparent to callers. 6 new tests in test_assemble_layout_no_bytes_copy_1756.py. 2123 existing geotiff tests pass; the 10 unrelated failures (test_no_georef_windowed_coords_1710, test_predictor2_big_endian_gpu_1517) reference the now-private read_to_array attribute (commit 8adb749, issue #1708) and predate this change. SAFE/IO-bound verdict holds. | Pass 7 (2026-05-12): re-audit identified 4 MEDIUM findings, all real, all backed by microbenches. (1) unpack_bits sub-byte loops for bps=2/4/12 in _compression.py:836-878 were 100-200x slower than vectorised numpy (filed #1713, fixed in this branch: bps=4 2M pixels drops from 165ms to 3ms = 55x; bps=2/12 similar). (2) _write_vrt_tiled at __init__.py:1708 uses scheduler='synchronous' on independent tile writes; measured 33% slowdown on 256-tile zstd write vs threads scheduler (filed #1714, no fix yet). (3) _nvcomp_batch_compress at _gpu_decode.py:2522-2526 still does per-tile cupy.get().tobytes() despite #1552 / #1659 fixing the same pattern elsewhere; measured 45% reduction with concat+single get on n=1024 (filed #1712, no fix yet). (4) _nvcomp_batch_compress at _gpu_decode.py:2457 uses per-tile cupy.empty allocations; 1024 tiles 16KB drops from 4.7ms to 1.0ms with single contiguous + views (bundled into #1712). Cat 6 OOM verdict: SAFE/IO-bound holds -- read_geotiff_dask caps task count at _MAX_DASK_CHUNKS=50_000 and per-chunk memory is bounded by chunk size. _inflate_tiles_kernel resource usage on Ampere: 67 regs/thread, 2896B local/thread, 8192B shared/block (LZW kernel: 29 regs, 24576B shared) -- register pressure under control; high local memory in inflate is unavoidable (LZ77 state) but only thread 0 in each block uses it. | Pass 4 (2026-05-10): re-audit after #1559 (centralise attrs across all read backends). New _populate_attrs_from_geo_info helper at __init__.py:301 runs once per read, not per-chunk -- no perf impact. Probe: 2560x2560 deflate-tiled file opened via read_geotiff_dask yields 400 tasks (4 tasks/chunk for 100 chunks), well under 1M cap. read_geotiff_gpu(1024x1024) returns cupy.ndarray end-to-end with no host round-trip (226ms incl. write+decode). No new HIGH/MEDIUM findings. SAFE/IO-bound holds. | Pass 3 (2026-05-10): SAFE/IO-bound. Audited 4 perf commits: #1558 (in-place NaN writes on uniquely-owned buffers correct), #1556 (fp-predictor ngjit ~297us/tile for 256x256 float32), #1552 (single cupy.concatenate + one .get() for batched D2H at _gpu_decode.py:870-913), #1551 (parallel decode threshold >=65536px engages 256x256 default at _reader.py:1121). Bench: 8192x8192 f32 deflate+pred2 256-tile write 782ms; 4096x4096 f32 deflate read 83ms with parallel decode. Deferred LOW (none filed, all <10% MEDIUM threshold): _writer.py:459/1109 redundant .copy() before predictor encode (~1% per tile), _compression.py:280 lzw_decompress dst[:n].copy() (~2% per LZW tile decode), _writer.py:1419 seg_np.copy() before in-place NaN substitution (negligible, conditional path), _CloudSource.read_range opens fresh fsspec handle per range (pre-existing, predates audit scope). nvCOMP per-tile D2H batching break-even confirmed (variable sizes need staging buffer, no win). | Pass 3 (2026-05-10): audited f157746,39322c3,f23ec8f,1aac3b7. All 5 commits correct. Redundant .copy() in _writer.py:459,1109 and _compression.py:280 (1-2% overhead, LOW). _CloudSource.read_range() per-call open is pre-existing arch issue. No HIGH/MEDIUM regressions. SAFE. | re-audit 2026-05-02: 6 commits since 2026-04-16 (predictor=3 CPU encode/decode, GPU predictor stride fix, validate_tile_layout, BigTIFF LONG8 offsets, AREA_OR_POINT VRT, per-tile alloc guard). 1M dask chunk cap intact at __init__.py:948; adler32 batch transfer intact at _gpu_decode.py:1825. New code is metadata validation and dispatcher logic with no extra materialization or per-tile sync points. No HIGH/MEDIUM regressions. | Pass 5 (2026-05-12): re-audit identified MEDIUM in _gpu_decode.py:1577 _try_nvcomp_from_device_bufs: per-tile cupy.empty + trailing cupy.concatenate doubled peak VRAM and added serial concat. Filed #1659 and fixed to single-buffer + pointer offsets (matches LZW/deflate/host-buffer patterns at L1847/L1878/L1114). Microbench (alloc+concat overhead only, not full nvCOMP latency): n=256 tile_bytes=65536 drops 3.66ms->0.69ms, n=256 tile_bytes=262144 drops 8.18ms->0.13ms. Tests: 5 new tests in test_nvcomp_from_device_bufs_single_alloc_1659.py (codec short-circuit, no-lib short-circuit, memory-guard contract, real ZSTD round-trip via nvCOMP, structural single-buffer check). 1458 existing geotiff tests pass, 3 unrelated matplotlib/py3.14 failures pre-existing. SAFE/IO-bound verdict holds. | Pass 6 (2026-05-12): re-audit on top of #1659. New HIGH in _try_kvikio_read_tiles at _gpu_decode.py:941: per-tile cupy.empty() + blocking IOFuture.get() inside loop serialised GDS reads to ~1 outstanding pread, missed parallelism the kvikio worker pool was designed for, paid per-tile cupy.empty setup (matches #1659 anti-pattern in nvCOMP path), and lacked _check_gpu_memory guard. Filed #1688 and fixed to single contiguous buffer + batched submit + guard. Microbench with 8-worker pool simulation: 256 tiles@1ms latency drops 256ms->38.7ms (~6.6x); single-thread simulation 256ms->28.5ms (9x). Tests: 9 new tests in test_kvikio_batched_pread_1688.py (kvikio-absent path, single-buffer pointer arithmetic, submit-before-get ordering, memory guard, partial-read fallback, round-trip data, zero-size/all-sparse tiles). All 1577 geotiff tests pass except pre-existing matplotlib/py3.14 failures." +geotiff,2026-07-27,SAFE,IO-bound,0,#3695,"Pass 19 (2026-07-27): 1 MEDIUM found and fixed, 0 CRITICAL/0 HIGH. Only ONE geotiff commit since Pass 18 (f1f81b8d #3643, _parse_rat pads category_names to pixel indices for sparse RATs) -- metadata-only, bounded by a new _MAX_CATEGORIES=1e6 fail-closed gate, runs once per open behind a local os.path.exists probe, no perf surface. Confirmed 5aa309e4 (datashader decoupling) touched ZERO geotiff files. MEDIUM (#3695 / PR #3696): to_geotiff(dask, cog=True, color_ramp=...) executed the source graph TWICE -- cog=True skips the streaming writer so the #3597 chunk_observer never fired and _write_sidecars fell back to _finite_stats -> _dask_finite_stats on the still-lazy _sym_data. Measured 33 chunk executions for a 16-chunk source (2.06x) vs 1.06x for every other cog/color_ramp combination; post-fix 17 (1.00x). Fixed by folding the already-materialised buffer into the same StreamingStats accumulator (snapshot taken BEFORE the sentinel restore so stats match _finite_stats semantics); sidecar bytes unchanged. 15 new tests in test_color_ramp_cog_single_pass_3695.py (7 fail without the fix) incl. dask+cupy and cupy legs run on-device; asv WriteCOGSymbology added (65.8ms -> 48.5ms on a cheap in-memory source, which is the FLOOR of the win -- expensive sources see ~2x). Full geotiff suite 6132 passed. This is the 4th recurrence of the 'write path re-executes the source graph' class (#3117/#3235/#3597/#3695); the streaming path was fixed each time and a sibling branch kept the second pass. REGISTER INTROSPECTION NOW RESOLVED (was 'unavailable' in passes 7-18): CUDADispatcher.compile(sig) + get_regs_per_thread() works on this numba. All 15 @cuda.jit kernels force-compiled and measured -- inflate 67 regs / 2896B local / 8192B shared / maxTPB 896; LZW 29 regs / 24576B shared; predictor decode u8-u64 23-25; predictor encode 24-27; fp predictor decode+encode 40; assemble/extract/byteswap 34. ZERO float64 locals in any kernel (all operate on uint8/16/32/64 byte buffers), and every launch is 1-D (inflate/LZW are [n_tiles, 32]), so the >20-float64-locals and >16x16-block rules cannot trigger. Probes this pass (PYTHONPATH=worktree, xrspatial.__file__ verified): dask read 2560x2560 chunks=256 = 400 tasks/100 chunks (4/chunk), max fan-in 10, 50k cap VERIFIED FIRING at chunks=8; streaming write peak scales with the BAND not the array (_stream_row_bands 2560 rows -> 4 computes, max band 768 rows, budget = streaming_buffer_bytes//bytes_per_src_row, default 256MB); eager GPU read cupy parity 0.0 in 6.5ms vs 24ms CPU (no host round trip); dask+cupy read lazy 4.02 tasks/chunk, corner cupy parity 0.0; cupy and dask+cupy WRITE round-trip parity 0.0. Cat5: no @jit without nopython anywhere in geotiff (21 @ngjit sites, ngjit=jit(nopython,nogil)); no np.zeros_like at all (all np.empty); _overview_kernels row-major, float64-stable accumulators. Prior deferred LOWs unchanged and NOT promoted (their files did not change since Pass 17). NEW deferred LOWs (documented, not fixed): (a) _decode.py:354/584 walk the WHOLE byte_counts tuple in Python per windowed read (once per dask chunk task) for the max_tile_bytes DoS guard -- O(n_tiles * n_chunks); max() would be the same check at C speed, ~10% on very large COGs only; same shape for _has_sparse at _decode.py:418/623 and _cog_http.py:636/879. (b) _writers/eager.py:1242 np.asarray(data) on a BARE dask array (not a DataArray) silently materialises -- out of the documented xr.DataArray|np.ndarray signature, sibling helpers test .compute() first. (c) _gpu_decode.py:1902/2200 _xp_byteswap allocates a full extra device image on big-endian GPU reads with predictor!=2, while the allocation-free in-place _swap_byte_lanes already exists and the predictor=2 path uses it. (d) _symbology.py:197-199 _eager_finite_stats does mask + arr[mask] compaction (~2.25x peak for float32) where nan-aware reductions would avoid the copy. (e) _overview.py:431 np.where on an already-owned float64 buffer (in-place assignment would do). (f) _compression.py:760/821 fp-predictor row kernels allocate np.empty(n) per row inside the ngjit body. (g) _pam.py _parse_rat can allocate two 1e6-element lists for a pathological sparse RAT (bounded, adversarial input only). SAFE/IO-bound holds. | Pass 18 (2026-07-05): re-audit, 0 CRITICAL/0 HIGH/0 MEDIUM. Only 2 commits since Pass 17: #3604 (f0275b61, O(1) ndim gate in _writers/eager.py before dispatch -- no perf surface) and #3605 (6b0887fb, tests only). Category scans across the whole subpackage came back matching Pass 17's known-clean state: Cat1 .compute()/.get() sites are the documented streaming writes (_writer.py row bands, _writers/gpu.py per-band stream #3166/#3117) and batched D2H; Cat2 50k-task cap intact at _backends/dask.py:580, #3597 StreamingStats single-pass intact, _dask_finite_stats fuses 5 reductions into one dask.compute; Cat3 all Device().synchronize() sites batch-end/error-recovery (#2212/#2107 fixes intact); Cat5 _overview_kernels.py @ngjit row-major float-stable. Probes on-device this pass (PYTHONPATH=worktree, verified xrspatial.__file__): dask read 2560x2560 deflate chunks=256 = 400 tasks/100 chunks (4/chunk); eager numpy parity 0.0; eager GPU read returns cupy, parity 0.0, 400ms incl warmup, no host round trip; dask+cupy lazy (401 tasks), corner chunk is cupy, parity 0.0. _gpu_decode.py untouched since Pass 17 so recorded register values stand (inflate 67 / LZW 29; introspection still unavailable via numba _Kernel API). Prior deferred LOWs unchanged (redundant .copy() in _writer predictor encode, identity pack arithmetic, twice-built IFD byte length, per-level device sync in _block_reduce_2d_gpu, _nvcomp_batch_compress adler32 staging, _CloudSource.read_range per-call open). No HIGH/MEDIUM => no /rockout, no new benchmark required. SAFE/IO-bound holds. | Pass 17 (2026-07-02): re-audit, 0 CRITICAL/0 HIGH/0 MEDIUM. Prior row logged as needing GPU re-validation; CUDA IS available on this host and all GPU decode/read paths were re-run on-device this pass (no cuda-unavailable token). Probes: dask read 2560x2560 chunks=256 = 400 tasks / 100 chunks (4/chunk, 50k cap intact at _backends/dask.py:580); eager numpy parity 0.0; eager GPU read of deflate-tiled file returns cupy, parity 0.0, 165ms, no host round trip; dask+cupy read lazy (401 tasks) producing cupy corner parity 0.0; GPU deflate+fp-predictor(predictor=3) decode parity 0.0 (routes through nvcomp GPU path on this host, numba _inflate_tiles fallback not specialized). Cat1: no .values/np.asarray/np.array on dask or cupy arrays in _backends/; the _writer.py:1221-1394 .compute() calls are the intentional streaming row-band writes (single-pass verified in #3117/#3235/#3597, not per-chunk re-execution). Cat3: all Device().synchronize() in _gpu_decode.py (1008/1029/1264/1758/2531/2824) are batch-end or error-recovery, not per-tile in loops -- #2212/#2107 fixes intact. Cat5: _overview_kernels.py are @ngjit CPU kernels (row-major inner loop, float-stable accumulators), not CUDA. Reviewed the 6 commits since last inspection (#3595/#3599 PAM overwrite cleanup, #3592/#3598 + #3593/#3596 docstrings, #3594 non-finite RAT guard, #3600 color_ramp streaming accumulation = Pass 16's own fix, #3588/#3589 isort): all correctness/docs, no new performance surface, no regressions. Register introspection unavailable in this numba's _Kernel API; kernels unchanged since Pass 7 (inflate 67 regs / LZW 29 regs). Prior deferred LOWs unchanged (redundant .copy() in _writer predictor encode, identity pack arithmetic, twice-built IFD byte length). No HIGH/MEDIUM => no /rockout, no new benchmark required. SAFE/IO-bound holds. | Pass 16 (2026-07-01): 1 MEDIUM found and fixed, 0 HIGH. to_geotiff(dask, color_ramp=...) executed the source graph twice: streaming write computed every chunk, then _write_sidecars -> _finite_stats ran dask.compute over the same source for the PAM/QML statistics (measured 32 chunk executions for a 16-chunk source; color_ramp_range escape hatch stays at 16). Filed #3597, fixed by threading a chunk_observer through _write_streaming's three materialisation sites (row band, wide-raster segment, strip band) into a new StreamingStats accumulator (_symbology.py; Chan mean/M2 combine, float64 accumulators, ddof=0, nodata/finite exclusion matches _finite_stats; all-valid buffers skip the boolean-index copy and reductions use dtype=float64 on the original buffer to avoid an astype copy). Post-fix executions 16/16; round-trip bench 8192x8192 f32 deflate write+stats 1.53s -> 1.44s and total source reads halve (the real win is expensive-to-recompute sources: HTTP COGs, cold storage, long pipelines). GPU writer and VRT paths keep the documented second pass (GPU writer fully materialises anyway); docstrings updated. 13 new tests in test_color_ramp_single_pass_3597.py incl. execution counters (row-band, strip, segmented wide path with full-width chunks), nodata/int/all-NaN/multiband gating, accumulator-vs-_finite_stats parity, and a dask+cupy gpu=False streaming leg (run on-device). Audited all 26 geotiff commits since 2026-06-11: pack range guards #3272/#3277 correctly defer per-chunk scans into the write's single compute (no #3235 regression); #3374 fixed Pass 14's deferred LOW (chunked GPU read now parses header/IFDs once); xarray engine #3375/#3377/#3380 is a thin wrapper (no eager compute); PAM sidecar read on open_geotiff is a local os.path.exists probe (no HTTP cost); _CloudSource cat_file change is neutral-to-better. Probes on-device this pass: dask read 400 tasks/100 chunks (4/chunk, cap intact), eager GPU read cupy parity 0.0 (222ms), dask+GPU lazy with cupy meta, symbology _finite_stats GPU/CPU parity 0.0. LOWs unchanged from prior passes. SAFE/IO-bound holds. | Pass 15 (2026-06-11): 1 MEDIUM found and fixed. _pack (_attrs.py:~1795) guarded the no-sentinel integer restore with an eager bool(out.isnull().any()), which executed the whole upstream dask graph at to_geotiff(pack=True) call time; the streaming writer then executed it again, so every source chunk computed twice (measured 32 decode-task executions for 16 chunks on a 512x512 int16 SCALE/OFFSET no-GDAL_NODATA source; 71->33 total task starts post-fix). Filed #3235, fixed by mapping a per-chunk NaN guard (_pack_guard_no_nan) into the graph for dask-backed data (raises from the write's single compute; numpy keeps the eager call-time check; meta= preserves cupy backing). 9 new tests in test_pack_lazy_nan_guard_3235.py incl. fusion-proof execution counter and cupy-chunk guard unit test (dask+cupy e2e still blocked upstream by #3112). Scrutinised all 16 commits since 2026-06-08 (pack/unpack series #3065/#3075/#3079/#3129/#3174/#3175, VRT placement #3135, compression_level gate #3176, streaming banding #3136, dask+cupy writer order fix #3171): no other regressions; #3171's get-then-asarray order is intentional D2H for gpu=False. GPU validated on-device this pass: eager GPU unpack returns cupy with exact parity (387ms incl warmup, only 0-d scalar .get()s -- no bulk host round trip), dask+GPU unpack lazy (112 tasks/16 chunks, cupy meta, compute returns cupy, parity 0.0), GDS fast path intact without unpack (4 tasks/chunk); unpack disqualifying GDS is documented intentional. Dask CPU probe 4 tasks/chunk, 50k-task cap intact. Note: #1714 (_write_vrt_tiled synchronous scheduler) is now FIXED+CLOSED (scheduler='threads' at _writers/eager.py:1517) -- drop from the open-issue list. LOW noted (no fix): _pack does identity (data-0.0)/1.0 arithmetic allocating two full-array temporaries when scale==1/offset==0 (masked_nodata-only pack); prior deferred LOWs unchanged. SAFE/IO-bound holds. | Pass 14 (2026-06-09): MEDIUM found and fixed -- _write_streaming ran one dask .compute() per 256-row tile-row/strip, so a source chunk taller than the band re-executed once per band it overlapped (measured 2x at chunks=512, 4x at chunks=1024, whole upstream graph re-runs for computed pipelines). Filed #3117, fixed via _stream_row_bands: consecutive tile-rows/strips group into row bands sized by the source chunk-row span (one-chunk halo, #3007 accounting) under streaming_buffer_bytes; each band computes once and tiles/strips are carved from the materialised band. Wide rasters needing column segmentation keep the per-tile-row path. Post-fix per-chunk executions == 1 on the default read->write round trip. 5 new tests (TestRowBandRecompute3117 + _stream_row_bands unit); write/integration/parity suites pass (2195). LOW deferred (no fix): _read_geotiff_gpu_chunked parses header+all IFDs twice at graph build (_backends/gpu.py ~1367-1419, cap check then GDS probe; build-time only). GPU paths validated on-device this pass: eager gpu read returns cupy with parity, dask+GPU chunked read lazy (17 tasks/4 chunks) with parity; GPU writer full materialisation is documented intentional (streaming_buffer_bytes no-op). Read path keeps 50k-task graph cap; dask read probe 4 tasks/chunk. SAFE/IO-bound holds. | Pass 13 (2026-05-20): 1 MEDIUM found and fixed. _nvjpeg_batch_encode (_gpu_decode.py:~L1560) and _nvjpeg2k_batch_encode (~L2958) called cupy.cuda.Device().synchronize() inside the per-tile encode loops, a whole-device fence that blocked every CUDA stream and serialised concurrent work (e.g. predictor encodes on other streams). The decode-side counterpart _try_nvjpeg_batch_decode already used cupy.cuda.Stream.null.synchronize() at L1442; the encoder side was inconsistent. Filed #2212 and fixed both encoders to use Stream.null.synchronize(), scoping the per-tile sync to the default stream the encode/retrieve calls were issued on. nvJPEG / nvJPEG2000 encoders maintain a single shared state per encoder so encodes within a batch are inherently serial; the fix removes the device-wide blocker without changing the API ordering contract. 5 new tests in test_nvjpeg_encode_stream_sync_2212.py (AST checks that neither encoder contains Device().synchronize() inside a for-loop, that both call Stream.null.synchronize() in the loop, and that the decoder reference pattern stays pinned). All 5 new tests + 19 existing related encode/decode tests pass. nvjpeg/nvjpeg2k shared libs not present on this host so end-to-end encode verification is gated; add cuda-unavailable-libs note to re-validate on a host with the RAPIDS conda env. SAFE/IO-bound verdict holds; no change in dask graph cost. Dask probe: 2560x2560 deflate-tiled file via read_geotiff_dask(chunks=256) yields 400 tasks for 100 chunks (4 tasks/chunk), well under the 50K cap. LOW deferred (no fix in this PR): _build_ifd called twice per IFD level in _assemble_standard_layout (_writer.py:1531+1543), _assemble_cog_layout (1582+1625), and the COG overview path (2519+2546+2740) -- the first call's bytes are discarded; only the overflow byte length is used to compute pixel_data_offset. Cost is bounded by IFD count (typically 1-5 overview levels) so absolute impact is minor. Pre-existing pattern. | Pass 12 (2026-05-18): 1 MEDIUM found and fixed. _try_nvjpeg2k_batch_decode at _gpu_decode.py:~L2725-2778 allocated per-tile per-component cupy.empty buffers (N*S round-trips through the cupy memory pool) and called cupy.cuda.Device().synchronize() once per tile, forcing default-stream serialisation that defeats nvJPEG2000's internal pipelining. Filed #2107 and fixed: pre-allocate a single d_comp_pool sized n_tiles*samples*tile_height*pitch under a _check_gpu_memory guard, derive per-tile/per-component views as slab offsets, and replace the per-tile sync with a single batch-end sync. Same pattern as #1659 (_try_nvcomp_from_device_bufs), #1688 (_try_kvikio_read_tiles), #1712 (_nvcomp_batch_compress). 7 new tests in test_nvjpeg2k_single_alloc_2107.py: AST-level structural assertions confirm no cupy.empty inside the for-loop and no Device().synchronize() inside the loop, plus pool/per_tile_comp_bytes presence and _check_gpu_memory guard checks; lib-absent short-circuit; unsupported-dtype cleanup contract; cupy-only pool slab-non-overlap test (gpu-marked). libnvjpeg2k.so not present on this host so the end-to-end nvJPEG2000 decode is gated -- note added to re-validate on a host with the RAPIDS conda env. All 30 jpeg2000/compression tests + 7 new tests pass. SAFE/IO-bound verdict holds (no change in dask graph cost). Dask probe: 4096x4096 deflate-tiled file via read_geotiff_dask(chunks=512) yields 256 tasks for 64 chunks (4 tasks/chunk), well under the 50K cap. | Pass 11 (2026-05-18): 1 MEDIUM found and fixed. _read_strips (_reader.py:~L1972) and _fetch_decode_cog_http_strips (_reader.py:~L2670) decoded strips sequentially in a Python for-loop while the tile counterparts (_read_tiles L2146, _fetch_decode_cog_http_tiles L2898) gated parallel decode on _PARALLEL_DECODE_PIXEL_THRESHOLD via ThreadPoolExecutor. Filed #2100 and fixed: both strip paths now collect jobs, parallel-decode when n_strips > 1 and strip_pixels >= 64K, then place sequentially. Measured (uint16, 4-core): 4096x4096 deflate 130ms->34ms (3.82x), 8192x8192 deflate 531ms->146ms (3.63x), 8192x8192 zstd 211ms->85ms (2.48x), uncompressed 25ms->22ms (1.14x). 5 new tests in test_parallel_strip_decode_2100.py (parallel/serial parity, pool-engaged on multi-strip, serial-path for single-strip, windowed cross-strip read, HTTP COG strip parity). 3998 tests pass; 8 pre-existing failures predating this change (predictor2 BE + size_param_validation_gpu_vrt reference now-private read_to_array attr). SAFE/IO-bound verdict holds. | Pass 10 (2026-05-15): 1 new MEDIUM found and fixed; 2 LOW noted. MEDIUM (_reader.py:2737): _fetch_decode_cog_http_tiles decoded tiles sequentially in a Python for-loop after the concurrent fetch landed (issue #1480). Local _read_tiles parallelises decode whenever tile_pixels >= 64K via ThreadPoolExecutor (_reader.py:2017); the HTTP path was structurally similar but never picked up the same gate, so wide windowed reads of multi-tile COGs left deflate/zstd decode single-threaded. Mirrored the local-path threshold + pool. 5 new tests in test_cog_http_parallel_decode_2026_05_15.py (parallel + serial round-trip correctness, pool-instantiation branch selection above the threshold, single-tile path skips the pool, structural _decode_strip_or_tile call count == n_tiles). All 262 COG/HTTP tests pass; 3162 of 3164 selected geotiff tests pass overall (2 pre-existing failures predating Pass 9 per prior notes -- test_predictor2_big_endian_gpu_1517 references the now-private read_to_array attr, and the test_size_param_validation_gpu_vrt_1776 tile_size=4 validator failure). LOW deferred (no fix in this PR): (1) _block_reduce_2d_gpu (_gpu_decode.py:3142/3163/3189) does bool(mask.any().item()) per overview level when nodata is set, paying one device sync per level; the alternative (unconditional cupy.putmask) always pays the work cost and the short-circuit is correct under the current API. (2) _nvcomp_batch_compress adler32 staging (_gpu_decode.py:2543-2546) issues n_tiles slice-assign kernels into a fresh contig buffer despite all callers passing slices of a single underlying d_tile_buf; an API refactor to accept the source buffer directly would skip the rebuild. SAFE/IO-bound verdict holds. Dask probe: 2560x2560 chunks=256 yields 400 tasks (4 per chunk), well under the 50000 cap. GPU probe: 1024x1024 float32 zstd read returns CuPy-backed in 236 ms with no host round-trip. | Rockout 2026-05-15: LOW filed #1934 -- _apply_nodata_mask_gpu used cupy.where (allocating); switched to cupy.putmask on the already-owned buffer (float path) and on the post-astype float64 buffer (int path). Saves one chunk-sized device allocation per call. 7 new tests in test_apply_nodata_mask_gpu_inplace_1934.py; 52 related nodata tests pass. | Pass 8 (2026-05-12): 1 new MEDIUM found and fixed. _assemble_standard_layout/_assemble_cog_layout returned bytes(bytearray), doubling peak memory transiently during eager writes. Filed #1756, fixed by returning the bytearray directly. Measured: 95 MB uint8 raster peak drops 202 MB -> 107 MB. _write_bytes / parse_header already accepted the buffer protocol so the change is transparent to callers. 6 new tests in test_assemble_layout_no_bytes_copy_1756.py. 2123 existing geotiff tests pass; the 10 unrelated failures (test_no_georef_windowed_coords_1710, test_predictor2_big_endian_gpu_1517) reference the now-private read_to_array attribute (commit 8adb749, issue #1708) and predate this change. SAFE/IO-bound verdict holds. | Pass 7 (2026-05-12): re-audit identified 4 MEDIUM findings, all real, all backed by microbenches. (1) unpack_bits sub-byte loops for bps=2/4/12 in _compression.py:836-878 were 100-200x slower than vectorised numpy (filed #1713, fixed in this branch: bps=4 2M pixels drops from 165ms to 3ms = 55x; bps=2/12 similar). (2) _write_vrt_tiled at __init__.py:1708 uses scheduler='synchronous' on independent tile writes; measured 33% slowdown on 256-tile zstd write vs threads scheduler (filed #1714, no fix yet). (3) _nvcomp_batch_compress at _gpu_decode.py:2522-2526 still does per-tile cupy.get().tobytes() despite #1552 / #1659 fixing the same pattern elsewhere; measured 45% reduction with concat+single get on n=1024 (filed #1712, no fix yet). (4) _nvcomp_batch_compress at _gpu_decode.py:2457 uses per-tile cupy.empty allocations; 1024 tiles 16KB drops from 4.7ms to 1.0ms with single contiguous + views (bundled into #1712). Cat 6 OOM verdict: SAFE/IO-bound holds -- read_geotiff_dask caps task count at _MAX_DASK_CHUNKS=50_000 and per-chunk memory is bounded by chunk size. _inflate_tiles_kernel resource usage on Ampere: 67 regs/thread, 2896B local/thread, 8192B shared/block (LZW kernel: 29 regs, 24576B shared) -- register pressure under control; high local memory in inflate is unavoidable (LZ77 state) but only thread 0 in each block uses it. | Pass 4 (2026-05-10): re-audit after #1559 (centralise attrs across all read backends). New _populate_attrs_from_geo_info helper at __init__.py:301 runs once per read, not per-chunk -- no perf impact. Probe: 2560x2560 deflate-tiled file opened via read_geotiff_dask yields 400 tasks (4 tasks/chunk for 100 chunks), well under 1M cap. read_geotiff_gpu(1024x1024) returns cupy.ndarray end-to-end with no host round-trip (226ms incl. write+decode). No new HIGH/MEDIUM findings. SAFE/IO-bound holds. | Pass 3 (2026-05-10): SAFE/IO-bound. Audited 4 perf commits: #1558 (in-place NaN writes on uniquely-owned buffers correct), #1556 (fp-predictor ngjit ~297us/tile for 256x256 float32), #1552 (single cupy.concatenate + one .get() for batched D2H at _gpu_decode.py:870-913), #1551 (parallel decode threshold >=65536px engages 256x256 default at _reader.py:1121). Bench: 8192x8192 f32 deflate+pred2 256-tile write 782ms; 4096x4096 f32 deflate read 83ms with parallel decode. Deferred LOW (none filed, all <10% MEDIUM threshold): _writer.py:459/1109 redundant .copy() before predictor encode (~1% per tile), _compression.py:280 lzw_decompress dst[:n].copy() (~2% per LZW tile decode), _writer.py:1419 seg_np.copy() before in-place NaN substitution (negligible, conditional path), _CloudSource.read_range opens fresh fsspec handle per range (pre-existing, predates audit scope). nvCOMP per-tile D2H batching break-even confirmed (variable sizes need staging buffer, no win). | Pass 3 (2026-05-10): audited f157746,39322c3,f23ec8f,1aac3b7. All 5 commits correct. Redundant .copy() in _writer.py:459,1109 and _compression.py:280 (1-2% overhead, LOW). _CloudSource.read_range() per-call open is pre-existing arch issue. No HIGH/MEDIUM regressions. SAFE. | re-audit 2026-05-02: 6 commits since 2026-04-16 (predictor=3 CPU encode/decode, GPU predictor stride fix, validate_tile_layout, BigTIFF LONG8 offsets, AREA_OR_POINT VRT, per-tile alloc guard). 1M dask chunk cap intact at __init__.py:948; adler32 batch transfer intact at _gpu_decode.py:1825. New code is metadata validation and dispatcher logic with no extra materialization or per-tile sync points. No HIGH/MEDIUM regressions. | Pass 5 (2026-05-12): re-audit identified MEDIUM in _gpu_decode.py:1577 _try_nvcomp_from_device_bufs: per-tile cupy.empty + trailing cupy.concatenate doubled peak VRAM and added serial concat. Filed #1659 and fixed to single-buffer + pointer offsets (matches LZW/deflate/host-buffer patterns at L1847/L1878/L1114). Microbench (alloc+concat overhead only, not full nvCOMP latency): n=256 tile_bytes=65536 drops 3.66ms->0.69ms, n=256 tile_bytes=262144 drops 8.18ms->0.13ms. Tests: 5 new tests in test_nvcomp_from_device_bufs_single_alloc_1659.py (codec short-circuit, no-lib short-circuit, memory-guard contract, real ZSTD round-trip via nvCOMP, structural single-buffer check). 1458 existing geotiff tests pass, 3 unrelated matplotlib/py3.14 failures pre-existing. SAFE/IO-bound verdict holds. | Pass 6 (2026-05-12): re-audit on top of #1659. New HIGH in _try_kvikio_read_tiles at _gpu_decode.py:941: per-tile cupy.empty() + blocking IOFuture.get() inside loop serialised GDS reads to ~1 outstanding pread, missed parallelism the kvikio worker pool was designed for, paid per-tile cupy.empty setup (matches #1659 anti-pattern in nvCOMP path), and lacked _check_gpu_memory guard. Filed #1688 and fixed to single contiguous buffer + batched submit + guard. Microbench with 8-worker pool simulation: 256 tiles@1ms latency drops 256ms->38.7ms (~6.6x); single-thread simulation 256ms->28.5ms (9x). Tests: 9 new tests in test_kvikio_batched_pread_1688.py (kvikio-absent path, single-buffer pointer arithmetic, submit-before-get ordering, memory guard, partial-read fallback, round-trip data, zero-size/all-sparse tiles). All 1577 geotiff tests pass except pre-existing matplotlib/py3.14 failures." glcm,2026-03-31T18:00:00Z,SAFE,compute-bound,0,,"Downgraded to MEDIUM. da.stack without rechunk is scheduling overhead, not OOM risk." gpu_rtx,2026-07-23,N/A,compute-bound,2,3691,"rtx-available (OptiX 9.1/A6000); 2 HIGH fixed via #3691 PR: mesh_utils datahash did full D2H copy (hash(str(data.get())), 27ms = 10% of hillshade_rtx/viewshed_gpu @4000x4000) -> edge+strided sample hash (0.14ms); _triangulate_terrain 100-block launch loop (38.2ms, 157 launches + NumbaPerformanceWarning) -> single launch (1.5ms). End-to-end 288->211ms hillshade, 278->196ms viewshed. No dask path -> OOM N/A; _memory.py guards device buffers at 120B/px vs 50% free. Register pressure verified fine (24-40 regs/thread @32x32). LOW not fixed: free_all_blocks() after mesh build forces re-cudaMalloc (intentional BVH headroom); fresh RTX() per call means mesh hash cache never hits (design); viewshed float64 astype adds 8B/px beyond budget but within headroom." hillshade,2026-04-16T12:00:00Z,SAFE,compute-bound,0,,"Re-audit after Horn's method rewrite (PR 1175): clean stencil, map_overlap depth=(1,1), no materialization. Zero findings." diff --git a/benchmarks/benchmarks/geotiff.py b/benchmarks/benchmarks/geotiff.py index b52260465..8914957da 100644 --- a/benchmarks/benchmarks/geotiff.py +++ b/benchmarks/benchmarks/geotiff.py @@ -101,6 +101,42 @@ def time_write_cog(self, nx, backend): compression="zstd", gpu=self.gpu) +# ------------------------------------------------------------------------- +# COG write with symbology from a dask source (#3695) +# +# cog=True skips the streaming writer, so the color_ramp statistics used to +# be taken by re-running the source graph after the pixels were already +# materialised: the caller's whole pipeline executed twice. The fix folds +# the materialised buffer into the same StreamingStats accumulator the +# streaming path uses. The dask parameter is the one that regresses; numpy +# is kept as the no-graph control, and the no-ramp case below is the +# baseline the ramp variant should stay close to. +# ------------------------------------------------------------------------- + +class WriteCOGSymbology: + params = ([512, 2048], ["numpy", "dask"]) + param_names = ("nx", "backend") + + def setup(self, nx, backend): + ny = nx // 2 + self.da = _make_dataarray(ny, nx, backend) + self.dir = tempfile.mkdtemp(prefix="asv_geotiff_cog_sym_") + self.path = os.path.join(self.dir, f"cogsym_{nx}_{backend}.tif") + + def teardown(self, nx, backend): + shutil.rmtree(getattr(self, "dir", ""), ignore_errors=True) + + def time_write_cog_color_ramp(self, nx, backend): + to_geotiff(self.da, self.path, cog=True, color_ramp="viridis", + overview_levels=[2, 4, 8], overview_resampling="mean", + compression="zstd") + + def time_write_cog_no_ramp(self, nx, backend): + to_geotiff(self.da, self.path, cog=True, + overview_levels=[2, 4, 8], overview_resampling="mean", + compression="zstd") + + # ------------------------------------------------------------------------- # Eager read path (open_geotiff) # diff --git a/xrspatial/geotiff/_writers/eager.py b/xrspatial/geotiff/_writers/eager.py index 2c66ba3cb..30a01f12a 100644 --- a/xrspatial/geotiff/_writers/eager.py +++ b/xrspatial/geotiff/_writers/eager.py @@ -445,11 +445,12 @@ def to_geotiff(data: xr.DataArray | np.ndarray, ``attrs['category_names']`` -- those get the RAT sidecar instead), a multiband array, a file-like destination, or data with no finite values. Computing the statistics is an extra reduction pass over the - data. The streaming dask write accumulates them from the buffers it - materialises anyway, so the source graph still runs once; the GPU - (``gpu=True``) and VRT (``.vrt``) write paths execute a dask source - a second time for the statistics (see ``color_ramp_range`` to skip - that). Ignored when ``pack=True``, whose on-disk packed values would + data. Both the streaming dask write and the ``cog=True`` write + accumulate them from the buffers they materialise anyway, so the + source graph still runs once; the GPU (``gpu=True``) and VRT + (``.vrt``) write paths execute a dask source a second time for the + statistics (see ``color_ramp_range`` to skip that). Ignored when + ``pack=True``, whose on-disk packed values would not match a ramp built from the logical values. Every string-path write refreshes the PAM ``.aux.xml``: a sidecar left by a previous write at the same path is removed and re-created only when this @@ -1224,6 +1225,30 @@ def _write_sidecars(): arr = raw.compute() # Dask -> numpy if hasattr(arr, 'get'): arr = arr.get() # Dask+CuPy -> numpy + # ``color_ramp`` statistics off the buffer we just + # materialised. ``cog=True`` skips the streaming writer, so + # the ``chunk_observer`` accumulation added for issue #3597 + # never runs and ``_write_sidecars`` would fall back to + # ``_finite_stats`` on the still-lazy ``_sym_data`` -- a + # second full execution of the caller's graph (issue #3695). + # Feed the accumulator the materialised array instead, + # before the sentinel restore below rewrites NaN, so the + # statistics describe the same logical values + # ``_finite_stats`` would have seen. + # + # ``StreamingStats`` accumulates the moments in float64 + # where ``_finite_stats`` accumulates at the input's native + # width, so a float32 source's mean / stddev shift by ~1e-7 + # relative. That is the accumulator being more accurate, and + # it makes a ``cog=True`` write agree with the ``cog=False`` + # streaming write on the same data rather than diverge. + # A file-like destination never reaches here with symbology + # pending: ``_sym_stops`` is only set for a string path. + if _sym_stops is not None and color_ramp_range is None: + from .._symbology import StreamingStats, _is_single_band + if _is_single_band(data): + _sym_stream_stats = StreamingStats(nodata=_sym_nodata) + _sym_stream_stats.update(arr) else: arr = np.asarray(raw) # Reject ambiguous 3D layouts. The validator runs diff --git a/xrspatial/geotiff/tests/write/test_color_ramp_cog_single_pass_3695.py b/xrspatial/geotiff/tests/write/test_color_ramp_cog_single_pass_3695.py new file mode 100644 index 000000000..cdb1db726 --- /dev/null +++ b/xrspatial/geotiff/tests/write/test_color_ramp_cog_single_pass_3695.py @@ -0,0 +1,324 @@ +"""``color_ramp`` statistics on the eager COG write path (#3695). + +#3597 stopped ``to_geotiff(dask_data, path, color_ramp=...)`` from executing +the source graph twice, but only on the streaming path. ``cog=True`` skips +the streaming writer (COG overviews need the full array), so it kept falling +through to ``_finite_stats`` on the still-lazy DataArray and ran the caller's +whole pipeline a second time. + +The fix folds the already-materialised buffer into a ``StreamingStats`` +accumulator, which ``write_symbology_sidecars`` prefers over the reduction. +These tests pin the single execution with a counting ``map_blocks`` layer and +check the resulting sidecar values against ``_finite_stats`` on every branch +the accumulator now covers. +""" +import os +import threading +import xml.etree.ElementTree as ET + +import numpy as np +import pytest +import xarray as xr + +from xrspatial.geotiff import to_geotiff +from xrspatial.geotiff._symbology import _eager_finite_stats + +from .._helpers.markers import requires_gpu + +pytest.importorskip("tifffile") + + +def _counting_da(base, chunks, dims=("y", "x"), coords=None, attrs=None): + """Dask-backed DataArray whose chunks count their own executions.""" + import dask.array as dsa + + counter = {"n": 0} + lock = threading.Lock() + + def _count(block): + if block.size: # skip dask's zero-size meta-inference call + with lock: + counter["n"] += 1 + return block + + arr = dsa.from_array(base, chunks=chunks).map_blocks( + _count, dtype=base.dtype) + if coords is None: + coords = {"y": np.arange(base.shape[0], dtype="float64"), + "x": np.arange(base.shape[1], dtype="float64")} + da = xr.DataArray( + arr, dims=dims, coords=coords, + attrs={"crs": 4326, **(attrs or {})}, + ) + return da, counter + + +def _aux_stats(path): + """Parse ``.aux.xml`` into a ``{STATISTICS_*: float}`` dict.""" + band = ET.parse(path + ".aux.xml").getroot().find(".//PAMRasterBand") + return {mdi.get("key"): float(mdi.text) + for mdi in band.findall("./Metadata/MDI")} + + +def _ref_stats(arr, nodata=None): + """(min, max, mean, population std) over finite non-nodata values. + + Thin wrapper over the production reduction the accumulator replaces, so + the assertions below compare against the real thing rather than a + reimplementation. ``_eager_finite_stats`` is used directly because + ``_finite_stats`` unwraps ``.data``, which on a bare numpy array is a + memoryview. + """ + return _eager_finite_stats(np.asarray(arr), nodata) + + +def _assert_aux_matches_rel(path, ref, rel): + stats = _aux_stats(path) + assert stats["STATISTICS_MINIMUM"] == pytest.approx(ref[0], rel=rel) + assert stats["STATISTICS_MAXIMUM"] == pytest.approx(ref[1], rel=rel) + assert stats["STATISTICS_MEAN"] == pytest.approx(ref[2], rel=rel) + assert stats["STATISTICS_STDDEV"] == pytest.approx(ref[3], rel=rel) + + +def _assert_aux_matches(path, ref): + _assert_aux_matches_rel(path, ref, rel=1e-9) + + +_RNG = np.random.default_rng(3695) +_BASE = _RNG.uniform(-50.0, 150.0, (64, 64)) +_BASE[3, 7] = np.nan +_BASE[40, 2] = np.nan + +_N_CHUNKS = 16 # 64x64 at chunks=(16, 16) + + +# -------------------------------------------------------------------------- +# single execution of the source graph +# -------------------------------------------------------------------------- + +def test_cog_color_ramp_executes_source_once(tmp_path): + """The regression: cog=True used to run the graph twice.""" + da, counter = _counting_da(_BASE, chunks=(16, 16)) + path = str(tmp_path / "cog_once_3695.tif") + to_geotiff(da, path, cog=True, color_ramp="viridis") + assert counter["n"] == _N_CHUNKS # was 2 * _N_CHUNKS before the fix + _assert_aux_matches(path, _ref_stats(_BASE)) + assert os.path.exists(str(tmp_path / "cog_once_3695.qml")) + + +def test_cog_without_color_ramp_still_single_pass(tmp_path): + """The no-symbology COG write was already single-pass; keep it that way.""" + da, counter = _counting_da(_BASE, chunks=(16, 16)) + path = str(tmp_path / "cog_plain_3695.tif") + to_geotiff(da, path, cog=True) + assert counter["n"] == _N_CHUNKS + assert not os.path.exists(path + ".aux.xml") + + +def test_streaming_path_unchanged(tmp_path): + """cog=False keeps the #3597 chunk_observer behaviour.""" + da, counter = _counting_da(_BASE, chunks=(16, 16)) + path = str(tmp_path / "stream_3695.tif") + to_geotiff(da, path, color_ramp="viridis") + assert counter["n"] == _N_CHUNKS + _assert_aux_matches(path, _ref_stats(_BASE)) + + +# -------------------------------------------------------------------------- +# the accumulated statistics must equal the reduction they replaced +# -------------------------------------------------------------------------- + +def test_cog_stats_match_finite_stats_across_chunkings(tmp_path): + """Chan combine must be chunking-invariant, so all layouts agree.""" + ref = _ref_stats(_BASE) + for i, chunks in enumerate([(16, 16), (64, 64), (8, 32), (7, 13)]): + da, _ = _counting_da(_BASE, chunks=chunks) + path = str(tmp_path / f"chunking_{i}_3695.tif") + to_geotiff(da, path, cog=True, color_ramp="viridis") + _assert_aux_matches(path, ref) + + +def test_cog_stats_exclude_nodata(tmp_path): + base = _BASE.copy() + base[10:14, 10:14] = -9999.0 + da, counter = _counting_da(base, chunks=(16, 16)) + path = str(tmp_path / "nodata_3695.tif") + to_geotiff(da, path, cog=True, color_ramp="viridis", nodata=-9999.0) + assert counter["n"] == _N_CHUNKS + _assert_aux_matches(path, _ref_stats(base, nodata=-9999.0)) + # The sentinel must not leak into the ramp bounds. + assert _aux_stats(path)["STATISTICS_MINIMUM"] > -9999.0 + + +def test_cog_integer_source_stats(tmp_path): + base = _RNG.integers(0, 500, (64, 64)).astype("int32") + da, counter = _counting_da(base, chunks=(16, 16)) + path = str(tmp_path / "int_3695.tif") + to_geotiff(da, path, cog=True, color_ramp="viridis") + assert counter["n"] == _N_CHUNKS + _assert_aux_matches(path, _ref_stats(base)) + + +def test_cog_integer_source_with_nodata(tmp_path): + base = _RNG.integers(0, 500, (64, 64)).astype("int32") + base[0, :] = -1 + da, counter = _counting_da(base, chunks=(16, 16)) + path = str(tmp_path / "int_nodata_3695.tif") + to_geotiff(da, path, cog=True, color_ramp="viridis", nodata=-1) + assert counter["n"] == _N_CHUNKS + _assert_aux_matches(path, _ref_stats(base, nodata=-1)) + + +# -------------------------------------------------------------------------- +# gating: the accumulator must not change who gets sidecars +# -------------------------------------------------------------------------- + +def test_cog_color_ramp_range_still_skips_stats(tmp_path): + """The escape hatch writes bounds only and never builds the accumulator.""" + da, counter = _counting_da(_BASE, chunks=(16, 16)) + path = str(tmp_path / "cog_rng_3695.tif") + to_geotiff(da, path, cog=True, color_ramp="viridis", + color_ramp_range=(0.0, 10.0)) + assert counter["n"] == _N_CHUNKS + stats = _aux_stats(path) + assert stats["STATISTICS_MINIMUM"] == pytest.approx(0.0) + assert stats["STATISTICS_MAXIMUM"] == pytest.approx(10.0) + assert "STATISTICS_MEAN" not in stats + assert "STATISTICS_STDDEV" not in stats + + +def test_cog_multiband_gets_no_symbology(tmp_path): + base = _RNG.uniform(0.0, 1.0, (32, 32, 3)) + da, counter = _counting_da( + base, chunks=(16, 16, 3), dims=("y", "x", "band"), + coords={"y": np.arange(32, dtype="float64"), + "x": np.arange(32, dtype="float64"), + "band": np.arange(3)}) + path = str(tmp_path / "multiband_3695.tif") + to_geotiff(da, path, cog=True, color_ramp="viridis") + assert counter["n"] == 4 + assert not os.path.exists(path + ".aux.xml") + assert not os.path.exists(str(tmp_path / "multiband_3695.qml")) + + +def test_cog_single_band_3d_gets_symbology(tmp_path): + """A 3D array with one band is still a single-band raster.""" + base = _BASE[:32, :32].reshape(32, 32, 1) + da, counter = _counting_da( + base, chunks=(16, 16, 1), dims=("y", "x", "band"), + coords={"y": np.arange(32, dtype="float64"), + "x": np.arange(32, dtype="float64"), + "band": np.arange(1)}) + path = str(tmp_path / "single3d_3695.tif") + to_geotiff(da, path, cog=True, color_ramp="viridis") + assert counter["n"] == 4 + _assert_aux_matches(path, _ref_stats(base)) + + +def test_cog_all_nan_writes_no_sidecar(tmp_path): + base = np.full((32, 32), np.nan) + da, counter = _counting_da(base, chunks=(16, 16)) + path = str(tmp_path / "allnan_3695.tif") + to_geotiff(da, path, cog=True, color_ramp="viridis") + assert counter["n"] == 4 + # No finite values -> no statistics, matching _finite_stats returning None. + assert _ref_stats(base) is None + assert not os.path.exists(path + ".aux.xml") + assert not os.path.exists(str(tmp_path / "allnan_3695.qml")) + + +def test_cog_constant_raster_writes_stats_but_no_qml(tmp_path): + """vmin == vmax is a degenerate ramp; stats still land.""" + base = np.full((32, 32), 7.5) + da, _ = _counting_da(base, chunks=(16, 16)) + path = str(tmp_path / "const_3695.tif") + to_geotiff(da, path, cog=True, color_ramp="viridis") + stats = _aux_stats(path) + assert stats["STATISTICS_MINIMUM"] == pytest.approx(7.5) + assert stats["STATISTICS_MAXIMUM"] == pytest.approx(7.5) + assert not os.path.exists(str(tmp_path / "const_3695.qml")) + + +def test_cog_float32_matches_streaming_write(tmp_path): + """float32 statistics agree with the cog=False write on the same data. + + ``StreamingStats`` accumulates in float64 while ``_finite_stats`` + accumulates at the input's native width, so a float32 source's mean and + stddev move by ~1e-7 relative against the old COG output. What matters is + that the two dask write paths now agree with each other, which they did + not before: cog=False has used the float64 accumulator since #3597. + """ + base = _RNG.uniform(-1e4, 1e4, (64, 64)).astype("float32") + base[5, 5] = np.nan + + cog_path = str(tmp_path / "f32_cog_3695.tif") + stream_path = str(tmp_path / "f32_stream_3695.tif") + da_cog, counter = _counting_da(base, chunks=(16, 16)) + da_stream, _ = _counting_da(base, chunks=(16, 16)) + to_geotiff(da_cog, cog_path, cog=True, color_ramp="viridis") + to_geotiff(da_stream, stream_path, color_ramp="viridis") + + assert counter["n"] == _N_CHUNKS + assert _aux_stats(cog_path) == _aux_stats(stream_path) + # Still the same numbers as the native-width reduction to float32 + # resolution, so the ramp bounds and stretch are unchanged in practice. + _assert_aux_matches_rel(cog_path, _ref_stats(base), rel=1e-6) + + +def test_eager_numpy_cog_unaffected(tmp_path): + """A numpy source has no graph to re-execute; sidecars stay the same.""" + da = xr.DataArray( + _BASE, dims=("y", "x"), + coords={"y": np.arange(64, dtype="float64"), + "x": np.arange(64, dtype="float64")}, + attrs={"crs": 4326}) + path = str(tmp_path / "numpy_3695.tif") + to_geotiff(da, path, cog=True, color_ramp="viridis") + _assert_aux_matches(path, _ref_stats(_BASE)) + + +# -------------------------------------------------------------------------- +# GPU +# -------------------------------------------------------------------------- + +@requires_gpu +def test_cog_dask_cupy_color_ramp_single_pass(tmp_path): + """dask+cupy reaches the same eager COG fallthrough (gpu=False writer).""" + import cupy + import dask.array as dsa + + counter = {"n": 0} + lock = threading.Lock() + + def _count(block): + if block.size: + with lock: + counter["n"] += 1 + return block + + arr = dsa.from_array(cupy.asarray(_BASE), chunks=(16, 16)).map_blocks( + _count, dtype=_BASE.dtype) + da = xr.DataArray( + arr, dims=("y", "x"), + coords={"y": np.arange(64, dtype="float64"), + "x": np.arange(64, dtype="float64")}, + attrs={"crs": 4326}) + path = str(tmp_path / "cog_gpu_3695.tif") + to_geotiff(da, path, cog=True, color_ramp="viridis", gpu=False) + assert counter["n"] == _N_CHUNKS + _assert_aux_matches(path, _ref_stats(_BASE)) + + +@requires_gpu +def test_cog_cupy_eager_color_ramp(tmp_path): + """A plain cupy source has no graph; the accumulator must not fire.""" + import cupy + + da = xr.DataArray( + cupy.asarray(_BASE), dims=("y", "x"), + coords={"y": np.arange(64, dtype="float64"), + "x": np.arange(64, dtype="float64")}, + attrs={"crs": 4326}) + path = str(tmp_path / "cupy_3695.tif") + to_geotiff(da, path, cog=True, color_ramp="viridis", gpu=False) + _assert_aux_matches(path, _ref_stats(_BASE))