diff --git a/.claude/sweep-accuracy-state.csv b/.claude/sweep-accuracy-state.csv index 9aeed5f8..bc03591e 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-11,1599,HIGH,2;5,"Pass 13 (2026-05-11): HIGH fixed -- issue #1599. write_geotiff_gpu (and to_geotiff gpu=True) emitted raw NaN bytes for missing pixels even when nodata= was supplied, while the CPU writer substituted NaN with the sentinel before encoding. xrspatial-only round-trips were unaffected (the reader masks both NaN and the sentinel), but external readers (rasterio/GDAL/QGIS) that mask only on the GDAL_NODATA tag saw NaN pixels as valid data -- rasterio reported 100% valid pixels on a 25-NaN file vs CPU's 25-invalid report. Root cause: __init__.py lines 2579-2587 jumped from shape/dtype resolution straight to compression, missing the equivalent of the CPU writer's NaN-to-sentinel rewrite at to_geotiff line ~1156. Fix: cupy.isnan + masked write on a defensive copy of arr, gated on np_dtype.kind=='f' and not np.isnan(float(nodata)). Caller's CuPy buffer preserved (copy before mutate). 7 regression tests in test_gpu_writer_nan_sentinel_1599.py: substitution lands as sentinel, CPU/GPU byte-equivalent, caller buffer not mutated, no-NaN no-op, NaN sentinel skips substitution, rasterio sees identical invalid count on CPU/GPU, multiband 3D path. All other GPU writer tests still pass (50 passed across band-first, attrs, nodata, dask+cupy, writer, nodata aliases). | Pass 12 (2026-05-11): HIGH fixed -- issue #1581. Reading a uint TIFF with a negative GDAL_NODATA sentinel (e.g. uint16 + -9999) raised OverflowError on every backend because the nodata-mask code did arr.dtype.type(int(nodata)) with no range check. Three identical cast sites in __init__.py (numpy eager, _apply_nodata_mask_gpu, _delayed_read_window) plus _resolve_masked_fill and _sparse_fill_value in _reader.py. Fix: _int_nodata_in_range helper gates the cast; out-of-range sentinels are a no-op for value matching (the file can never contain that value), file dtype is preserved, attrs['nodata'] still surfaces the original sentinel so write round-trips keep the GDAL_NODATA tag intact. Matches rasterio behavior. 8 regression tests in test_nodata_out_of_range_1581.py cover the helper, both eager and dask read paths, in-range sentinel non-regression, and GPU helper (cupy-gated). | 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." +geotiff,2026-05-11,1611,HIGH,2;5,"Pass 14 (2026-05-11): HIGH fixed -- issue #1611. read_vrt(band=None) on a multi-band integer VRT with per-band tags only masks band 0's sentinel. __init__.py lines 2795-2809 in read_vrt apply vrt.bands[0].nodata to the full ndim==3 array; bands 1+ keep their integer sentinels as literal finite values (e.g. 65000 surfaces as 65000.0 after the dtype=float64 cast, not NaN). Float-VRT path masks per-band correctly in _vrt._read_data lines 296-297 + 347-351. PR #1602 fixed the single-band band=N case for issue #1598; the band=None multi-band case is the same class of bug. Repro: 2-band uint16 VRT with NoDataValue 65535 / 65000 returns r.values[1,1,1] == 65000.0 instead of NaN; r.values[1,1,0] is NaN (band 0 sentinel masked). Fix scope: in read_vrt, when band is None, iterate over vrt.bands and mask each arr[..., i] slice against its own (gated by the same _int_nodata_in_range guard PR #1583 introduced). Severity HIGH (Cat 2 NaN propagation + Cat 5 backend inconsistency: identical input semantics produce different masking outcomes based on dtype, with finite garbage values where NaN expected). Fix in PR #1612: walks vrt.bands when band is None and ndim==3, masks each arr[..., i] slice against its own via the refactored _sentinel_for_dtype helper (reuses PR #1583's range guard so out-of-range/non-finite/fractional sentinels are a no-op). attrs['nodata'] still carries band 0's sentinel for band=None reads (documented contract). 7 regression tests in test_vrt_multiband_int_nodata_1611.py: uint16 per-band, int32 negative, mixed presence, dtype preservation when no sentinel hit, out-of-range gating, band=N non-regression, attrs contract. 135 existing vrt/nodata geotiff tests still pass. | Pass 13 (2026-05-11): HIGH fixed -- issue #1599. write_geotiff_gpu (and to_geotiff gpu=True) emitted raw NaN bytes for missing pixels even when nodata= was supplied, while the CPU writer substituted NaN with the sentinel before encoding. xrspatial-only round-trips were unaffected (the reader masks both NaN and the sentinel), but external readers (rasterio/GDAL/QGIS) that mask only on the GDAL_NODATA tag saw NaN pixels as valid data -- rasterio reported 100% valid pixels on a 25-NaN file vs CPU's 25-invalid report. Root cause: __init__.py lines 2579-2587 jumped from shape/dtype resolution straight to compression, missing the equivalent of the CPU writer's NaN-to-sentinel rewrite at to_geotiff line ~1156. Fix: cupy.isnan + masked write on a defensive copy of arr, gated on np_dtype.kind=='f' and not np.isnan(float(nodata)). Caller's CuPy buffer preserved (copy before mutate). 7 regression tests in test_gpu_writer_nan_sentinel_1599.py: substitution lands as sentinel, CPU/GPU byte-equivalent, caller buffer not mutated, no-NaN no-op, NaN sentinel skips substitution, rasterio sees identical invalid count on CPU/GPU, multiband 3D path. All other GPU writer tests still pass (50 passed across band-first, attrs, nodata, dask+cupy, writer, nodata aliases). | Pass 12 (2026-05-11): HIGH fixed -- issue #1581. Reading a uint TIFF with a negative GDAL_NODATA sentinel (e.g. uint16 + -9999) raised OverflowError on every backend because the nodata-mask code did arr.dtype.type(int(nodata)) with no range check. Three identical cast sites in __init__.py (numpy eager, _apply_nodata_mask_gpu, _delayed_read_window) plus _resolve_masked_fill and _sparse_fill_value in _reader.py. Fix: _int_nodata_in_range helper gates the cast; out-of-range sentinels are a no-op for value matching (the file can never contain that value), file dtype is preserved, attrs['nodata'] still surfaces the original sentinel so write round-trips keep the GDAL_NODATA tag intact. Matches rasterio behavior. 8 regression tests in test_nodata_out_of_range_1581.py cover the helper, both eager and dask read paths, in-range sentinel non-regression, and GPU helper (cupy-gated). | 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. diff --git a/xrspatial/geotiff/__init__.py b/xrspatial/geotiff/__init__.py index e66e6dc9..dac78a0a 100644 --- a/xrspatial/geotiff/__init__.py +++ b/xrspatial/geotiff/__init__.py @@ -2792,21 +2792,65 @@ def read_vrt(source: str, *, dtype=None, window=None, # carried the literal sentinel value, breaking the convention that # downstream code follows (``attrs['nodata']`` is present iff the # array has already been NaN-masked). - if nodata is not None and arr.dtype.kind in ('u', 'i'): - # VRT nodata is parsed as float, so reject fractional, non-finite, or - # out-of-dtype-range values rather than truncate/wrap them into a - # sentinel that could mask real pixels. - nodata_f = float(nodata) - info = np.iinfo(arr.dtype) - if ( - np.isfinite(nodata_f) - and nodata_f.is_integer() - and info.min <= nodata_f <= info.max - ): - mask = arr == arr.dtype.type(int(nodata_f)) - if mask.any(): - arr = arr.astype(np.float64) - arr[mask] = np.nan + # + # For multi-band reads (``band is None`` and ``arr.ndim == 3``), each + # band can declare its own ````. The float-VRT path masks + # per-band inline in ``_vrt._read_data``; mirror that here by walking + # ``vrt.bands`` and masking each ``arr[..., i]`` slice against its own + # sentinel. Before this branch, only band 0's sentinel was applied and + # bands 1+ left their integer sentinels as literal finite values in + # the returned float64 array. See issue #1611. + def _sentinel_for_dtype(nodata_val, dtype): + """Return ``dtype``-cast sentinel for ``nodata_val`` or ``None`` + if the value can't be represented in ``dtype`` (non-integer + dtype, out-of-range, non-finite, or fractional). Mirrors the + gating PR #1583 added to other read paths via + ``_int_nodata_in_range``. + """ + if nodata_val is None or dtype.kind not in ('u', 'i'): + return None + try: + nodata_f = float(nodata_val) + except (TypeError, ValueError): + return None + info = np.iinfo(dtype) + if not (np.isfinite(nodata_f) and nodata_f.is_integer() + and info.min <= nodata_f <= info.max): + return None + return dtype.type(int(nodata_f)) + + if arr.dtype.kind in ('u', 'i'): + if arr.ndim == 3 and band is None and vrt.bands: + # Per-band masking: walk ``vrt.bands`` once and stream each + # band's mask. The first band with a sentinel hit promotes + # ``arr`` to float64 in place; ``int_arr`` keeps the original + # integer view alive so subsequent bands still compare against + # the exact sentinel dtype (the post-promotion float64 view + # works too, but staying on the integer dtype avoids any + # rounding edge case on extreme sentinels). Peak boolean-mask + # memory is O(H * W), not O(bands * H * W) like the earlier + # collect-then-apply implementation. + int_arr = arr + int_dtype = int_arr.dtype + for i, vrt_band in enumerate(vrt.bands): + if i >= int_arr.shape[-1]: + break + sentinel = _sentinel_for_dtype(vrt_band.nodata, int_dtype) + if sentinel is None: + continue + mask = int_arr[..., i] == sentinel + if not mask.any(): + continue + if arr.dtype != np.float64: + arr = arr.astype(np.float64) + arr[..., i][mask] = np.nan + elif nodata is not None: + sentinel = _sentinel_for_dtype(nodata, arr.dtype) + if sentinel is not None: + mask = arr == sentinel + if mask.any(): + arr = arr.astype(np.float64) + arr[mask] = np.nan # Surface the source GeoTransform in the same rasterio ordering used # by open_geotiff: (pixel_width, 0, origin_x, 0, pixel_height, origin_y). diff --git a/xrspatial/geotiff/tests/test_vrt_multiband_int_nodata_1611.py b/xrspatial/geotiff/tests/test_vrt_multiband_int_nodata_1611.py new file mode 100644 index 00000000..b392c630 --- /dev/null +++ b/xrspatial/geotiff/tests/test_vrt_multiband_int_nodata_1611.py @@ -0,0 +1,229 @@ +"""Regression tests for issue #1611. + +``read_vrt(path)`` on a multi-band integer VRT with per-band +```` tags used to only mask band 0's sentinel. Bands 1 +and up retained their integer sentinels as literal finite values in +the returned float64 array, breaking the convention that +``attrs['nodata']`` is present iff sentinel pixels are already NaN. + +The float-VRT path masks per-band correctly in ``_vrt._read_data`` +(lines 296-297, 347-351). For integer rasters the masking moved to +``__init__.py:read_vrt`` and used ``vrt.bands[0].nodata`` for every +band. The fix walks ``vrt.bands`` when ``band is None`` and masks +each ``arr[..., i]`` slice against its own ````. + +This file mirrors test_vrt_band_nodata_1598.py for the multi-band +``band=None`` path. PR #1602 fixed ``band=N`` single-band selection. +""" +from __future__ import annotations + +import numpy as np + +from xrspatial.geotiff import read_vrt +from xrspatial.geotiff._writer import write + + +def _write_two_band_per_band_nodata_vrt(tmp_path, *, dtype_str="UInt16", + np_dtype=np.uint16, + band0_sentinel=65535, + band1_sentinel=65000, + band0_other=(1, 2, 3), + band1_other=(7, 8, 9)): + """Two single-band integer sources, each with a distinct nodata + sentinel, exposed as bands 1 and 2 of a hand-rolled VRT. + + Used to be band 0's sentinel was the only one masked. Now every + band gets its own sentinel. + """ + b0_arr = np.array([[band0_other[0], band0_other[1]], + [band0_other[2], band0_sentinel]], dtype=np_dtype) + b1_arr = np.array([[band1_other[0], band1_other[1]], + [band1_other[2], band1_sentinel]], dtype=np_dtype) + p0 = str(tmp_path / 'vrt_b0_1611.tif') + p1 = str(tmp_path / 'vrt_b1_1611.tif') + write(b0_arr, p0, nodata=band0_sentinel, compression='none', + tiled=False) + write(b1_arr, p1, nodata=band1_sentinel, compression='none', + tiled=False) + + vrt_path = str(tmp_path / 'two_band_per_band_nodata_1611.vrt') + vrt_xml = f""" + 0.0, 1.0, 0.0, 0.0, 0.0, -1.0 + + {band0_sentinel} + + {p0} + 1 + + + + + + {band1_sentinel} + + {p1} + 1 + + + + +""" + with open(vrt_path, 'w') as f: + f.write(vrt_xml) + return vrt_path + + +def test_multiband_uint16_per_band_sentinel_each_masked(tmp_path): + """The previously-broken case: every band's sentinel must be NaN. + + Before the fix this returned dtype=float64 with band 0's (1,1) cell + as NaN but band 1's (1,1) cell as the literal 65000.0. + """ + vrt_path = _write_two_band_per_band_nodata_vrt(tmp_path) + r = read_vrt(vrt_path) + assert r.shape == (2, 2, 2) + assert r.dtype == np.float64, ( + f"expected float64 promotion, got {r.dtype}" + ) + assert np.isnan(r.values[1, 1, 0]), ( + "band 0's sentinel pixel was not NaN-masked." + ) + assert np.isnan(r.values[1, 1, 1]), ( + "band 1's sentinel pixel was not NaN-masked; " + "the regression from issue #1611 has returned." + ) + # Non-sentinel pixels survive unchanged + assert r.values[0, 0, 0] == 1 + assert r.values[0, 0, 1] == 7 + assert r.values[1, 0, 0] == 3 + assert r.values[1, 0, 1] == 9 + + +def test_multiband_int32_negative_per_band_sentinel(tmp_path): + """Negative sentinels in a signed integer VRT also mask per-band. + + The original bug was dtype-independent: any integer dtype with + per-band would have hit it. Cover int32 + negative + sentinels to make sure the helper handles signed types and the + range guard accepts negatives. + """ + vrt_path = _write_two_band_per_band_nodata_vrt( + tmp_path, dtype_str="Int32", np_dtype=np.int32, + band0_sentinel=-9999, band1_sentinel=-7777, + band0_other=(10, 20, 30), band1_other=(40, 50, 60)) + r = read_vrt(vrt_path) + assert r.dtype == np.float64 + assert np.isnan(r.values[1, 1, 0]) + assert np.isnan(r.values[1, 1, 1]) + assert r.values[0, 0, 0] == 10 + assert r.values[0, 0, 1] == 40 + + +def test_multiband_only_one_band_has_sentinel_present(tmp_path): + """If only one band's sentinel actually appears in the data, only + that band should change. The non-hitting band stays the same float64 + value (no spurious NaN introduced). + + Force band 1's sentinel never to appear by writing 99 instead. + """ + vrt_path = _write_two_band_per_band_nodata_vrt( + tmp_path, + band0_sentinel=65535, band1_sentinel=65000, + band1_other=(7, 8, 9)) + # Overwrite the band 1 source so the sentinel value 65000 is NOT + # present in band 1's actual data. + b1_no_sentinel = np.array([[7, 8], [9, 99]], dtype=np.uint16) + import os + p1 = os.path.join(os.path.dirname(vrt_path), 'vrt_b1_1611.tif') + write(b1_no_sentinel, p1, nodata=65000, compression='none', + tiled=False) + + r = read_vrt(vrt_path) + assert r.dtype == np.float64, ( + "Even when only band 0 has a present sentinel, the array still " + "needs promotion so band 0's NaN can be expressed." + ) + assert np.isnan(r.values[1, 1, 0]) + assert r.values[1, 1, 1] == 99.0 # band 1 sentinel absent, no NaN + + +def test_multiband_no_sentinel_present_anywhere_keeps_int_dtype(tmp_path): + """When no band actually contains its declared sentinel, skip + promotion entirely. Avoids a needless float64 cast on integer data. + """ + vrt_path = _write_two_band_per_band_nodata_vrt( + tmp_path, + band0_sentinel=65535, band1_sentinel=65000, + band0_other=(1, 2, 3), band1_other=(7, 8, 9)) + # Replace BOTH source files so neither contains its sentinel + import os + b0 = np.array([[1, 2], [3, 4]], dtype=np.uint16) + b1 = np.array([[7, 8], [9, 10]], dtype=np.uint16) + p0 = os.path.join(os.path.dirname(vrt_path), 'vrt_b0_1611.tif') + p1 = os.path.join(os.path.dirname(vrt_path), 'vrt_b1_1611.tif') + write(b0, p0, nodata=65535, compression='none', tiled=False) + write(b1, p1, nodata=65000, compression='none', tiled=False) + + r = read_vrt(vrt_path) + # Sentinels not present -> integer dtype preserved (matches the + # eager open_geotiff fast-path which also skips promotion when the + # mask is empty). + assert r.dtype == np.uint16 + assert r.values[1, 1, 0] == 4 + assert r.values[1, 1, 1] == 10 + + +def test_multiband_per_band_out_of_range_sentinel_is_no_op(tmp_path): + """A sentinel out of the integer dtype's range should be a no-op + for that band rather than raising. Mirrors PR #1583's behaviour + (#1581): the helper ``_int_nodata_in_range`` gates the cast. + """ + # uint16 cannot represent -9999; the helper should skip band 1. + vrt_path = _write_two_band_per_band_nodata_vrt( + tmp_path, + dtype_str="UInt16", np_dtype=np.uint16, + band0_sentinel=65535, band1_sentinel=10, # placeholder; rewrite below + band0_other=(1, 2, 3), band1_other=(7, 8, 9)) + + # Hand-rewrite the VRT so band 1's NoDataValue is the out-of-range -9999. + with open(vrt_path) as f: + xml = f.read() + xml = xml.replace("10", + "-9999") + with open(vrt_path, 'w') as f: + f.write(xml) + + # Should not raise and should still mask band 0. + r = read_vrt(vrt_path) + assert np.isnan(r.values[1, 1, 0]) # band 0 sentinel still masked + # Band 1's sentinel (-9999) is out of uint16 range; the value 10 + # in band1[1,1] survives unchanged. + assert r.values[1, 1, 1] == 10.0 or r.values[1, 1, 1] == 10 + + +def test_multiband_band_kwarg_still_per_band_post_pr1602(tmp_path): + """Non-regression check that PR #1602's band=N path still works. + + The fix here only changes the ``band is None`` branch; ``band=N`` + must still route through the single-band masking with its own + sentinel. + """ + vrt_path = _write_two_band_per_band_nodata_vrt(tmp_path) + r0 = read_vrt(vrt_path, band=0) + r1 = read_vrt(vrt_path, band=1) + assert r0.dtype == np.float64 + assert r1.dtype == np.float64 + assert r0.attrs.get('nodata') == 65535.0 + assert r1.attrs.get('nodata') == 65000.0 + assert np.isnan(r0.values[1, 1]) + assert np.isnan(r1.values[1, 1]) + + +def test_multiband_attrs_nodata_still_band0(tmp_path): + """``attrs['nodata']`` for band=None reads is documented as band + 0's sentinel (the canonical attr cannot encode per-band values). + The pixel-level fix must not change that contract. + """ + vrt_path = _write_two_band_per_band_nodata_vrt(tmp_path) + r = read_vrt(vrt_path) + assert r.attrs.get('nodata') == 65535.0