From 5635fd39ed72b0b1de7f75fd39bc9a141e2d360b Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Sun, 10 May 2026 22:04:49 -0700 Subject: [PATCH 1/6] Record geotiff security sweep pass 11 (2026-05-10) Audited the one new commit since pass 10: - #1559 (769bcf7): centralise GeoTIFF attrs population across all read backends via a single _populate_attrs_from_geo_info helper. Pure metadata-handling refactor. The helper dedupes attrs population across the eager numpy, dask, and two GPU read paths so all four backends emit the same key set. No new allocations, no file I/O changes, no kernel changes, no dtype handling changes -- the helper only writes to an in-memory attrs dict from already-validated geo_info fields. Re-verified all prior guards intact: - MAX_PIXELS_DEFAULT=1e9 at _reader.py:46 - MAX_IFD_ENTRY_COUNT=100_000 at _header.py:34 - MAX_IFDS=256 at _header.py:45 - MAX_TILE_BYTES_DEFAULT=256 MiB at _reader.py:68 - _MAX_DASK_CHUNKS=50_000 at __init__.py:1428 - realpath canonicalisation at _reader.py:161 and _vrt.py:160 Cat 1-6 all clean. No PRs opened. --- .claude/sweep-security-state.csv | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.claude/sweep-security-state.csv b/.claude/sweep-security-state.csv index 257f25a59..c91c1c745 100644 --- a/.claude/sweep-security-state.csv +++ b/.claude/sweep-security-state.csv @@ -18,7 +18,7 @@ fire,2026-04-25,,,,,"Clean. Despite the module's size hint, fire.py is purely pe flood,2026-05-03,1437,MEDIUM,3,,Re-audit 2026-05-03. MEDIUM Cat 3 fixed in PR #1438 (travel_time and flood_depth_vegetation now validate mannings_n DataArray values are finite and strictly positive via _validate_mannings_n_dataarray helper). No remaining unfixed findings. Other categories clean: every allocation is same-shape as input; no flat index math; NaN propagation explicit in every backend; tan_slope clamped by _TAN_MIN; no CUDA kernels; no file I/O; every public API calls _validate_raster on DataArray inputs. focal,2026-04-27,1284,HIGH,1,,"HIGH (fixed PR #1286): apply(), focal_stats(), and hotspots() accepted unbounded user-supplied kernels via custom_kernel(), which only checks shape parity. The kernel-size guard from #1241 (_check_kernel_memory) only ran inside circle_kernel/annulus_kernel, so a (50001, 50001) custom kernel on a 10x10 raster allocated ~10 GB on the kernel itself plus a much larger padded raster before any work -- same shape as the bilateral DoS in #1236. Fixed by adding _check_kernel_vs_raster_memory in focal.py and wiring it into apply(), focal_stats(), and hotspots() after custom_kernel() validation. All 134 focal tests + 19 bilateral tests pass. No other findings: 10 CUDA kernels all have proper bounds + stencil guards; _validate_raster called on every public entry point; hotspots already raises ZeroDivisionError on constant-value rasters; _focal_variety_cuda uses a fixed-size local buffer (silent truncation but bounded); _focal_std_cuda/_focal_var_cuda clamp the catastrophic-cancellation case via if var < 0.0: var = 0.0; no file I/O." geodesic,2026-04-27,1283,HIGH,1,,"HIGH (fixed PR #1285): slope(method='geodesic') and aspect(method='geodesic') stack a (3, H, W) float64 array (data, lat, lon) before dispatch with no memory check. A large lat/lon-tagged raster passed to either function would OOM. Fixed by adding _check_geodesic_memory(rows, cols) in xrspatial/geodesic.py (mirrors morphology._check_kernel_memory): budgets 56 bytes/cell (24 stacked float64 + 4 float32 output + 24 padded copy + slack) and raises MemoryError when > 50% of available RAM; called from slope.py and aspect.py inside the geodesic branch before dispatch. No other findings: 6 CUDA kernels all have bounds guards (e.g. _run_gpu_geodesic_aspect at geodesic.py:395), custom 16x16 thread blocks avoid register spill, no shared memory, _validate_raster runs upstream in slope/aspect, all backends cast to float32, slope_mag < 1e-7 flat threshold prevents arctan2 NaN propagation, curvature correction uses hardcoded WGS84 R." -geotiff,2026-05-10,,,,,"Pass 10 (2026-05-10): CLEAN. Audited 8 commits since pass 9: #1558 (drop nodata copies - buffers uniquely owned), #1556 (fp-predictor ngjit row loop - bit-exact), #1552 (batched GDS-fallback D2H - _check_gpu_memory guards new staging), #1551 (parallel decode at tile_size=256 - pool bounded by min(n_tiles, cpu_count), mmap thread-safe read-only), #1549 (nvjpeg constants - actually closed CRITICAL: prior off-by-two passed PLANAR where INTERLEAVED expected, dereferenced NULL channel ptrs), #1542 (GPU nodata mask - called once per path, stripped fallback returns early), #1543 (combined backend tests-only), #1544 (__all__/deprecation no security surface). Re-verified guards intact: MAX_PIXELS_DEFAULT, MAX_IFD_ENTRY_COUNT, MAX_IFDS, MAX_TILE_BYTES_DEFAULT (256 MiB), _MAX_DASK_CHUNKS, _check_gpu_memory, _MmapCache TOCTOU, VRT realpath, LZW (sp<4096) and Deflate (overflow array) bounds. No issues opened. | Pass 2 (2026-05-09): Found unbounded compressed-byte-count fetch in HTTP COG path (_fetch_decode_cog_http_tiles passes attacker-controlled TileByteCounts to read_ranges_coalesced). Verified exploit: 4 tiles x 100 MB byte_count -> single 100 MB Range GET. Fixed via per-tile cap (MAX_TILE_BYTES_DEFAULT=256 MiB, XRSPATIAL_COG_MAX_TILE_BYTES env override). Other recent PRs (#1530 IFD chain cap, #1527 IFD count, #1535 LERC mask, #1532 PlanarConfig=2, #1534 HTTP coalesce, #1531 parallel write, #1528 nvCOMP batch) audited clean. | Prior: Re-audit 2026-04-28: only #7984f7a since prior pass (predictor=3 multi-band CPU fix, LERC dedup, VRT nodata dtype, AREA_OR_POINT, BigTIFF LONG8 offsets). LZW/inflate GPU kernels gate every write on out_pos single 100 MB Range GET. Fixed via per-tile cap (MAX_TILE_BYTES_DEFAULT=256 MiB, XRSPATIAL_COG_MAX_TILE_BYTES env override). Other recent PRs (#1530 IFD chain cap, #1527 IFD count, #1535 LERC mask, #1532 PlanarConfig=2, #1534 HTTP coalesce, #1531 parallel write, #1528 nvCOMP batch) audited clean. | Prior: Re-audit 2026-04-28: only #7984f7a since prior pass (predictor=3 multi-band CPU fix, LERC dedup, VRT nodata dtype, AREA_OR_POINT, BigTIFF LONG8 offsets). LZW/inflate GPU kernels gate every write on out_pos= 3 and distance only as >= 1, with no upper bound on either. _glcm_numba_kernel iterates range(r-half, r+half+1) for every pixel, so window_size=1_000_001 on a 10x10 raster ran ~10^14 loop iterations with all neighbors failing the interior bounds check (CPU DoS). On the dask backends depth = window_size // 2 + distance drove map_overlap padding, so a huge window also caused oversize per-chunk allocations (memory DoS). Fixed by adding max_val caps in the public entrypoint: window_size <= max(3, min(rows, cols)) and distance <= max(1, window_size // 2). One cap covers every backend because cupy and dask+cupy call through to the CPU kernel after cupy.asnumpy. No other HIGH findings: levels is already capped at 256 so the per-pixel np.zeros((levels, levels)) matrix in the kernel is bounded to 512 KB. No CUDA kernels. No file I/O. Quantization clips to [0, levels-1] before the kernel and NaN maps to -1 which the kernel filters with i_val >= 0. Entropy log(p) and correlation p / (std_i * std_j) are both guarded. All four backends use _validate_raster and cast to float64 before quantizing. MEDIUM (unfixed, Cat 1): the per-pixel np.zeros((levels, levels)) allocation inside the hot loop is a perf issue (levels=256 -> 512 KB alloc+free per pixel) but not a security issue because levels is bounded. Could be hoisted out of the loop or replaced with an in-place clear, but that is an efficiency concern, not security." gpu_rtx,2026-04-29,1308,HIGH,1,,"HIGH (fixed #1308 / PR #1310): hillshade_rtx (gpu_rtx/hillshade.py:184) and viewshed_gpu (gpu_rtx/viewshed.py:269) allocated cupy device buffers sized by raster shape with no memory check. create_triangulation (mesh_utils.py:23-24) adds verts (12 B/px) + triangles (24 B/px) = 36 B/px; hillshade_rtx adds d_rays(32) + d_hits(16) + d_aux(12) + d_output(4) = 64 B/px (100 B/px total); viewshed_gpu adds d_rays(32) + d_hits(16) + d_visgrid(4) + d_vsrays(32) = 84 B/px (120 B/px total). A 30000x30000 raster asked for 90-108 GB of VRAM before cupy surfaced an opaque allocator error. Fixed by adding gpu_rtx/_memory.py with _available_gpu_memory_bytes() and _check_gpu_memory(func_name, h, w) helpers (cost_distance #1262 / sky_view_factor #1299 pattern, 120 B/px budget covers worst case, raises MemoryError when required > 50% of free VRAM, skips silently when memGetInfo() unavailable). Wired into both entry points after the cupy.ndarray type check and before create_triangulation. 9 new tests in test_gpu_rtx_memory.py (5 helper-unit + 4 end-to-end gated on has_rtx). All 81 existing hillshade/viewshed tests still pass. Cat 4 clean: all CUDA kernels (hillshade.py:25/62/106, viewshed.py:32/74/116, mesh_utils.py:50) have bounds guards; no shared memory, no syncthreads needed. MEDIUM not fixed (Cat 6): hillshade_rtx and viewshed_gpu do not call _validate_raster directly but parent hillshade() (hillshade.py:252) and viewshed() (viewshed.py:1707) already validate, so input validation runs before the gpu_rtx entry point - defense-in-depth, not exploitable. MEDIUM not fixed (Cat 2): mesh_utils.py:64-68 cast mesh_map_index to int32 in the triangle index buffer; overflows at H*W > 2.1B vertices (~46341x46341+) but the new memory guard rejects rasters that large first - documentation/clarity item rather than exploitable. MEDIUM not fixed (Cat 3): mesh_utils.py:19 scale = maxDim / maxH divides by zero on an all-zero raster, propagating inf/NaN into mesh vertex z-coords; separate follow-up. LOW not fixed (Cat 5): mesh_utils.write() opens user-supplied path without canonicalization but its only call site (mesh_utils.py:38-39) sits behind if False: in create_triangulation, not reachable in production." hillshade,2026-04-27,,,,,"Clean. Cat 1: only allocation is the output np.empty(data.shape) at line 32 (cupy at line 165) and a _pad_array with hardcoded depth=1 (line 62) -- bounded by caller, no user-controlled amplifier. Azimuth/altitude are scalars and don't drive size. Cat 2: numba kernel uses range(1, rows-1) with simple (y, x) indexing; numba range loops promote to int64. Cat 3: math.sqrt(1.0 + xx_plus_yy) is always >= 1.0 (no neg sqrt, no div-by-zero); NaN elevation propagates correctly through dz_dx/dz_dy -> shaded -> output (the shaded < 0.0 / shaded > 1.0 clamps don't fire on NaN). Azimuth validated to [0, 360], altitude to [0, 90]. Cat 4: _gpu_calc_numba (line 107) guards both grid bounds and 3x3 stencil reads via i > 0 and i < shape[0]-1 and j > 0 and j < shape[1]-1; no shared memory. Cat 5: no file I/O. Cat 6: hillshade() calls _validate_raster (line 252) and _validate_scalar for both azimuth (253) and angle_altitude (254); all four backend paths cast to float32; tests parametrize int32/int64/float32/float64." From c6f66ef6358b29e9d9d9d6619a48e818023478fd Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Sun, 10 May 2026 22:05:49 -0700 Subject: [PATCH 2/6] Record geotiff perf sweep pass 4 (2026-05-10) Re-audit after #1559 (centralise attrs across all read backends). No new HIGH/MEDIUM findings. SAFE/IO-bound holds. - New _populate_attrs_from_geo_info helper runs once per read, not per chunk -- no perf impact on hot paths. - Probe: read_geotiff_dask(2560x2560, chunks=256) yields 400 tasks (4 tasks/chunk for 100 chunks), well under 1M cap. - Probe: read_geotiff_gpu(1024x1024) returns cupy.ndarray end-to-end with no host round-trip (226ms incl. write+decode). --- .claude/sweep-performance-state.csv | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.claude/sweep-performance-state.csv b/.claude/sweep-performance-state.csv index 4b3ff160f..e332fcd79 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-03-31T18:00:00Z,SAFE,compute-bound,0,, focal,2026-03-31T18:00:00Z,SAFE,compute-bound,0,, geodesic,2026-03-31T18:00:00Z,N/A,compute-bound,0,, -geotiff,2026-05-10,SAFE,IO-bound,0,,"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." +geotiff,2026-05-10,SAFE,IO-bound,0,,"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." glcm,2026-03-31T18:00:00Z,SAFE,compute-bound,0,,"Downgraded to MEDIUM. da.stack without rechunk is scheduling overhead, not OOM risk." 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." hydro,2026-05-01,RISKY,memory-bound,0,1416,"Fixed-in-tree 2026-05-01: hand_mfd._hand_mfd_dask now assembles via da.map_blocks instead of eager da.block of pre-computed tiles (matches hand_dinf pattern). Remaining MEDIUM: sink_d8 CCL fully materializes labels (inherently global), flow_accumulation_mfd frac_bdry held in driver dict instead of memmap-backed BoundaryStore. D8 iterative paths (flow_accum/fill/watershed/basin/stream_*) use serial-tile sweep with memmap-backed boundary store -- per-tile RAM bounded but driver iterates O(diameter) times. flow_direction_*, flow_path/snap_pour_point/twi/hand_d8/hand_dinf are SAFE." From f2fc74b379e0fdeaefcb661cecbbef960bc6d13a Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Sun, 10 May 2026 22:05:52 -0700 Subject: [PATCH 3/6] Record geotiff accuracy sweep pass 11 (2026-05-10): clean Audited the one geotiff commit added since pass 10 -- #1559 (PR 1548), which centralises attrs population across the four read backends. The change is a pure metadata refactor; no data-path arithmetic was touched. Windowed-origin math and cross-backend attrs/data parity verified. --- .claude/sweep-accuracy-state.csv | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.claude/sweep-accuracy-state.csv b/.claude/sweep-accuracy-state.csv index a8872d7fa..fbf2f0804 100644 --- a/.claude/sweep-accuracy-state.csv +++ b/.claude/sweep-accuracy-state.csv @@ -12,7 +12,7 @@ emerging_hotspots,2026-04-30,,MEDIUM,2;3,MEDIUM: threshold_90 uses int() (trunca fire,2026-04-30,,,,All ops per-pixel (no accumulation/stencil/projected distance). NaN handled via x!=x; CUDA bounds use strict <; rdnbr and ros divisions guarded; CPU/GPU/dask paths algorithmically identical. No accuracy issues found. flood,2026-04-30,,MEDIUM,2;5,"MEDIUM (not fixed): dask backend preserves float32 input dtype while numpy promotes to float64 in flood_depth and curve_number_runoff; DataArray inputs for curve_number, mannings_n bypass scalar > 0 (and CN <= 100) range validation, silently producing NaN/garbage." focal,2026-03-30T13:00:00Z,1092,,, -geotiff,2026-05-10,,,,"Pass 10 (2026-05-10): CLEAN. Audited 5 recent commits: #1558 drop-defensive-copies (frombuffer path still .copy()s before in-place predictor decode at _reader.py:778), #1556 fp-predictor ngjit (writer pre-ravels so 1-D slice arg is correct, float32/64 LE+BE bit-exact), #1552 batched D2H (OOM guard fires before cupy.concatenate, host_buf offsets correct), #1551 parallel-decode gate (>= vs > sends 256x256 default to parallel path, no value diff confirmed via partial-tile parity), #1549 nvjpeg constants (gray + RGB GPU JPEG decode pixel-identical to Pillow CPU, max diff = 0). Cross-backend parity re-verified clean: numpy/dask+numpy/cupy/dask+cupy equal .data/.dtype/.coords/nodata/NaN-mask on deflate+predictor3+nodata; orientations 1-8 numpy==GPU; partial edge tiles 100x150, 257x383, 512x257 numpy==GPU==dask; predictor2 LE/BE round-trip uint8/int16/uint16/int32/uint32 pass; predictor3 LE/BE float32/64 pass. Deferred LOW (pre-existing, not opened): float16 (bps=16, SampleFormat=3) absent from tiff_dtype_to_numpy map - writer never emits, asymmetric but unreachable. | Pass 9 (2026-05-09): TWO HIGH fixed -- (a) PR #1539 closes #1537: TIFF Orientation tag 2/3/4 (mirror flips) on georeferenced files left y/x coords computed from the un-flipped transform, so xarray label lookups returned the wrong pixel even though _apply_orientation flipped the buffer. PR #1521 only updated the transform for the 5-8 axis-swap branch. Fix updates origin and pixel-scale signs along whichever axes were flipped, for both PixelIsArea (origin shifts by N*step) and PixelIsPoint (shifts by (N-1)*step). 10 new tests in test_orientation.py. (b) PR #1546 closes #1540: read_geotiff_gpu ignored Orientation tag completely; CPU correctly applied 2-8 (PR #1521) but GPU returned the raw stored buffer. Cross-backend disagreement on every non-default orientation. Fix adds _apply_orientation_gpu (cupy slicing mirror of the CPU helper) and _apply_orientation_geo_info, threads them into the tiled GPU pipeline, reuses CPU-fallback geo_info for the stripped path to avoid double-applying. 28 new tests in test_orientation_gpu.py (every orientation, single-band tiled, single-band stripped, 3-band tiled, mirror-flip sel-fidelity, default no-tag passthrough). Re-confirmed clean: HTTP coalesce_ranges with overlapping ranges and zero-length ranges, parallel streaming write thread-safety (each tile gets independent buffer via copy or padded zeros), planar=2 + chunky GPU LERC mask propagation matches CPU, IFD chain cap MAX_IFDS=256, max_z_error round-trip on tiled write, _resolve_masked_fill float vs integer dtype semantics. Deferred LOW: per-sample LERC mask (3D mask (h,w,samples)) collapsed to per-pixel ""any sample invalid"" on GPU while CPU honours per-sample; LERC implementations rarely emit 3D masks (verified: lerc.encode with 2D mask on 3-band returns 2D mask). Documented planar=2 + LERC + GPU silently drops mask (rare in practice, source comment acknowledges). | Pass 8 (2026-05-07): HIGH fixed in fix-jpeg-tiff-disable -- to_geotiff(compression='jpeg') wrote files that no external reader can decode. The writer tags compression=7 (new-style JPEG) but emits a self-contained JFIF stream per tile/strip and never writes the JPEGTables tag (347) that the TIFF spec requires for that codec. libtiff/GDAL/rasterio all reject the file with TIFFReadEncodedStrip() failed; our reader round-trips because Pillow decodes the standalone JFIF, hiding the break. Pass-4 notes flagged the read side of the same JPEGTables gap and deferred it; pass-8 covers the write side. Fix: reject compression='jpeg' at the to_geotiff entry with a clear ValueError pointing at deflate/zstd/lzw. The internal _writer.write is untouched so the existing self-decoding tests still cover the codec; re-enabling the public path needs a JPEGTables-aware encoder. PR diffs reviewed but not merged: #1512 (BytesIO source) and #1513 (LERC max_z_error) -- both look correct; #1512 file-like read path goes through read_all() once so the per-call BytesIOSource lock is theoretical, and #1513 forwards max_z_error through every overview/tile/strip/streaming path including _write_vrt_tiled and _compress_block. No regressions found in either open PR. Other surfaces audited clean: predictor=3 with float16 (writer auto-promotes to float32 on both eager and streaming paths, value-exact round-trip); planar=2 multi-tile read uses band_idx*tiles_per_band offset so no cross-contamination between planes; _header.py multi-byte tag parsing uses bo (byte_order) consistently; Pillow YCbCr-vs-tagged-RGB photometric mismatch becomes moot once JPEG is disabled. Deferred (LOW/MEDIUM, not filed): JPEG2000 writer accepts arbitrary dtype with no validation (rare codec, narrow risk); float16 dtype not in tiff_dtype_to_numpy decode map (writer never emits it - asymmetric but unreachable); Orientation tag (274) still ignored on read (pass-4 deferral). | Pass 7 (2026-05-07): HIGH fixed in fix-mmap-cache-refcount-after-replace -- _MmapCache.release() looked up the cache entry by realpath, so a holder that acquired the OLD mmap before an os.replace and released it AFTER another caller had acquired the post-replace entry would decrement the new holder's refcount. Subsequent eviction (cache full, or another acquire) closed the still-in-use mmap, breaking reads with 'mmap closed or invalid'. Real exposure: any concurrent reader/writer pattern where to_geotiff replaces a file that another reader had just opened via open_geotiff with chunks= or via _FileSource. PR #1506 added stale-replacement detection but did not fix the refcount confusion across the pop. Fix: acquire returns an opaque entry token; release takes the token and decrements that exact entry, regardless of cache state. Orphaned (popped) entries close their fh+mmap when their own refcount hits zero. _FileSource updated to pass the token. Regression test test_release_after_path_replacement_does_not_clobber_new_holder added. All 665 geotiff tests pass; GPU path verified. | Pass 6 (2026-05-07) PR #1507: BE pred2 numba TypingError. | Pass 5 (2026-05-06) PR #1506: mmap cache stale after file replace. | Pass 4 (2026-05-06) PR #1501: sparse COG tiles. | Pass 3 (2026-05-06) PR #1500: predictor=3 byte order. | Pass 2 (2026-05-05) PR #1498: predictor=2 sample-wise. | Pass 1 (2026-04-23) PR #1247. Re-confirmed clean over passes 2-7: items 2 (writer always emits LE TIFFs - hardcoded b'II'), 3 (RowsPerStrip default = height when missing), 4 (StripByteCounts missing raises clear ValueError), 5 (TileWidth without TileLength caught by 'tw <= 0 or th <= 0' check at _reader.py:688), 9 (read determinism on compressed+tiled+multiband), 11 (predictor=2 with awkward sample stride round-trips), 18 (compression_level=99 raises ValueError 'out of range for deflate (valid: 1-9)'), 21 (concurrent writes serialize correctly via mkstemp+os.replace), 24 (uint16 dtype preserved on numpy backend, dask honors chunks param), 26 (chunks rounds correctly with remainder chunk for non-tile-aligned). Deferred: item 8 (BytesIO/file-like sources are not supported, source.lower() error) - documented as 'str' parameter, not a bug; item 19 (LERC max_z_error not user-exposed by to_geotiff) - missing feature, not a bug." +geotiff,2026-05-10,,,,"Pass 11 (2026-05-10): CLEAN. Audited the one additional commit since pass 10 -- #1559 (PR 1548, Centralise GeoTIFF attrs population across all read backends). Refactor extracts _populate_attrs_from_geo_info helper and routes eager numpy, dask, GPU stripped, GPU tiled read paths through it; before the fix dask only emitted crs/transform/raster_type/nodata while numpy emitted the full attrs set including x/y_resolution, resolution_unit, image_description, extra_samples, GDAL metadata, and the CRS-description fields. No data-path arithmetic touched; only attrs dict population. Windowed origin math (origin_x + c0*pixel_width, origin_y + r0*pixel_height) verified to produce -98.0 / 48.75 origin for window=(10,20,50,70) on a (0.1,-0.125) pixel-size raster, with PixelIsArea half-pixel offset preserved on coord lookups (-97.95, 48.6875). Cross-backend attrs parity re-verified: numpy/dask/cupy all emit identical key set on deflate+predictor3+nodata round-trip (crs, crs_wkt, nodata, transform, x_resolution, y_resolution). Data bit-parity re-verified across numpy/dask/cupy on same payload (np.array_equal with equal_nan=True). test_attrs_parity_1548.py (5 tests), test_reader.py/test_writer.py/test_dask_cupy_combined.py (25 tests), GPU orientation/predictor2-BE/LERC-mask/nodata/byteswap suites (65 tests) all green. No accuracy or backend-divergence findings. | Pass 10 (2026-05-10): CLEAN. Audited 5 recent commits: #1558 drop-defensive-copies (frombuffer path still .copy()s before in-place predictor decode at _reader.py:778), #1556 fp-predictor ngjit (writer pre-ravels so 1-D slice arg is correct, float32/64 LE+BE bit-exact), #1552 batched D2H (OOM guard fires before cupy.concatenate, host_buf offsets correct), #1551 parallel-decode gate (>= vs > sends 256x256 default to parallel path, no value diff confirmed via partial-tile parity), #1549 nvjpeg constants (gray + RGB GPU JPEG decode pixel-identical to Pillow CPU, max diff = 0). Cross-backend parity re-verified clean: numpy/dask+numpy/cupy/dask+cupy equal .data/.dtype/.coords/nodata/NaN-mask on deflate+predictor3+nodata; orientations 1-8 numpy==GPU; partial edge tiles 100x150, 257x383, 512x257 numpy==GPU==dask; predictor2 LE/BE round-trip uint8/int16/uint16/int32/uint32 pass; predictor3 LE/BE float32/64 pass. Deferred LOW (pre-existing, not opened): float16 (bps=16, SampleFormat=3) absent from tiff_dtype_to_numpy map - writer never emits, asymmetric but unreachable. | Pass 9 (2026-05-09): TWO HIGH fixed -- (a) PR #1539 closes #1537: TIFF Orientation tag 2/3/4 (mirror flips) on georeferenced files left y/x coords computed from the un-flipped transform, so xarray label lookups returned the wrong pixel even though _apply_orientation flipped the buffer. PR #1521 only updated the transform for the 5-8 axis-swap branch. Fix updates origin and pixel-scale signs along whichever axes were flipped, for both PixelIsArea (origin shifts by N*step) and PixelIsPoint (shifts by (N-1)*step). 10 new tests in test_orientation.py. (b) PR #1546 closes #1540: read_geotiff_gpu ignored Orientation tag completely; CPU correctly applied 2-8 (PR #1521) but GPU returned the raw stored buffer. Cross-backend disagreement on every non-default orientation. Fix adds _apply_orientation_gpu (cupy slicing mirror of the CPU helper) and _apply_orientation_geo_info, threads them into the tiled GPU pipeline, reuses CPU-fallback geo_info for the stripped path to avoid double-applying. 28 new tests in test_orientation_gpu.py (every orientation, single-band tiled, single-band stripped, 3-band tiled, mirror-flip sel-fidelity, default no-tag passthrough). Re-confirmed clean: HTTP coalesce_ranges with overlapping ranges and zero-length ranges, parallel streaming write thread-safety (each tile gets independent buffer via copy or padded zeros), planar=2 + chunky GPU LERC mask propagation matches CPU, IFD chain cap MAX_IFDS=256, max_z_error round-trip on tiled write, _resolve_masked_fill float vs integer dtype semantics. Deferred LOW: per-sample LERC mask (3D mask (h,w,samples)) collapsed to per-pixel ""any sample invalid"" on GPU while CPU honours per-sample; LERC implementations rarely emit 3D masks (verified: lerc.encode with 2D mask on 3-band returns 2D mask). Documented planar=2 + LERC + GPU silently drops mask (rare in practice, source comment acknowledges). | Pass 8 (2026-05-07): HIGH fixed in fix-jpeg-tiff-disable -- to_geotiff(compression='jpeg') wrote files that no external reader can decode. The writer tags compression=7 (new-style JPEG) but emits a self-contained JFIF stream per tile/strip and never writes the JPEGTables tag (347) that the TIFF spec requires for that codec. libtiff/GDAL/rasterio all reject the file with TIFFReadEncodedStrip() failed; our reader round-trips because Pillow decodes the standalone JFIF, hiding the break. Pass-4 notes flagged the read side of the same JPEGTables gap and deferred it; pass-8 covers the write side. Fix: reject compression='jpeg' at the to_geotiff entry with a clear ValueError pointing at deflate/zstd/lzw. The internal _writer.write is untouched so the existing self-decoding tests still cover the codec; re-enabling the public path needs a JPEGTables-aware encoder. PR diffs reviewed but not merged: #1512 (BytesIO source) and #1513 (LERC max_z_error) -- both look correct; #1512 file-like read path goes through read_all() once so the per-call BytesIOSource lock is theoretical, and #1513 forwards max_z_error through every overview/tile/strip/streaming path including _write_vrt_tiled and _compress_block. No regressions found in either open PR. Other surfaces audited clean: predictor=3 with float16 (writer auto-promotes to float32 on both eager and streaming paths, value-exact round-trip); planar=2 multi-tile read uses band_idx*tiles_per_band offset so no cross-contamination between planes; _header.py multi-byte tag parsing uses bo (byte_order) consistently; Pillow YCbCr-vs-tagged-RGB photometric mismatch becomes moot once JPEG is disabled. Deferred (LOW/MEDIUM, not filed): JPEG2000 writer accepts arbitrary dtype with no validation (rare codec, narrow risk); float16 dtype not in tiff_dtype_to_numpy decode map (writer never emits it - asymmetric but unreachable); Orientation tag (274) still ignored on read (pass-4 deferral). | Pass 7 (2026-05-07): HIGH fixed in fix-mmap-cache-refcount-after-replace -- _MmapCache.release() looked up the cache entry by realpath, so a holder that acquired the OLD mmap before an os.replace and released it AFTER another caller had acquired the post-replace entry would decrement the new holder's refcount. Subsequent eviction (cache full, or another acquire) closed the still-in-use mmap, breaking reads with 'mmap closed or invalid'. Real exposure: any concurrent reader/writer pattern where to_geotiff replaces a file that another reader had just opened via open_geotiff with chunks= or via _FileSource. PR #1506 added stale-replacement detection but did not fix the refcount confusion across the pop. Fix: acquire returns an opaque entry token; release takes the token and decrements that exact entry, regardless of cache state. Orphaned (popped) entries close their fh+mmap when their own refcount hits zero. _FileSource updated to pass the token. Regression test test_release_after_path_replacement_does_not_clobber_new_holder added. All 665 geotiff tests pass; GPU path verified. | Pass 6 (2026-05-07) PR #1507: BE pred2 numba TypingError. | Pass 5 (2026-05-06) PR #1506: mmap cache stale after file replace. | Pass 4 (2026-05-06) PR #1501: sparse COG tiles. | Pass 3 (2026-05-06) PR #1500: predictor=3 byte order. | Pass 2 (2026-05-05) PR #1498: predictor=2 sample-wise. | Pass 1 (2026-04-23) PR #1247. Re-confirmed clean over passes 2-7: items 2 (writer always emits LE TIFFs - hardcoded b'II'), 3 (RowsPerStrip default = height when missing), 4 (StripByteCounts missing raises clear ValueError), 5 (TileWidth without TileLength caught by 'tw <= 0 or th <= 0' check at _reader.py:688), 9 (read determinism on compressed+tiled+multiband), 11 (predictor=2 with awkward sample stride round-trips), 18 (compression_level=99 raises ValueError 'out of range for deflate (valid: 1-9)'), 21 (concurrent writes serialize correctly via mkstemp+os.replace), 24 (uint16 dtype preserved on numpy backend, dask honors chunks param), 26 (chunks rounds correctly with remainder chunk for non-tile-aligned). Deferred: item 8 (BytesIO/file-like sources are not supported, source.lower() error) - documented as 'str' parameter, not a bug; item 19 (LERC max_z_error not user-exposed by to_geotiff) - missing feature, not a bug." glcm,2026-05-01,1408,HIGH,2,"angle=None averaged NaN as 0, masking no-valid-pairs as zero texture; fixed via nanmean-style averaging" hillshade,2026-04-10T12:00:00Z,,,,"Horn's method correct. All backends consistent. NaN propagation correct. float32 adequate for [0,1] output." hydro,2026-04-30,,LOW,1,Only LOW: twi log(0)=-inf if fa=0 (out-of-contract); MFD weighted sum no Kahan (negligible). No CRIT/HIGH issues. From 8cecca8876e016f841d0f4540229717ddb9d2232 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Sun, 10 May 2026 22:07:33 -0700 Subject: [PATCH 4/6] Record geotiff api-consistency sweep pass: 3 issues filed Files issues 1560 (gpu kwarg type drift, HIGH), 1561 (backend signature parity gap, HIGH), and 1562 (writer docstring drift, MEDIUM). The prior orphan-__all__ HIGH from issue 1541 already landed in PR 1544. Categories 1, 3, 5 found. --- .claude/sweep-api-consistency-state.csv | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.claude/sweep-api-consistency-state.csv b/.claude/sweep-api-consistency-state.csv index fe34f9646..f72d1c1bf 100644 --- a/.claude/sweep-api-consistency-state.csv +++ b/.claude/sweep-api-consistency-state.csv @@ -1,2 +1,2 @@ module,last_inspected,issue,severity_max,categories_found,notes -geotiff,2026-05-09,1541,HIGH,1;2;3;4;5,"PR adds 4 prod fns to __all__, real DeprecationWarning on plot_geotiff. Other findings (gpu kwarg type drift, sibling param asymmetry, write/dask_data rename, chunks default drift, write docstring drift) listed in issue #1541 for follow-up." +geotiff,2026-05-10,1560;1561;1562,HIGH,1;3;5,"Filed gpu kwarg type drift (1560 HIGH), backend signature parity gap (1561 HIGH), writer docstring drift (1562 MEDIUM). Prior orphan-__all__ HIGH fixed in PR 1544. Cat 1 MEDIUM (write/dask_data internal asymmetry) and Cat 4 MEDIUM (chunks default) reviewed and judged defensible per private-helper / different-semantics scope." From 415b62c0dfa53014e2e5bb1923c0ac443714e3b4 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Sun, 10 May 2026 22:09:59 -0700 Subject: [PATCH 5/6] Close geotiff test-coverage gaps deferred from 2026-05-09 sweep The previous test-coverage sweep flagged four HIGH/MEDIUM gaps as deferred. Add tests pinning each path so a regression cannot ship undetected: - dask + planar multiband (planar=1 contig, planar=2 separate; stripped + tiled; bands 3 and 4; uint8 + uint16) - dask + overview_level (levels 0, 1, 2, None over a 3-IFD pyramid) - predictor=2 with int8 samples (little- and big-endian; stripped and tiled; CPU and GPU read paths) - dask chunk size misaligned with on-disk tile size (chunks smaller than tile, larger than tile, doubly-unaligned corner cell) All 29 new tests pass locally. GPU tests use the existing cupy/CUDA-availability guard. --- .claude/sweep-test-coverage-state.csv | 4 +- .../test_dask_chunk_tile_misalignment.py | 115 +++++++++++++++ .../geotiff/tests/test_dask_overview_level.py | 115 +++++++++++++++ .../tests/test_dask_planar_multiband.py | 121 +++++++++++++++ .../geotiff/tests/test_predictor2_int8.py | 138 ++++++++++++++++++ 5 files changed, 491 insertions(+), 2 deletions(-) create mode 100644 xrspatial/geotiff/tests/test_dask_chunk_tile_misalignment.py create mode 100644 xrspatial/geotiff/tests/test_dask_overview_level.py create mode 100644 xrspatial/geotiff/tests/test_dask_planar_multiband.py create mode 100644 xrspatial/geotiff/tests/test_predictor2_int8.py diff --git a/.claude/sweep-test-coverage-state.csv b/.claude/sweep-test-coverage-state.csv index ef5bc45fa..510c4d364 100644 --- a/.claude/sweep-test-coverage-state.csv +++ b/.claude/sweep-test-coverage-state.csv @@ -1,2 +1,2 @@ -module,last_inspected,issue,severity_max,categories_found,notes -geotiff,2026-05-09,1543,MEDIUM,1,"Filed #1543 / PR adds dask+cupy combined backend tests; remaining HIGH gaps (realistic dask chunk-boundary, planar dask, orientation coord-flip, pred=2 int8) deferred" +module,last_inspected,issue,severity_max,categories_found,notes +geotiff,2026-05-10,,HIGH,1;4,"Closed remaining HIGH gaps: dask + planar multiband, dask + overview_level, predictor=2 int8 (CPU+GPU), dask chunk-vs-tile misalignment" diff --git a/xrspatial/geotiff/tests/test_dask_chunk_tile_misalignment.py b/xrspatial/geotiff/tests/test_dask_chunk_tile_misalignment.py new file mode 100644 index 000000000..cc5808b09 --- /dev/null +++ b/xrspatial/geotiff/tests/test_dask_chunk_tile_misalignment.py @@ -0,0 +1,115 @@ +"""``read_geotiff_dask`` chunk boundaries misaligned with TIFF tile size. + +``read_geotiff_dask`` builds chunks of size ``chunks`` (default 512) +regardless of the underlying file's ``TileWidth``/``TileLength`` tags. +When the requested chunk size does not align with the on-disk tile +grid, the per-window reader must re-tile its decoded tile buffer into +the requested window before returning it to dask. Existing dask tests +all use chunk sizes that line up with the tile boundary; this module +covers the misaligned case so a regression in the windowed re-tile +path (off-by-one cropping, wrong row stride at a tile-spanning chunk, +band-axis misalignment) does not ship undetected. + +Three flavours of misalignment are exercised: + + * Chunk smaller than tile (e.g. ``chunks=11`` on a 16-tile file): a + single tile must be diced into multiple chunks. + * Chunk larger than tile and not a multiple (e.g. ``chunks=23`` on + a 16-tile file): a single chunk must stitch fragments from + multiple tiles. + * Final chunk that crops both axes simultaneously (image size not a + multiple of chunk size, and chunk size not a multiple of tile + size). Catches the corner cell where every boundary is partial. +""" +from __future__ import annotations + +import numpy as np +import pytest + + +tifffile = pytest.importorskip("tifffile") +dask_array = pytest.importorskip("dask.array") + + +def _write_tiled(path: str, arr: np.ndarray, tile: int = 16) -> None: + """Write *arr* as a tiled TIFF with the requested tile size.""" + tifffile.imwrite(str(path), arr, tile=(tile, tile), + photometric="minisblack", compression="deflate") + + +@pytest.fixture(scope="module") +def _arr_64x96(): + """Deterministic 64x96 uint16 raster reused across chunk-size cases.""" + rng = np.random.RandomState(0xC4AE) + return rng.randint(0, 60_000, size=(64, 96), dtype=np.uint16) + + +def test_chunk_smaller_than_tile(tmp_path, _arr_64x96): + """``chunks=11`` on a 16x16-tile file: tile is subdivided across chunks. + + With image 64x96 and chunks=11 the dask layout is ceil(64/11)=6 row + blocks and ceil(96/11)=9 column blocks (54 chunks total). Each + chunk straddles a tile boundary -- if the window-to-tile mapping + is off by one row or column, the computed value will differ from + the source. + """ + from xrspatial.geotiff import read_geotiff_dask + + path = tmp_path / "tiled_misalign_small.tif" + _write_tiled(path, _arr_64x96, tile=16) + + da_arr = read_geotiff_dask(str(path), chunks=11) + assert isinstance(da_arr.data, dask_array.Array) + # 11 < 16: every tile is dispersed across at least 2 chunks. + assert da_arr.data.chunksize[:2] == (11, 11) + np.testing.assert_array_equal(da_arr.compute().values, _arr_64x96) + + +def test_chunk_larger_than_tile_nonmultiple(tmp_path, _arr_64x96): + """``chunks=23`` on a 16x16-tile file: each chunk stitches partial tiles. + + 23 % 16 == 7, so every chunk pulls bytes from a partial tile on at + least one side. If the reader rounds the requested window down to + the nearest tile boundary, the chunk shape comes out wrong; if it + rounds up, the values shift. + """ + from xrspatial.geotiff import read_geotiff_dask + + path = tmp_path / "tiled_misalign_large.tif" + _write_tiled(path, _arr_64x96, tile=16) + + da_arr = read_geotiff_dask(str(path), chunks=23) + assert isinstance(da_arr.data, dask_array.Array) + assert da_arr.data.chunksize[:2] == (23, 23) + np.testing.assert_array_equal(da_arr.compute().values, _arr_64x96) + + +def test_chunk_tuple_doubly_unaligned(tmp_path): + """Image not a multiple of chunk, chunk not a multiple of tile. + + Image 50x70, tile 16x16, chunks (17, 19). The final row chunk and + final column chunk both crop, and neither chunk dimension is + aligned with the tile grid. This is the corner-cell case. + """ + from xrspatial.geotiff import read_geotiff_dask + + rng = np.random.RandomState(0xDCED) + arr = rng.randint(0, 256, size=(50, 70), dtype=np.uint8) + + path = tmp_path / "tiled_corner_misalign.tif" + _write_tiled(path, arr, tile=16) + + da_arr = read_geotiff_dask(str(path), chunks=(17, 19)) + assert da_arr.shape == (50, 70) + # Last block in each axis is the trimmed remainder. + block_h = da_arr.data.chunks[0] + block_w = da_arr.data.chunks[1] + assert block_h == (17, 17, 16), ( + f"row chunks should be 17,17,16 (50-pixel image, chunks=17), " + f"got {block_h}" + ) + assert block_w == (19, 19, 19, 13), ( + f"col chunks should be 19,19,19,13 (70-pixel image, chunks=19), " + f"got {block_w}" + ) + np.testing.assert_array_equal(da_arr.compute().values, arr) diff --git a/xrspatial/geotiff/tests/test_dask_overview_level.py b/xrspatial/geotiff/tests/test_dask_overview_level.py new file mode 100644 index 000000000..103bc2c33 --- /dev/null +++ b/xrspatial/geotiff/tests/test_dask_overview_level.py @@ -0,0 +1,115 @@ +"""Dask read of a TIFF with overview levels (COG pyramid). + +``read_geotiff_dask`` accepts an ``overview_level`` kwarg that forwards +to ``_read_geo_info`` for IFD selection and to ``_delayed_read_window`` +for per-chunk decoding. Before this module landed, no test exercised +that combination, so a regression where the dask path silently ignored +the overview level (returning full-resolution chunks) or sampled the +wrong IFD would not be caught. + +This module writes a small COG-style file with two overview levels +(half- and quarter-resolution) and asserts that: + + * The returned ``DataArray`` shape matches the overview level's + dimensions, not the full-resolution dimensions. + * The per-chunk windowed reader pulls bytes from the correct IFD + (the computed values agree with a non-dask reference at the same + overview level). + * ``overview_level=None`` (default) still returns the full-resolution + image, so the new code path does not change default behaviour. +""" +from __future__ import annotations + +import numpy as np +import pytest + + +tifffile = pytest.importorskip("tifffile") +dask_array = pytest.importorskip("dask.array") + + +def _write_cog_with_overviews(path: str, data: np.ndarray) -> None: + """Write *data* as a tiled TIFF with two precomputed overview IFDs. + + Writes the primary IFD followed by half- and quarter-resolution + overview IFDs, each tagged ``subfiletype=1`` so the reader treats + them as a pyramid (matching how ``_write_normal_cog`` in + ``test_overview_filter.py`` builds COG fixtures). This mirrors what + GDAL's ``gdaladdo`` emits. + """ + half = data[::2, ::2] + quart = data[::4, ::4] + with tifffile.TiffWriter(path) as tw: + tw.write(data, tile=(32, 32), photometric="minisblack") + tw.write(half, tile=(32, 32), photometric="minisblack", + subfiletype=1) + tw.write(quart, tile=(32, 32), photometric="minisblack", + subfiletype=1) + + +def test_dask_overview_level_zero_matches_full_res(tmp_path): + """``overview_level=0`` returns full resolution (the base IFD).""" + from xrspatial.geotiff import read_geotiff_dask + + rng = np.random.RandomState(0xD0E) + arr = rng.randint(0, 256, size=(128, 192), dtype=np.uint8) + path = str(tmp_path / "cog_dask_ov.tif") + _write_cog_with_overviews(path, arr) + + da_arr = read_geotiff_dask(path, chunks=32, overview_level=0) + assert da_arr.shape == arr.shape + np.testing.assert_array_equal(da_arr.compute().values, arr) + + +def test_dask_overview_level_one_returns_half_res(tmp_path): + """``overview_level=1`` materialises the half-resolution overview.""" + from xrspatial.geotiff import read_geotiff_dask + from xrspatial.geotiff import open_geotiff + + rng = np.random.RandomState(0xD0E) + arr = rng.randint(0, 256, size=(128, 192), dtype=np.uint8) + path = str(tmp_path / "cog_dask_ov1.tif") + _write_cog_with_overviews(path, arr) + + # Eager reference at the same overview level -- the dask path should + # pull the same bytes from the same IFD. + eager = open_geotiff(path, overview_level=1) + + da_arr = read_geotiff_dask(path, chunks=16, overview_level=1) + assert da_arr.shape == eager.shape, ( + f"dask returned {da_arr.shape} but eager returned {eager.shape} " + "at overview_level=1" + ) + assert isinstance(da_arr.data, dask_array.Array) + np.testing.assert_array_equal(da_arr.compute().values, eager.values) + + +def test_dask_overview_level_two_returns_quarter_res(tmp_path): + """``overview_level=2`` materialises the quarter-resolution overview.""" + from xrspatial.geotiff import read_geotiff_dask + from xrspatial.geotiff import open_geotiff + + rng = np.random.RandomState(0xD0E) + arr = rng.randint(0, 256, size=(128, 192), dtype=np.uint8) + path = str(tmp_path / "cog_dask_ov2.tif") + _write_cog_with_overviews(path, arr) + + eager = open_geotiff(path, overview_level=2) + + da_arr = read_geotiff_dask(path, chunks=8, overview_level=2) + assert da_arr.shape == eager.shape + np.testing.assert_array_equal(da_arr.compute().values, eager.values) + + +def test_dask_overview_level_none_returns_full_res(tmp_path): + """``overview_level=None`` keeps default behaviour: full resolution.""" + from xrspatial.geotiff import read_geotiff_dask + + rng = np.random.RandomState(0xD0E) + arr = rng.randint(0, 256, size=(128, 192), dtype=np.uint8) + path = str(tmp_path / "cog_dask_ov_none.tif") + _write_cog_with_overviews(path, arr) + + da_arr = read_geotiff_dask(path, chunks=32, overview_level=None) + assert da_arr.shape == arr.shape + np.testing.assert_array_equal(da_arr.compute().values, arr) diff --git a/xrspatial/geotiff/tests/test_dask_planar_multiband.py b/xrspatial/geotiff/tests/test_dask_planar_multiband.py new file mode 100644 index 000000000..197aa7434 --- /dev/null +++ b/xrspatial/geotiff/tests/test_dask_planar_multiband.py @@ -0,0 +1,121 @@ +"""Dask read of multi-band planar TIFF files. + +``read_geotiff_dask`` advertises multi-band support through the ``n_bands`` +branch in ``read_geotiff_dask`` -- when ``samples_per_pixel > 1`` the +returned ``DataArray`` is shaped ``(y, x, band)``. Until this module +landed, no test exercised that branch through the dask reader, so a +regression in the underlying ``read_to_array(window=...)`` call for +``PlanarConfiguration=2`` (planar=separate, the COG-friendly layout) +would ship undetected. + +This module pins: + + * ``read_geotiff_dask`` returns the expected ``(y, x, band)`` shape + and dtype for both planar=1 (contig, chunky) and planar=2 + (separate) source files. + * The computed values match the original numpy buffer pixel-for-pixel + after the lazy dask graph is materialised. + * Chunks tuple (row_chunk, col_chunk) is honoured on the y/x axes + while the band axis stays a single contiguous chunk. + +Both stripped and tiled file layouts are covered because the +``_decode_strip_or_tile`` path branches on layout and the planar +handling differs in each branch. +""" +from __future__ import annotations + +import numpy as np +import pytest + + +tifffile = pytest.importorskip("tifffile") +dask_array = pytest.importorskip("dask.array") + + +def _write_planar_tiff(path: str, data: np.ndarray, *, + planar: str, tiled: bool) -> None: + """Write *data* shaped ``(bands, height, width)`` with chosen layout. + + tifffile expects ``(bands, h, w)`` for ``planarconfig='separate'`` and + ``(h, w, bands)`` for ``planarconfig='contig'``. This helper centralises + the transpose so the test bodies stay focused on the assertion. + """ + kwargs: dict = {"photometric": "minisblack"} + if data.shape[0] == 3: + kwargs["photometric"] = "rgb" + if tiled: + kwargs["tile"] = (32, 32) + if planar == "separate": + kwargs["planarconfig"] = "separate" + tifffile.imwrite(path, data, **kwargs) + elif planar == "contig": + kwargs["planarconfig"] = "contig" + tifffile.imwrite(path, np.transpose(data, (1, 2, 0)), **kwargs) + else: + raise ValueError(f"unknown planar={planar!r}") + + +def _make_data(bands: int, height: int, width: int, dtype) -> np.ndarray: + rng = np.random.RandomState(0xD45C + bands * 100 + height) + info = np.iinfo(dtype) + high = min(int(info.max), 60_000) + 1 + return rng.randint(0, high, size=(bands, height, width)).astype(dtype) + + +@pytest.mark.parametrize("planar", ["separate", "contig"]) +@pytest.mark.parametrize("tiled", [True, False]) +@pytest.mark.parametrize("bands", [3, 4]) +@pytest.mark.parametrize("dtype", [np.uint8, np.uint16]) +def test_dask_planar_multiband_matches_numpy( + tmp_path, planar, tiled, bands, dtype +): + """``read_geotiff_dask`` returns ``(y, x, band)`` matching the source.""" + from xrspatial.geotiff import read_geotiff_dask + + height, width = 96, 128 + data = _make_data(bands, height, width, dtype) + # On disk the file stores ``(bands, h, w)`` but the reader returns + # the xarray convention ``(y, x, band)``. + expected = np.transpose(data, (1, 2, 0)) + + path = str(tmp_path + / f"dask_planar_{planar}_{'tile' if tiled else 'strip'}_" + f"b{bands}_{np.dtype(dtype).name}.tif") + _write_planar_tiff(path, data, planar=planar, tiled=tiled) + + da_arr = read_geotiff_dask(path, chunks=32) + + assert isinstance(da_arr.data, dask_array.Array), ( + f"expected dask Array, got {type(da_arr.data).__name__}" + ) + assert da_arr.shape == (height, width, bands), ( + f"shape mismatch: {da_arr.shape} vs {(height, width, bands)}" + ) + assert da_arr.dtype == np.dtype(dtype) + assert list(da_arr.dims) == ["y", "x", "band"] + + materialised = da_arr.compute().values + np.testing.assert_array_equal(materialised, expected) + + +def test_dask_planar_separate_chunks_tuple(tmp_path): + """Tuple chunks ``(ch_h, ch_w)`` honoured; band axis stays single chunk.""" + from xrspatial.geotiff import read_geotiff_dask + + bands, height, width = 3, 80, 120 + data = _make_data(bands, height, width, np.uint8) + expected = np.transpose(data, (1, 2, 0)) + + path = str(tmp_path / "dask_planar_chunktuple.tif") + _write_planar_tiff(path, data, planar="separate", tiled=True) + + da_arr = read_geotiff_dask(path, chunks=(40, 60)) + + # ``read_geotiff_dask`` builds row-major chunks of (ch_h, ch_w, n_bands). + # With height=80, width=120, chunks=(40, 60) the expected layout is + # 2 row blocks x 2 col blocks x 1 band block. + assert da_arr.data.chunksize[:2] == (40, 60) + # The band axis is concatenated as one block (n_bands shape). + assert da_arr.data.chunksize[2] == bands + + np.testing.assert_array_equal(da_arr.compute().values, expected) diff --git a/xrspatial/geotiff/tests/test_predictor2_int8.py b/xrspatial/geotiff/tests/test_predictor2_int8.py new file mode 100644 index 000000000..05dd10184 --- /dev/null +++ b/xrspatial/geotiff/tests/test_predictor2_int8.py @@ -0,0 +1,138 @@ +"""Predictor=2 with signed 8-bit (``int8``) samples. + +The predictor=2 (horizontal differencing) path historically tested +``uint8`` (byte-wise kernel) and 16/32-bit signed and unsigned integers +(multi-byte kernel). ``int8`` slipped through both buckets: + + * It has a single-byte stride so the reader takes the byte-wise + path, but the sample format flag is ``2`` (signed) -- a different + branch in the sample-format decoder than the well-trodden + ``uint8`` path (sample format ``1``). + * Signed 8-bit values exercise the integer-promotion semantics of + the cumulative-sum reconstruction. If the kernel accidentally + treated the byte stream as unsigned, the decode would produce + wraparound errors at every transition through zero. + +This module pins the read path: + + * Little-endian and big-endian predictor=2 ``int8`` files decode + back to the original signed values byte-for-byte. + * The GPU read path matches the CPU baseline. + * Stripped and tiled layouts both succeed (separate code paths). +""" +from __future__ import annotations + +import importlib.util + +import numpy as np +import pytest + + +tifffile = pytest.importorskip("tifffile") + + +def _gpu_available() -> bool: + if importlib.util.find_spec("cupy") is None: + return False + try: + import cupy + return bool(cupy.cuda.is_available()) + except Exception: + return False + + +_HAS_GPU = _gpu_available() +_gpu_only = pytest.mark.skipif( + not _HAS_GPU, + reason="cupy + CUDA required", +) + + +def _signed_int8_grid(h: int = 16, w: int = 16) -> np.ndarray: + """Deterministic ``int8`` raster spanning the full -128..127 range. + + The pattern is chosen to walk through the signed wraparound (going + from positive into negative and back) on every row, so a bug that + decodes the byte stream as unsigned would produce a different + cumulative sum than the signed reference. + """ + info = np.iinfo(np.int8) + rng = np.random.RandomState(0x1781) + return rng.randint(info.min, info.max + 1, + size=(h, w)).astype(np.int8) + + +@pytest.mark.parametrize("byteorder", ["<", ">"]) +def test_cpu_predictor2_int8_round_trip(tmp_path, byteorder): + """``int8`` predictor=2 round-trips through the CPU reader. + + Byte order is irrelevant for an 8-bit sample but the writer still + accepts the parameter; both forms must decode identically. + """ + from xrspatial.geotiff._reader import read_to_array + + arr = _signed_int8_grid() + path = tmp_path / f"int8_pred2_{'be' if byteorder == '>' else 'le'}.tif" + tifffile.imwrite(str(path), arr, byteorder=byteorder, predictor=2, + compression="deflate") + + out, _ = read_to_array(str(path)) + assert out.dtype == np.int8, f"dtype changed: {out.dtype}" + np.testing.assert_array_equal(out, arr) + + +def test_cpu_predictor2_int8_tiled(tmp_path): + """Tiled layout (separate ``_decode_tile`` path) handles int8 + pred=2.""" + from xrspatial.geotiff._reader import read_to_array + + arr = _signed_int8_grid(32, 48) + path = tmp_path / "int8_pred2_tiled.tif" + tifffile.imwrite(str(path), arr, predictor=2, compression="deflate", + tile=(16, 16)) + + out, _ = read_to_array(str(path)) + assert out.dtype == np.int8 + np.testing.assert_array_equal(out, arr) + + +@_gpu_only +def test_gpu_predictor2_int8_tiled_matches_cpu(tmp_path): + """The GPU tiled decoder matches the CPU baseline for int8 + pred=2.""" + import cupy + + from xrspatial.geotiff import read_geotiff_gpu + from xrspatial.geotiff._reader import read_to_array + + arr = _signed_int8_grid(32, 48) + path = tmp_path / "int8_pred2_tiled_gpu.tif" + tifffile.imwrite(str(path), arr, predictor=2, compression="deflate", + tile=(16, 16)) + + cpu, _ = read_to_array(str(path)) + np.testing.assert_array_equal(cpu, arr) + + gpu_da = read_geotiff_gpu(str(path)) + assert isinstance(gpu_da.data, cupy.ndarray) + assert gpu_da.data.dtype == np.int8 + np.testing.assert_array_equal(gpu_da.data.get(), cpu) + + +@_gpu_only +def test_gpu_predictor2_int8_stripped_matches_cpu(tmp_path): + """Stripped layout: GPU falls back to CPU; the result must still match.""" + import cupy + + from xrspatial.geotiff import read_geotiff_gpu + from xrspatial.geotiff._reader import read_to_array + + arr = _signed_int8_grid() + path = tmp_path / "int8_pred2_stripped_gpu.tif" + tifffile.imwrite(str(path), arr, predictor=2, compression="deflate") + + cpu, _ = read_to_array(str(path)) + np.testing.assert_array_equal(cpu, arr) + + gpu_da = read_geotiff_gpu(str(path)) + assert isinstance(gpu_da.data, cupy.ndarray) + assert gpu_da.data.dtype == np.int8 + np.testing.assert_array_equal(gpu_da.data.get(), cpu) From f48706a16534ea10657cfbd8a75eaebbab2aac54 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Sun, 10 May 2026 22:10:30 -0700 Subject: [PATCH 6/6] Record geotiff test-coverage sweep pass 2026-05-10: PR #1565 --- .claude/sweep-test-coverage-state.csv | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.claude/sweep-test-coverage-state.csv b/.claude/sweep-test-coverage-state.csv index 510c4d364..6ebeb5e2b 100644 --- a/.claude/sweep-test-coverage-state.csv +++ b/.claude/sweep-test-coverage-state.csv @@ -1,2 +1,2 @@ module,last_inspected,issue,severity_max,categories_found,notes -geotiff,2026-05-10,,HIGH,1;4,"Closed remaining HIGH gaps: dask + planar multiband, dask + overview_level, predictor=2 int8 (CPU+GPU), dask chunk-vs-tile misalignment" +geotiff,2026-05-10,1565,HIGH,1;3;4,"PR #1565 closes deferred HIGH gaps: dask+planar multiband, dask+overview_level, predictor=2 int8 CPU+GPU, dask chunk-vs-tile misalignment"