From 008fcb123e9aaccf2fe6c5c3ae567b38d3da3243 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Mon, 11 May 2026 10:28:48 -0700 Subject: [PATCH 1/3] Fix read_vrt(band=N) using band 0's nodata sentinel (#1598) read_vrt sourced attrs['nodata'] and the integer-with-nodata float64 promotion mask from vrt.bands[0].nodata regardless of which band was selected. A multi-band VRT with per-band tags would silently return band N's pixels with band 0's sentinel reported in attrs and band N's actual sentinel left as literal integers in the array (no NaN, no float64 promotion). Use vrt.bands[band].nodata when a band is selected, fall back to band 0 when no band kwarg was given (preserves the prior multi-band "first band wins" contract). Tests cover band=0 sanity, the broken band=1 case, and the multi-band fallback. --- .claude/sweep-metadata-state.csv | 2 +- xrspatial/geotiff/__init__.py | 17 ++- .../tests/test_vrt_band_nodata_1598.py | 107 ++++++++++++++++++ 3 files changed, 122 insertions(+), 4 deletions(-) create mode 100644 xrspatial/geotiff/tests/test_vrt_band_nodata_1598.py diff --git a/.claude/sweep-metadata-state.csv b/.claude/sweep-metadata-state.csv index a454a620..b5c3b65b 100644 --- a/.claude/sweep-metadata-state.csv +++ b/.claude/sweep-metadata-state.csv @@ -1,3 +1,3 @@ module,last_inspected,issue,severity_max,categories_found,notes -geotiff,2026-05-11,1580;1582,HIGH,1;3;4;5,"write_geotiff_gpu silently mis-ordered (band, y, x) DataArrays (#1580). to_geotiff ignored rioxarray nodatavals and CF _FillValue attrs (#1582). Both fixed in same branch." +geotiff,2026-05-11,1598,MEDIUM,4,"read_vrt(band=N) used vrt.bands[0].nodata for attrs and integer-promotion mask, mis-masking band N reads with per-band sentinels (#1598)." reproject,2026-05-10,1572;1573,HIGH,1;3;4,geoid_height_raster dropped input attrs and used dims[-2:] for 3D inputs (#1572). reproject/merge ignored nodatavals (rasterio convention) when rioxarray absent (#1573). Fixed in same branch. diff --git a/xrspatial/geotiff/__init__.py b/xrspatial/geotiff/__init__.py index 3e2a65b2..581ebc7b 100644 --- a/xrspatial/geotiff/__init__.py +++ b/xrspatial/geotiff/__init__.py @@ -2732,11 +2732,22 @@ def read_vrt(source: str, *, dtype=None, window=None, attrs['crs_wkt'] = vrt.crs_wkt if vrt.raster_type == 'point': attrs['raster_type'] = 'point' + # When a specific band is selected, source its nodata from that + # band's instead of band 0's. Otherwise multi-band + # VRTs with per-band sentinels would mis-mask the read: attrs would + # advertise band 0's sentinel, the integer-promotion block below + # would mask against band 0's sentinel, and band N's actual nodata + # pixels would survive as literal integers. See issue #1598. nodata = None if vrt.bands: - nodata = vrt.bands[0].nodata - if nodata is not None: - attrs['nodata'] = nodata + band_idx_for_nodata = band if band is not None else 0 + # ``_vrt.read_vrt`` already validates ``band`` is in range; the + # extra guard keeps a clearer message if a future refactor + # widens the public range without updating the internal reader. + if 0 <= band_idx_for_nodata < len(vrt.bands): + nodata = vrt.bands[band_idx_for_nodata].nodata + if nodata is not None: + attrs['nodata'] = nodata # Mirror the integer-with-nodata promotion that open_geotiff / # read_geotiff_dask / read_geotiff_gpu apply post-decode. The VRT diff --git a/xrspatial/geotiff/tests/test_vrt_band_nodata_1598.py b/xrspatial/geotiff/tests/test_vrt_band_nodata_1598.py new file mode 100644 index 00000000..a228d598 --- /dev/null +++ b/xrspatial/geotiff/tests/test_vrt_band_nodata_1598.py @@ -0,0 +1,107 @@ +"""Regression tests for issue #1598. + +``read_vrt(path, band=N)`` used to always source the nodata sentinel +from ``vrt.bands[0]`` rather than the requested band, so a multi-band +VRT with per-band ```` would mis-mask reads of any band +other than band 0: + +* ``attrs['nodata']`` advertised band 0's sentinel (wrong). +* The integer-to-float64 promotion ran against band 0's sentinel, so + band N's actual nodata pixels stayed as literal integers. +* The returned dtype was integer when it should have been float64. + +The fix uses ``vrt.bands[band].nodata`` when a band is selected. +""" +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): + """Two single-band uint16 sources, each with a distinct nodata + sentinel, exposed as bands 1 and 2 of a hand-rolled VRT. + """ + band0 = np.array([[1, 2], [3, 65535]], dtype=np.uint16) + band1 = np.array([[7, 8], [9, 65000]], dtype=np.uint16) + p0 = str(tmp_path / 'vrt_band0_1598.tif') + p1 = str(tmp_path / 'vrt_band1_1598.tif') + write(band0, p0, nodata=65535, compression='none', tiled=False) + write(band1, p1, nodata=65000, compression='none', tiled=False) + + vrt_path = str(tmp_path / 'two_band_per_band_nodata_1598.vrt') + vrt_xml = f""" + 0.0, 1.0, 0.0, 0.0, 0.0, -1.0 + + 65535 + + {p0} + 1 + + + + + + 65000 + + {p1} + 1 + + + + +""" + with open(vrt_path, 'w') as f: + f.write(vrt_xml) + return vrt_path + + +def test_read_vrt_band0_uses_band0_nodata(tmp_path): + """Sanity check the band-0 selection still works after the fix. + + Confirms the refactor did not flip the index. + """ + vrt_path = _write_two_band_per_band_nodata_vrt(tmp_path) + r = read_vrt(vrt_path, band=0) + assert r.dtype == np.float64 + assert r.attrs.get('nodata') == 65535.0 + assert np.isnan(r.values[1, 1]) + assert r.values[0, 0] == 1 + + +def test_read_vrt_band1_uses_band1_nodata(tmp_path): + """The previously-broken case: band=1 must use band 1's sentinel. + + Before the fix this returned dtype=uint16 with values=[[7,8], + [9,65000]] and attrs['nodata']=65535. + """ + vrt_path = _write_two_band_per_band_nodata_vrt(tmp_path) + r = read_vrt(vrt_path, band=1) + assert r.dtype == np.float64, ( + "band=1 read kept uint16 dtype; per-band nodata regression." + ) + assert r.attrs.get('nodata') == 65000.0, ( + f"attrs['nodata'] was {r.attrs.get('nodata')}, " + f"expected 65000 from band 1's ." + ) + assert np.isnan(r.values[1, 1]), ( + "band 1's sentinel pixel was not NaN-masked; " + "promotion ran against the wrong sentinel." + ) + assert r.values[0, 0] == 7 + assert r.values[1, 0] == 9 + + +def test_read_vrt_no_band_keeps_band0_nodata_attr(tmp_path): + """Unselected reads still surface band 0's sentinel. + + Multi-band VRTs with mixed sentinels return all bands stacked, and + the canonical attr cannot encode per-band values; advertising + band 0's sentinel matches the prior behavior and the documented + "first band wins" contract for multi-band reads. + """ + vrt_path = _write_two_band_per_band_nodata_vrt(tmp_path) + r = read_vrt(vrt_path) + assert r.attrs.get('nodata') == 65535.0 From 6a999af94f3685e6d0cf13a99f47e5da4ac97175 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Mon, 11 May 2026 11:24:33 -0700 Subject: [PATCH 2/3] Address PR #1602 review: validate band upstream in _vrt.read_vrt - Move the ``band`` range/type check into the internal ``_vrt.read_vrt`` reader where ``len(vrt.bands)`` is already in scope. Reject negative indices, non-int values, and out-of-range positive values with a clear ``ValueError`` that names the available band count. Previously these would either silently return a different band (Python negative indexing), raise an opaque ``IndexError`` deep in the read path, or raise ``TypeError`` on the list index. - Drop the now-redundant post-parse guard in the public ``read_vrt`` and update the comment to reflect that ``band`` is validated upstream. - Add three regression tests covering the negative, out-of-range, and non-int rejection paths. The original bug from #1598 (band=N reads pairing with band 0's nodata sentinel) was masked by these silent acceptances on the negative-band edge: ``read_vrt(path, band=-1)`` previously returned the last band's data with ``attrs['nodata']`` unset entirely. --- xrspatial/geotiff/__init__.py | 12 +++---- xrspatial/geotiff/_vrt.py | 17 ++++++++++ .../tests/test_vrt_band_nodata_1598.py | 34 +++++++++++++++++++ 3 files changed, 56 insertions(+), 7 deletions(-) diff --git a/xrspatial/geotiff/__init__.py b/xrspatial/geotiff/__init__.py index 581ebc7b..89245d9e 100644 --- a/xrspatial/geotiff/__init__.py +++ b/xrspatial/geotiff/__init__.py @@ -2738,16 +2738,14 @@ def read_vrt(source: str, *, dtype=None, window=None, # advertise band 0's sentinel, the integer-promotion block below # would mask against band 0's sentinel, and band N's actual nodata # pixels would survive as literal integers. See issue #1598. + # ``band`` has already been validated by ``_vrt.read_vrt`` as + # 0 <= band < len(vrt.bands), so a simple lookup is safe here. nodata = None if vrt.bands: band_idx_for_nodata = band if band is not None else 0 - # ``_vrt.read_vrt`` already validates ``band`` is in range; the - # extra guard keeps a clearer message if a future refactor - # widens the public range without updating the internal reader. - if 0 <= band_idx_for_nodata < len(vrt.bands): - nodata = vrt.bands[band_idx_for_nodata].nodata - if nodata is not None: - attrs['nodata'] = nodata + nodata = vrt.bands[band_idx_for_nodata].nodata + if nodata is not None: + attrs['nodata'] = nodata # Mirror the integer-with-nodata promotion that open_geotiff / # read_geotiff_dask / read_geotiff_gpu apply post-decode. The VRT diff --git a/xrspatial/geotiff/_vrt.py b/xrspatial/geotiff/_vrt.py index 7afba86e..a29ae767 100644 --- a/xrspatial/geotiff/_vrt.py +++ b/xrspatial/geotiff/_vrt.py @@ -242,6 +242,23 @@ def read_vrt(vrt_path: str, *, window=None, vrt_dir = os.path.dirname(os.path.abspath(vrt_path)) vrt = parse_vrt(xml_str, vrt_dir) + # Validate ``band`` against the parsed band count. Python list + # indexing would silently accept negative values (``vrt.bands[-1]`` + # returns the last band) and raise an opaque ``IndexError`` for + # out-of-range positive values, neither of which is the contract the + # public API documents (``band`` is a 0-based positive index). Reject + # both up front with a clear ``ValueError`` so callers do not get + # band-N data paired with band-0's nodata sentinel or a downstream + # IndexError from deep in the read path. + if band is not None: + if not isinstance(band, (int, np.integer)) or isinstance(band, bool): + raise ValueError( + f"band must be a non-negative int, got {band!r}") + if band < 0 or band >= len(vrt.bands): + raise ValueError( + f"band index {band} out of range for VRT with " + f"{len(vrt.bands)} band(s)") + if window is not None: r0, c0, r1, c1 = window r0 = max(0, r0) diff --git a/xrspatial/geotiff/tests/test_vrt_band_nodata_1598.py b/xrspatial/geotiff/tests/test_vrt_band_nodata_1598.py index a228d598..44fdea84 100644 --- a/xrspatial/geotiff/tests/test_vrt_band_nodata_1598.py +++ b/xrspatial/geotiff/tests/test_vrt_band_nodata_1598.py @@ -15,6 +15,7 @@ from __future__ import annotations import numpy as np +import pytest from xrspatial.geotiff import read_vrt from xrspatial.geotiff._writer import write @@ -105,3 +106,36 @@ def test_read_vrt_no_band_keeps_band0_nodata_attr(tmp_path): vrt_path = _write_two_band_per_band_nodata_vrt(tmp_path) r = read_vrt(vrt_path) assert r.attrs.get('nodata') == 65535.0 + + +def test_read_vrt_negative_band_raises(tmp_path): + """Negative band indices used to be silently accepted via Python + list indexing (``vrt.bands[-1]`` returned the last band) while the + public reader's nodata lookup rejected them, producing band-N data + with no nodata sentinel. They are now a clear ValueError up front. + """ + vrt_path = _write_two_band_per_band_nodata_vrt(tmp_path) + with pytest.raises(ValueError, match="band"): + read_vrt(vrt_path, band=-1) + + +def test_read_vrt_out_of_range_band_raises(tmp_path): + """Out-of-range band indices used to raise IndexError from deep in + the read path. They are now a ValueError that names the available + band count. + """ + vrt_path = _write_two_band_per_band_nodata_vrt(tmp_path) + with pytest.raises(ValueError, match="out of range"): + read_vrt(vrt_path, band=5) + + +def test_read_vrt_non_integer_band_raises(tmp_path): + """A non-int ``band`` would previously have raised TypeError on the + list index. ValueError here matches the rest of the input + validation surface. + """ + vrt_path = _write_two_band_per_band_nodata_vrt(tmp_path) + with pytest.raises(ValueError, match="band"): + read_vrt(vrt_path, band="1") + with pytest.raises(ValueError, match="band"): + read_vrt(vrt_path, band=True) From 865de2051bc0a65800af94ab124efdfc4b03b54c Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Mon, 11 May 2026 18:25:27 +0000 Subject: [PATCH 3/3] Merge origin/main into branch (resolve merge) Agent-Logs-Url: https://github.com/xarray-contrib/xarray-spatial/sessions/5e4ca757-8a1b-4021-aa2e-a9e7eab80367 Co-authored-by: brendancol <433221+brendancol@users.noreply.github.com> --- .claude/sweep-accuracy-state.csv | 2 +- .claude/sweep-metadata-state.csv | 1 + CONTRIBUTING.md | 43 +++- xrspatial/geotiff/__init__.py | 39 +++- .../tests/test_dask_int_nodata_chunks_1597.py | 116 ++++++++++ .../test_gpu_writer_nan_sentinel_1599.py | 213 ++++++++++++++++++ .../tests/test_streaming_codecs_2026_05_11.py | 8 +- 7 files changed, 414 insertions(+), 8 deletions(-) create mode 100644 xrspatial/geotiff/tests/test_dask_int_nodata_chunks_1597.py create mode 100644 xrspatial/geotiff/tests/test_gpu_writer_nan_sentinel_1599.py diff --git a/.claude/sweep-accuracy-state.csv b/.claude/sweep-accuracy-state.csv index 0ff28b72..9aeed5f8 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,1581,HIGH,2,"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,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." 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/.claude/sweep-metadata-state.csv b/.claude/sweep-metadata-state.csv index b5c3b65b..6f8ebcf6 100644 --- a/.claude/sweep-metadata-state.csv +++ b/.claude/sweep-metadata-state.csv @@ -1,3 +1,4 @@ module,last_inspected,issue,severity_max,categories_found,notes geotiff,2026-05-11,1598,MEDIUM,4,"read_vrt(band=N) used vrt.bands[0].nodata for attrs and integer-promotion mask, mis-masking band N reads with per-band sentinels (#1598)." +geotiff,2026-05-11,1597,HIGH,4;5,read_geotiff_dask dropped nodata mask when sentinel only hit non-first chunks (per-chunk dtype divergence; #1597). reproject,2026-05-10,1572;1573,HIGH,1;3;4,geoid_height_raster dropped input attrs and used dims[-2:] for 3D inputs (#1572). reproject/merge ignored nodatavals (rasterio convention) when rioxarray absent (#1573). Fixed in same branch. diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 14331987..c593a670 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -1,6 +1,6 @@ -# Contributing to Xarray-Spatial +# Contributing to xarray-spatial -As stated in [Xarray Spatial code of conduct](https://github.com/xarray-contrib/xarray-spatial/blob/main/CODE_OF_CONDUCT.md), a primary goal of Xarray Spatial is to be inclusive to the largest number of contributors. However, we do have some requests for how contributions should be made. Please read these guidelines before contributing to have a most positive experience with Xarray Spatial. +As stated in the [xarray-spatial code of conduct](https://github.com/xarray-contrib/xarray-spatial/blob/main/CODE_OF_CONDUCT.md), a primary goal of xarray-spatial is to be inclusive to the largest number of contributors. However, we do have some requests for how contributions should be made. Please read these guidelines before contributing to have a most positive experience with xarray-spatial. ### Getting Started @@ -28,7 +28,7 @@ In order to avoid duplication of effort, it's always a good idea to comment on a 1. Make sure that there is a corresponding issue for your change first. If there isn't yet, create one. -2. Create a fork of the Xarray Spatial repository on GitHub (this is only done before *first*) contribution. +2. Create a fork of the xarray-spatial repository on GitHub (this is only done before *first*) contribution. 3. Create a branch off the `main` branch with a meaningful name. Preferably include issue number and a few keywords, so that we will have a rough idea what the branch refers to, without looking up the issue. @@ -43,6 +43,43 @@ In order to avoid duplication of effort, it's always a good idea to comment on a 8. We will review your PR as time permits. Reviewers may comment on your contributions, ask you questions regarding the implementation or request changes. If changes are requested, push new commits to the existing branch. Do *NOT* rebase, amend, or cherry-pick published commits. Any of those actions will make us start the review from scratch. If you need updates from `main`, just merge it into your branch. +### AI-Assisted Contribution and Review + +xarray-spatial welcomes responsible AI-assisted development and review when it helps contributors move faster, improve quality, and reduce maintainer burden. The project has limited dedicated development capacity, so contributors are encouraged to use AI tools for testing, debugging, feature work, documentation, deployment, and release preparation. + +Attribution and responsibility belong to human contributors. Do not add AI tools, coding agents, or AI companies as commit authors, co-authors, reviewers, or attribution lines in commits, pull requests, changelogs, or release notes. + +Please avoid attribution such as: + + Co-authored-by: Claude <...> + Co-authored-by: ChatGPT <...> + Generated with Claude Code + Generated by Copilot + Reviewed-by: + +If AI-assisted code introduces a bug, security issue, regression, or incorrect behavior, the responsible party is the human contributor who submitted and reviewed the change. AI tools may assist the work, but they do not own the work. + +#### Review Expectations + +Before submitting a pull request, contributors should review and use the relevant commands in: + + .claude/commands/ + +These command markdown files reflect current project expectations for performance, accuracy, security, testing, maintainability, deployment, and release readiness. New contributors should read through them to understand the project's quality standards. + +Contributors who do not use Claude Code can adapt the command markdown files into prompts or checklists for their preferred AI tools and review workflows. The important part is the quality review they represent, not the specific tool used to run them. + +Where applicable, run the relevant sweep or review commands before requesting maintainer review. If a command flags an issue, either address it or clearly explain why it is acceptable for the current change. + +Copilot code review may run on pull requests as automated review assistance. Its feedback is useful, but it does not replace human responsibility or maintainer judgment. + +#### Current Priority + +The next sprint is focused on stabilizing xarray-spatial for its first major release. Contributors should prioritize bug fixes, reliability, tests, release blockers, and performance, accuracy, or security issues. New features should avoid adding significant release risk unless they are essential. + +The preferred contribution style is fast, careful, and human-accountable: use AI aggressively to improve velocity and quality, but keep ownership and attribution with the people contributing to the project. + + ### Attribution Portions of text derived from [Bokeh CONTRIBUTING file]: (https://github.com/bokeh/bokeh/blob/branch-3.4/.github/CONTRIBUTING.md) diff --git a/xrspatial/geotiff/__init__.py b/xrspatial/geotiff/__init__.py index 89245d9e..e66e6dc9 100644 --- a/xrspatial/geotiff/__init__.py +++ b/xrspatial/geotiff/__init__.py @@ -1641,13 +1641,23 @@ def read_geotiff_dask(source: str, *, dtype=None, chunks: int | tuple = 512, # Translate window-relative chunk coords back to file-relative # coords for ``read_to_array``. ``win_r0`` / ``win_c0`` are 0 # when no window was requested. + # Always thread ``target_dtype`` so each delayed chunk lands + # in the same dtype that the dask array declared. Without + # this, an integer raster with an in-range nodata sentinel + # would have ``effective_dtype=float64`` declared on the + # dask graph but per-chunk arrays only promoted when a + # chunk happened to contain a sentinel pixel. Dask + # concatenation then preallocates from the first chunk's + # actual dtype (uint16), silently casting later float64 + # chunks back to int and converting their NaNs to 0. See + # issue #1597. block = da.from_delayed( _delayed_read_window(source, r0 + win_r0, c0 + win_c0, r1 + win_r0, c1 + win_c0, overview_level, nodata, band_arg, - target_dtype=target_dtype if dtype is not None else None, + target_dtype=target_dtype, http_meta_key=http_meta_key, max_pixels=max_pixels), shape=block_shape, @@ -2583,6 +2593,33 @@ def write_geotiff_gpu(data, path: str, *, samples = arr.shape[2] if arr.ndim == 3 else 1 np_dtype = np.dtype(str(arr.dtype)) # cupy dtype -> numpy dtype + # Mirror the CPU writer's NaN-to-sentinel substitution (issue #1599). + # Without this step the GPU writer emits raw NaN bytes interleaved + # with valid data even when ``nodata=`` is supplied; the + # GDAL_NODATA tag still advertises the sentinel but external readers + # (rasterio / GDAL / QGIS) mask only on the sentinel value and + # therefore see the NaN pixels as valid data. The CPU writer does + # the equivalent rewrite at ``to_geotiff`` (lines around + # ``arr.copy(); arr[nan_mask] = arr.dtype.type(nodata)``); both + # paths must produce byte-equivalent files for the same input. + # We always copy before the in-place sentinel write. Some upstream + # branches above already produce a fresh buffer (``cupy.asarray`` + # from numpy/dask, ``ascontiguousarray`` from the band-first + # moveaxis); others (a CuPy-backed DataArray taking the no-moveaxis + # path, or a plain CuPy positional ``data``) hand ``arr`` back as + # the caller's buffer. Rather than tracking provenance across that + # branch tree, copy unconditionally when we are about to mutate -- + # the cost is one GPU array allocation, only on the NaN-present + # path, and it guarantees the CPU writer's defensive-copy semantics + # in every case. + if (nodata is not None + and np_dtype.kind == 'f' + and not np.isnan(float(nodata))): + nan_mask = cupy.isnan(arr) + if bool(nan_mask.any()): + arr = arr.copy() + arr[nan_mask] = np_dtype.type(nodata) + comp_tag = _compression_tag(compression) pred_val = normalize_predictor(predictor, np_dtype, comp_tag) diff --git a/xrspatial/geotiff/tests/test_dask_int_nodata_chunks_1597.py b/xrspatial/geotiff/tests/test_dask_int_nodata_chunks_1597.py new file mode 100644 index 00000000..d86e1789 --- /dev/null +++ b/xrspatial/geotiff/tests/test_dask_int_nodata_chunks_1597.py @@ -0,0 +1,116 @@ +"""Regression tests for issue #1597. + +``read_geotiff_dask`` on an integer raster with an in-range nodata +sentinel used to silently lose the mask when the sentinel only appeared +in non-first chunks. Per-chunk dtype divergence (uint16 vs float64) +caused dask concatenation to preallocate from the first chunk's actual +dtype, casting float64 chunks back to int and converting NaN to 0. + +The fix threads the resolved ``target_dtype`` (the dask graph's +declared dtype) unconditionally through ``_delayed_read_window`` so +every chunk lands as float64 regardless of whether its mask hit. +""" +from __future__ import annotations + +import numpy as np +import pytest + +from xrspatial.geotiff import open_geotiff +from xrspatial.geotiff._writer import write + + +@pytest.fixture +def uint16_with_sentinel_only_in_corner(tmp_path): + """Write a uint16 8x8 TIFF whose nodata sentinel is in the + bottom-right 2x2 quadrant. With ``chunks=4`` the top-left chunk + never sees a sentinel and used to keep its uint16 dtype. + """ + arr = np.arange(64, dtype=np.uint16).reshape(8, 8) + 1 + arr[6:8, 6:8] = 65535 + path = str(tmp_path / 'uint16_corner_sentinel_1597.tif') + write(arr, path, nodata=65535, compression='none', tiled=False) + return path, arr + + +def test_eager_promotes_to_float64_and_masks(uint16_with_sentinel_only_in_corner): + """Baseline: the eager path produces float64 with 4 NaNs.""" + path, _ = uint16_with_sentinel_only_in_corner + eager = open_geotiff(path) + assert eager.dtype == np.float64 + assert np.isnan(eager.values).sum() == 4 + assert np.isnan(eager.values[6:8, 6:8]).all() + + +def test_dask_chunks_4_matches_eager(uint16_with_sentinel_only_in_corner): + """The dask compute result matches the eager path bit-for-bit. + + Before the fix this returned a uint16 array with 0s where the + sentinel had been, because dask coerced the late-arriving float64 + chunk back to uint16 at concat time. + """ + path, _ = uint16_with_sentinel_only_in_corner + eager = open_geotiff(path) + dk = open_geotiff(path, chunks=4) + assert dk.dtype == np.float64 + computed = dk.compute() + assert computed.dtype == np.float64 + np.testing.assert_array_equal(np.isnan(computed.values), + np.isnan(eager.values)) + finite = ~np.isnan(eager.values) + np.testing.assert_array_equal(computed.values[finite], + eager.values[finite]) + + +def test_dask_chunks_2_per_chunk_dtype_uniform( + uint16_with_sentinel_only_in_corner): + """Every dask chunk returns float64 regardless of mask hit. + + Iterates the delayed blocks and asserts each one computes to + float64; the regression had the first chunk's actual data come back + as uint16 because the mask never matched there. + """ + path, _ = uint16_with_sentinel_only_in_corner + dk = open_geotiff(path, chunks=2) + blocks = dk.data.to_delayed().flatten() + for i, block in enumerate(blocks): + chunk = block.compute() + assert chunk.dtype == np.float64, ( + f"chunk {i} computed as {chunk.dtype}, expected float64; " + f"per-chunk dtype divergence is the #1597 regression." + ) + + +def test_dask_keeps_dtype_for_out_of_range_sentinel(tmp_path): + """Out-of-range sentinels (uint16 + nodata=-9999) stay uint16. + + The fix should not regress #1581: when the sentinel cannot match + any pixel, no float64 promotion is needed and the dask path keeps + the file's native dtype. + """ + arr = np.array([[1, 2, 3, 4]] * 4, dtype=np.uint16) + path = str(tmp_path / 'uint16_out_of_range_1597.tif') + write(arr, path, nodata=-9999, compression='none', tiled=False) + + dk = open_geotiff(path, chunks=2) + assert dk.dtype == np.uint16 + result = dk.compute() + assert result.dtype == np.uint16 + np.testing.assert_array_equal(result.values, arr) + + +def test_dask_float_input_with_sentinel_in_one_chunk(tmp_path): + """Float rasters with sentinel in non-first chunk also stay float. + + The float path doesn't promote dtype, but it does in-place NaN + substitution. Verify the substitution holds for chunks with and + without the sentinel. + """ + arr = np.arange(64, dtype=np.float32).reshape(8, 8) + 1 + arr[6:8, 6:8] = -9999.0 + path = str(tmp_path / 'float_corner_sentinel_1597.tif') + write(arr, path, nodata=-9999, compression='none', tiled=False) + + eager = open_geotiff(path) + dk = open_geotiff(path, chunks=4).compute() + np.testing.assert_array_equal(np.isnan(dk.values), + np.isnan(eager.values)) diff --git a/xrspatial/geotiff/tests/test_gpu_writer_nan_sentinel_1599.py b/xrspatial/geotiff/tests/test_gpu_writer_nan_sentinel_1599.py new file mode 100644 index 00000000..8b420379 --- /dev/null +++ b/xrspatial/geotiff/tests/test_gpu_writer_nan_sentinel_1599.py @@ -0,0 +1,213 @@ +"""Regression tests for issue #1599. + +``write_geotiff_gpu`` (and ``to_geotiff(..., gpu=True)``) used to emit +raw NaN bytes for missing pixels even when ``nodata=`` was +supplied, while the CPU writer substituted NaN with the sentinel value +before encoding. The mismatch was invisible on xrspatial-only round +trips (the reader masks both NaN and the sentinel) but external +readers that mask only on the GDAL_NODATA value -- rasterio, GDAL, +QGIS -- treated the NaN pixels in GPU-written files as valid data. + +The GPU writer now mirrors the CPU writer's NaN-to-sentinel rewrite so +both backends produce byte-equivalent files for the same input. +""" +from __future__ import annotations + +import importlib.util + +import numpy as np +import pytest +import xarray as xr + +from xrspatial.geotiff import to_geotiff, write_geotiff_gpu +from xrspatial.geotiff._reader import read_to_array + + +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 _float_raster_with_nans(dtype=np.float32): + h, w = 50, 60 + arr = np.arange(h * w, dtype=dtype).reshape(h, w) + # Two disconnected NaN patches so single-pixel guards don't accidentally + # mask the bug. + arr[10:15, 20:25] = np.nan + arr[40, 50] = np.nan + y = np.arange(h, dtype=np.float64) * -0.1 + 50.0 + x = np.arange(w, dtype=np.float64) * 0.1 - 100.0 + return arr, y, x + + +@_gpu_only +def test_gpu_writer_substitutes_nan_with_sentinel(tmp_path): + """GPU-written file stores the sentinel where the input held NaN.""" + import cupy as cp + + arr, y, x = _float_raster_with_nans() + da = xr.DataArray(cp.asarray(arr.copy()), dims=['y', 'x'], + coords={'y': y, 'x': x}) + p = str(tmp_path / "gpu_nan_sentinel.tif") + write_geotiff_gpu(da, p, crs=4326, nodata=-9999) + + raw, _ = read_to_array(p) + # The raw decoded bytes should contain the sentinel, not NaN. + assert not np.isnan(raw).any(), \ + "GPU writer left NaN in file bytes despite nodata=-9999" + nan_locations = np.isnan(arr) + assert (raw[nan_locations] == -9999.0).all() + # And the valid pixels round-trip untouched. + assert np.array_equal(raw[~nan_locations], arr[~nan_locations]) + + +@_gpu_only +def test_gpu_and_cpu_writers_byte_equivalent_on_nan_input(tmp_path): + """For a float NaN input + finite sentinel, GPU and CPU writers + must produce identical pixel data.""" + import cupy as cp + + arr, y, x = _float_raster_with_nans() + da_cpu = xr.DataArray(arr.copy(), dims=['y', 'x'], + coords={'y': y, 'x': x}) + da_gpu = xr.DataArray(cp.asarray(arr.copy()), dims=['y', 'x'], + coords={'y': y, 'x': x}) + + p_cpu = str(tmp_path / "cpu.tif") + p_gpu = str(tmp_path / "gpu.tif") + to_geotiff(da_cpu, p_cpu, crs=4326, nodata=-9999) + write_geotiff_gpu(da_gpu, p_gpu, crs=4326, nodata=-9999) + + raw_cpu, _ = read_to_array(p_cpu) + raw_gpu, _ = read_to_array(p_gpu) + assert np.array_equal(raw_cpu, raw_gpu) + + +@_gpu_only +def test_gpu_writer_preserves_caller_cupy_buffer(tmp_path): + """The NaN-to-sentinel rewrite must NOT mutate the caller's CuPy + array. The CPU writer takes a defensive copy at the matching step; + the GPU writer mirrors that behaviour.""" + import cupy as cp + + arr, y, x = _float_raster_with_nans() + cp_arr = cp.asarray(arr.copy()) + da = xr.DataArray(cp_arr, dims=['y', 'x'], coords={'y': y, 'x': x}) + + p = str(tmp_path / "gpu_preserve.tif") + write_geotiff_gpu(da, p, crs=4326, nodata=-9999) + + after = cp.asnumpy(cp_arr) + # NaNs must still be present at the original locations. + assert np.isnan(after[10:15, 20:25]).all() + assert np.isnan(after[40, 50]) + # And non-NaN cells must be unchanged. + np.testing.assert_array_equal( + after[~np.isnan(arr)], arr[~np.isnan(arr)]) + + +@_gpu_only +def test_gpu_writer_no_rewrite_when_no_nans(tmp_path): + """A NaN-free input must round-trip bit-exact with the GPU writer + regardless of the nodata kwarg.""" + import cupy as cp + + arr, y, x = _float_raster_with_nans() + # Replace NaNs with finite values so the array has no NaN at all. + arr = np.where(np.isnan(arr), 1.5, arr).astype(np.float32) + assert not np.isnan(arr).any() + da = xr.DataArray(cp.asarray(arr.copy()), dims=['y', 'x'], + coords={'y': y, 'x': x}) + p = str(tmp_path / "gpu_no_nans.tif") + write_geotiff_gpu(da, p, crs=4326, nodata=-9999) + + raw, _ = read_to_array(p) + assert np.array_equal(raw, arr) + + +@_gpu_only +def test_gpu_writer_nan_nodata_skips_substitution(tmp_path): + """If the requested sentinel is NaN itself, no rewrite is needed. + The file keeps the NaN bytes verbatim (CPU and GPU agree).""" + import cupy as cp + + arr, y, x = _float_raster_with_nans() + da = xr.DataArray(cp.asarray(arr.copy()), dims=['y', 'x'], + coords={'y': y, 'x': x}) + p = str(tmp_path / "gpu_nan_sentinel_nan.tif") + write_geotiff_gpu(da, p, crs=4326, nodata=float('nan')) + + raw, _ = read_to_array(p) + # NaN pixels remain NaN; finite pixels remain finite. + nan_locations = np.isnan(arr) + assert np.isnan(raw[nan_locations]).all() + np.testing.assert_array_equal(raw[~nan_locations], arr[~nan_locations]) + + +@_gpu_only +def test_gpu_writer_external_reader_sees_correct_nodata_mask(tmp_path): + """rasterio (and any other GDAL_NODATA-strict reader) must see the + same valid-pixel set on CPU and GPU outputs. This is the bug from + #1599: the GPU file used to report 100% valid pixels because the + sentinel was never written into NaN positions.""" + rasterio = pytest.importorskip("rasterio") + import cupy as cp + + arr, y, x = _float_raster_with_nans() + n_nans = int(np.isnan(arr).sum()) + da_cpu = xr.DataArray(arr.copy(), dims=['y', 'x'], + coords={'y': y, 'x': x}) + da_gpu = xr.DataArray(cp.asarray(arr.copy()), dims=['y', 'x'], + coords={'y': y, 'x': x}) + + p_cpu = str(tmp_path / "cpu.tif") + p_gpu = str(tmp_path / "gpu.tif") + # Use deflate rather than the default zstd: some rasterio/GDAL builds + # ship without ZSTD codec support and would fail this round-trip + # for environment reasons unrelated to the nodata mask under test. + to_geotiff(da_cpu, p_cpu, crs=4326, nodata=-9999, compression='deflate') + write_geotiff_gpu( + da_gpu, p_gpu, crs=4326, nodata=-9999, compression='deflate') + + with rasterio.open(p_cpu) as src: + cpu_masked = src.read(1, masked=True) + with rasterio.open(p_gpu) as src: + gpu_masked = src.read(1, masked=True) + + expected_invalid = n_nans + assert cpu_masked.size - cpu_masked.count() == expected_invalid + assert gpu_masked.size - gpu_masked.count() == expected_invalid + + +@_gpu_only +def test_gpu_writer_multiband_nan_substitution(tmp_path): + """The substitution must work for 3D (y, x, band) inputs as well.""" + import cupy as cp + + h, w, b = 30, 40, 3 + arr = np.arange(h * w * b, dtype=np.float32).reshape(h, w, b) + arr[5:10, 8:12, 1] = np.nan + arr[20, 30, 0] = np.nan + y = np.arange(h, dtype=np.float64) * -0.1 + 50.0 + x = np.arange(w, dtype=np.float64) * 0.1 - 100.0 + da = xr.DataArray(cp.asarray(arr.copy()), dims=['y', 'x', 'band'], + coords={'y': y, 'x': x, 'band': np.arange(b)}) + + p = str(tmp_path / "gpu_mb.tif") + write_geotiff_gpu(da, p, crs=4326, nodata=-9999) + + raw, _ = read_to_array(p) + nan_locations = np.isnan(arr) + assert not np.isnan(raw).any() + assert (raw[nan_locations] == -9999.0).all() + np.testing.assert_array_equal( + raw[~nan_locations], arr[~nan_locations]) diff --git a/xrspatial/geotiff/tests/test_streaming_codecs_2026_05_11.py b/xrspatial/geotiff/tests/test_streaming_codecs_2026_05_11.py index 5339b432..25bda4ad 100644 --- a/xrspatial/geotiff/tests/test_streaming_codecs_2026_05_11.py +++ b/xrspatial/geotiff/tests/test_streaming_codecs_2026_05_11.py @@ -9,8 +9,10 @@ * Dask streaming write + LERC (lossless and lossy ``max_z_error``) produces identical output to the eager writer. * Dask streaming write + LZ4 round-trips a float32 raster. -* Dask streaming write + packbits round-trips a uint8 raster (the only - dtype packbits supports in this writer). +* Dask streaming write + packbits round-trips a uint8 raster. packbits + operates on the raw byte stream so any dtype is technically supported; + uint8 is the variant exercised here and the one packbits was designed + for (run-length encoding of byte runs). * COG output with ``overview_resampling='cubic'`` round-trips through scipy.ndimage.zoom (the only overview method that takes a separate code path in ``_block_reduce_2d``). @@ -139,7 +141,7 @@ def test_streaming_matches_eager(self, float_raster, dask_float_raster, # --------------------------------------------------------------------------- -# packbits: round-trip on uint8 (the only supported dtype) +# packbits: round-trip on uint8 (the dtype packbits is designed for) # --------------------------------------------------------------------------- class TestStreamingPackbits: