Apply nodata mask in read_geotiff_gpu (#1542)

xrspatial/geotiff/__init__.py

@@ -563,6 +563,41 @@ def _is_gpu_data(data) -> bool:
     return isinstance(data, _cupy_type)
 
 
+def _apply_nodata_mask_gpu(arr_gpu, nodata):
+    """Replace nodata sentinel values with NaN on a CuPy array.
+
+    Mirrors the CPU eager path in :func:`open_geotiff` (lines around the
+    ``# Apply nodata mask`` comment). Float arrays get NaN substituted in
+    place of the sentinel; integer arrays are promoted to float64 with NaN
+    where the sentinel was. NaN nodata on a float array is a no-op (the
+    sentinel already matches NaN-aware code).
+
+    Returns the (possibly promoted, possibly nodata-masked) CuPy array.
+    The caller is responsible for setting ``attrs['nodata']`` so the
+    sentinel is still discoverable downstream.
+    """
+    import cupy
+
+    if nodata is None:
+        return arr_gpu
+    arr_dtype = np.dtype(str(arr_gpu.dtype))
+    if arr_dtype.kind == 'f':
+        if not np.isnan(nodata):
+            sentinel = arr_dtype.type(nodata)
+            arr_gpu = cupy.where(arr_gpu == sentinel,
+                                 arr_dtype.type('nan'), arr_gpu)
+        return arr_gpu
+    if arr_dtype.kind in ('u', 'i'):
+        nodata_int = int(nodata)
+        sentinel = arr_dtype.type(nodata_int)
+        mask = arr_gpu == sentinel
+        if bool(mask.any().item()):
+            arr_gpu = arr_gpu.astype(cupy.float64)
+            arr_gpu = cupy.where(mask, cupy.float64('nan'), arr_gpu)
+        return arr_gpu
+    return arr_gpu
+
+
 _LEVEL_RANGES = {
     'deflate': (1, 9),
     'zstd': (1, 22),
@@ -1665,6 +1700,14 @@ def read_geotiff_gpu(source: str, *,
             t_tuple = _transform_tuple(geo_info)
             if t_tuple is not None:
                 attrs['transform'] = t_tuple
+            # Apply nodata mask + record sentinel so the GPU read agrees
+            # with the CPU eager path. Without this, integer rasters keep
+            # the literal sentinel value and float rasters keep the
+            # sentinel rather than NaN -- a silent backend divergence.
+            nodata = geo_info.nodata
+            if nodata is not None:
+                attrs['nodata'] = nodata
+                arr_gpu = _apply_nodata_mask_gpu(arr_gpu, nodata)
             if dtype is not None:
                 target = np.dtype(dtype)
                 _validate_dtype_cast(np.dtype(str(arr_gpu.dtype)), target)
@@ -1866,6 +1909,16 @@ def _read_once():
                 f"({height}, {width}, {samples})"
             )
 
+    # Apply nodata mask + record sentinel so the GPU read agrees with the
+    # CPU eager path (issue #1542). Without this, integer rasters keep the
+    # literal sentinel value and float rasters keep the sentinel rather
+    # than NaN -- a silent backend divergence. Apply before the optional
+    # dtype cast so the float promotion for masked integer rasters doesn't
+    # surprise a user-supplied dtype.
+    nodata = geo_info.nodata
+    if nodata is not None:
+        arr_gpu = _apply_nodata_mask_gpu(arr_gpu, nodata)
+
     if dtype is not None:
         target = np.dtype(dtype)
         _validate_dtype_cast(np.dtype(str(arr_gpu.dtype)), target)
@@ -1886,6 +1939,8 @@ def _read_once():
     t_tuple = _transform_tuple(geo_info)
     if t_tuple is not None:
         attrs['transform'] = t_tuple
+    if nodata is not None:
+        attrs['nodata'] = nodata
 
     if arr_gpu.ndim == 3:
         dims = ['y', 'x', 'band']

xrspatial/geotiff/tests/test_gpu_nodata_1542.py

@@ -0,0 +1,197 @@
+"""GPU read backend nodata propagation tests for issue #1542.
+
+Before the fix, ``read_geotiff_gpu`` silently differed from the CPU
+eager path on rasters with a declared nodata sentinel:
+
+* The returned DataArray's ``attrs`` had no ``nodata`` key, even when
+  the file declared one via the GDAL_NODATA tag.
+* The pixel data still carried the raw sentinel: integer rasters were
+  not promoted to float64 with NaN, and float rasters kept the
+  sentinel rather than NaN.
+
+These tests pin the contract that the GPU read agrees with the CPU
+read on dtype, NaN positions, and ``attrs['nodata']`` for the
+combinations the audit found broken.
+"""
+from __future__ import annotations
+
+import importlib.util
+import os
+import tempfile
+
+import numpy as np
+import pytest
+
+
+def _gpu_available() -> bool:
+    """True if cupy is importable and CUDA is initialised."""
+    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",
+)
+
+
+@_gpu_only
+def test_gpu_uint16_nodata_promoted_and_masked_tiled(tmp_path):
+    """uint16 + nodata sentinel -> GPU returns float64 with NaN, attrs[nodata] set."""
+    from xrspatial.geotiff import open_geotiff
+    from xrspatial.geotiff._writer import write
+
+    arr = np.array([[1, 2, 3], [65535, 5, 6]], dtype=np.uint16)
+    path = str(tmp_path / 'gpu_u16_nodata_1542_tiled.tif')
+    write(arr, path, nodata=65535, compression='deflate',
+          tiled=True, tile_size=16)
+
+    cpu = open_geotiff(path)
+    gpu = open_geotiff(path, gpu=True)
+
+    assert cpu.dtype == gpu.dtype == np.float64
+    assert cpu.attrs.get('nodata') == 65535.0
+    assert gpu.attrs.get('nodata') == 65535.0
+    cpu_arr = cpu.values
+    gpu_arr = gpu.data.get()
+    np.testing.assert_array_equal(np.isnan(cpu_arr), np.isnan(gpu_arr))
+    finite = ~np.isnan(cpu_arr)
+    np.testing.assert_array_equal(cpu_arr[finite], gpu_arr[finite])
+
+
+@_gpu_only
+def test_gpu_uint16_nodata_promoted_and_masked_stripped(tmp_path):
+    """Stripped fallback path also promotes + masks integer nodata."""
+    from xrspatial.geotiff import open_geotiff
+    from xrspatial.geotiff._writer import write
+
+    arr = np.array([[1, 2, 3], [65535, 5, 6]], dtype=np.uint16)
+    path = str(tmp_path / 'gpu_u16_nodata_1542_stripped.tif')
+    write(arr, path, nodata=65535, compression='deflate', tiled=False)
+
+    cpu = open_geotiff(path)
+    gpu = open_geotiff(path, gpu=True)
+
+    assert cpu.dtype == gpu.dtype == np.float64
+    assert gpu.attrs.get('nodata') == 65535.0
+    np.testing.assert_array_equal(np.isnan(cpu.values),
+                                  np.isnan(gpu.data.get()))
+
+
+@_gpu_only
+def test_gpu_float32_sentinel_replaced_with_nan(tmp_path):
+    """float32 + finite sentinel -> GPU returns float32 with NaN at sentinel positions."""
+    from xrspatial.geotiff import open_geotiff
+    from xrspatial.geotiff._writer import write
+
+    arr = np.arange(24, dtype=np.float32).reshape(4, 6)
+    arr[0, 0] = -9999.0
+    arr[2, 3] = -9999.0
+    path = str(tmp_path / 'gpu_f32_sentinel_1542.tif')
+    write(arr, path, nodata=-9999.0, compression='deflate',
+          tiled=True, tile_size=16)
+
+    cpu = open_geotiff(path)
+    gpu = open_geotiff(path, gpu=True)
+
+    assert cpu.dtype == gpu.dtype == np.float32
+    assert gpu.attrs.get('nodata') == -9999.0
+    cpu_arr = cpu.values
+    gpu_arr = gpu.data.get()
+    np.testing.assert_array_equal(np.isnan(cpu_arr), np.isnan(gpu_arr))
+    finite = ~np.isnan(cpu_arr)
+    np.testing.assert_array_equal(cpu_arr[finite], gpu_arr[finite])
+
+
+@_gpu_only
+def test_gpu_no_nodata_keeps_dtype(tmp_path):
+    """No nodata declared -> GPU keeps source dtype, no nodata attr added."""
+    from xrspatial.geotiff import open_geotiff
+    from xrspatial.geotiff._writer import write
+
+    arr = np.array([[1, 2, 3], [4, 5, 6]], dtype=np.uint16)
+    path = str(tmp_path / 'gpu_u16_no_nodata_1542.tif')
+    write(arr, path, compression='deflate', tiled=True, tile_size=16)
+
+    gpu = open_geotiff(path, gpu=True)
+    assert gpu.dtype == np.uint16
+    assert gpu.attrs.get('nodata') is None
+    np.testing.assert_array_equal(gpu.data.get(), arr)
+
+
+@_gpu_only
+def test_gpu_nan_nodata_passes_through(tmp_path):
+    """nodata=NaN on float data -> GPU returns NaN positions intact, attrs[nodata]=nan."""
+    from xrspatial.geotiff import open_geotiff
+    from xrspatial.geotiff._writer import write
+
+    arr = np.array([[1.0, 2.0, 3.0], [np.nan, 5.0, 6.0]], dtype=np.float32)
+    path = str(tmp_path / 'gpu_f32_nan_1542.tif')
+    write(arr, path, nodata=float('nan'), compression='deflate',
+          tiled=True, tile_size=16)
+
+    gpu = open_geotiff(path, gpu=True)
+    assert np.isnan(gpu.attrs.get('nodata'))
+    gpu_arr = gpu.data.get()
+    assert np.isnan(gpu_arr[1, 0])
+    finite = ~np.isnan(gpu_arr)
+    np.testing.assert_array_equal(gpu_arr[finite], arr[~np.isnan(arr)])
+
+
+@_gpu_only
+def test_gpu_all_four_backends_agree_on_nodata(tmp_path):
+    """numpy / dask+numpy / cupy / dask+cupy all agree on dtype + nodata + NaN."""
+    from xrspatial.geotiff import open_geotiff
+    from xrspatial.geotiff._writer import write
+
+    arr = np.array([[1, 2, 3], [65535, 5, 6]], dtype=np.uint16)
+    path = str(tmp_path / 'gpu_4backends_1542.tif')
+    write(arr, path, nodata=65535, compression='deflate',
+          tiled=True, tile_size=16)
+
+    da_np = open_geotiff(path)
+    da_dask = open_geotiff(path, chunks=512)
+    da_gpu = open_geotiff(path, gpu=True)
+    da_gpu_dask = open_geotiff(path, gpu=True, chunks=512)
+
+    # dtype
+    for label, da in [('np', da_np), ('dask+np', da_dask),
+                      ('gpu', da_gpu), ('gpu+dask', da_gpu_dask)]:
+        assert da.dtype == np.float64, f"{label}: dtype={da.dtype}"
+        assert da.attrs.get('nodata') == 65535.0, (
+            f"{label}: missing nodata attr (got {da.attrs.get('nodata')!r})"
+        )
+
+    # NaN positions
+    np_arr = da_np.values
+    dask_arr = da_dask.compute().values
+    gpu_arr = da_gpu.data.get()
+    gpu_dask_arr = da_gpu_dask.compute().data.get()
+    np.testing.assert_array_equal(np.isnan(np_arr), np.isnan(dask_arr))
+    np.testing.assert_array_equal(np.isnan(np_arr), np.isnan(gpu_arr))
+    np.testing.assert_array_equal(np.isnan(np_arr), np.isnan(gpu_dask_arr))
+
+
+@_gpu_only
+def test_gpu_int16_negative_nodata(tmp_path):
+    """Signed integer with negative nodata: also promoted to float64 + NaN."""
+    from xrspatial.geotiff import open_geotiff
+    from xrspatial.geotiff._writer import write
+
+    arr = np.array([[-9999, 10], [20, -9999]], dtype=np.int16)
+    path = str(tmp_path / 'gpu_i16_nodata_1542.tif')
+    write(arr, path, nodata=-9999, compression='deflate',
+          tiled=True, tile_size=16)
+
+    cpu = open_geotiff(path)
+    gpu = open_geotiff(path, gpu=True)
+    assert cpu.dtype == gpu.dtype == np.float64
+    assert gpu.attrs.get('nodata') == -9999.0
+    np.testing.assert_array_equal(np.isnan(cpu.values),
+                                  np.isnan(gpu.data.get()))

xrspatial/geotiff/tests/test_lerc_valid_mask_gpu.py

@@ -80,7 +80,17 @@ def _patched(data, width, height, samples=1,
 
 
 def _read_cpu_gpu(path):
-    """Read *path* with both readers and return ``(cpu_array, gpu_host_array)``."""
+    """Read *path* with both readers and return ``(cpu_array, gpu_host_array)``.
+
+    Uses the low-level ``read_to_array`` for CPU so that nodata sentinels
+    stay as the literal value (this module checks LERC mask preservation,
+    not the higher-level NaN promotion that ``open_geotiff`` performs).
+
+    The GPU reader (``read_geotiff_gpu``) applies the same nodata masking
+    that ``open_geotiff`` does (PR #1542), so its output uses NaN where
+    the sentinel was. Callers that want a bit-for-bit comparison should
+    pass ``raw_gpu=True`` to skip the high-level masking on the GPU side.
+    """
     from xrspatial.geotiff import read_geotiff_gpu
     from xrspatial.geotiff._reader import read_to_array
 
@@ -90,6 +100,17 @@ def _read_cpu_gpu(path):
     return cpu, gpu_host
 
 
+def _restore_sentinel(arr, nodata):
+    """Replace NaN positions in *arr* with *nodata* so high-level GPU
+    reads compare bit-exactly against low-level CPU reads (which keep
+    the sentinel value verbatim)."""
+    if nodata is None or arr.dtype.kind != 'f' or np.isnan(nodata):
+        return arr
+    out = arr.copy()
+    out[np.isnan(out)] = arr.dtype.type(nodata)
+    return out
+
+
 @_gpu_only
 class TestGpuLercValidMask:
     """End-to-end TIFF round-trips comparing GPU vs CPU output."""
@@ -141,9 +162,15 @@ def invalid_pred(a):
               nodata=-9999.0)
 
         cpu, gpu = _read_cpu_gpu(path)
-        np.testing.assert_array_equal(cpu, gpu)
+        # ``read_geotiff_gpu`` applies the high-level nodata mask (#1542),
+        # so masked pixels come back as NaN. ``read_to_array`` keeps the
+        # sentinel verbatim. Restore the sentinel on the GPU side so the
+        # bit-for-bit comparison still pins LERC mask preservation.
+        gpu_with_sentinel = _restore_sentinel(gpu, -9999.0)
+        np.testing.assert_array_equal(cpu, gpu_with_sentinel)
         for (r, c) in invalid_positions:
-            assert gpu[r, c] == np.float32(-9999.0)
+            assert np.isnan(gpu[r, c])
+            assert gpu_with_sentinel[r, c] == np.float32(-9999.0)
 
     def test_uint16_sentinel_nodata(self, tmp_path, lerc_writer_with_mask):
         """Uint16 LERC + sentinel nodata (65535): GPU matches CPU."""
@@ -164,9 +191,20 @@ def invalid_pred(a):
               nodata=65535)
 
         cpu, gpu = _read_cpu_gpu(path)
-        np.testing.assert_array_equal(cpu, gpu)
+        # ``read_geotiff_gpu`` applies the high-level nodata mask on
+        # integer rasters (#1542): the array is promoted to float64 with
+        # NaN where the sentinel was. ``read_to_array`` keeps uint16 with
+        # the sentinel literal. Restore the sentinel + dtype on the GPU
+        # side so the bit-for-bit comparison still pins LERC mask
+        # preservation. Replace NaN before the uint16 cast to avoid
+        # numpy's "invalid value encountered in cast" warning.
+        assert gpu.dtype == np.float64
+        gpu_no_nan = np.where(np.isnan(gpu), 65535.0, gpu)
+        gpu_u16 = gpu_no_nan.astype(np.uint16)
+        np.testing.assert_array_equal(cpu, gpu_u16)
         for (r, c) in invalid_positions:
-            assert gpu[r, c] == np.uint16(65535)
+            assert np.isnan(gpu[r, c])
+            assert gpu_u16[r, c] == np.uint16(65535)
 
     def test_no_mask_roundtrip_bitexact(self, tmp_path):
         """All-valid LERC (no encoded mask): GPU and CPU agree bit-exact."""

xrspatial/geotiff/tests/test_sparse_cog.py

@@ -105,10 +105,15 @@ def test_sparse_tile_gpu_round_trip(self, tmp_path):
         _write_sparse_tiled(path, nodata=0)
 
         arr = open_geotiff(path, gpu=True)
-        # GPU read returns raw values (no nodata->NaN promotion); sparse
-        # tiles become the nodata sentinel (0 here).
+        # GPU read now applies the high-level nodata mask (#1542): the
+        # source uint16 raster is promoted to float64 and sentinel
+        # values become NaN, matching the CPU eager path.
         host = arr.data.get()
-        assert host.dtype == np.uint16
-        assert host[:64, :64].sum() == 64 * 64 * 100
-        assert host[:64, 64:].sum() == 0
-        assert host[64:, :].sum() == 0
+        assert host.dtype == np.float64
+        assert arr.attrs.get('nodata') == 0.0
+        # Top-left tile holds the real data (100 in every cell).
+        np.testing.assert_array_equal(host[:64, :64],
+                                      np.full((64, 64), 100.0))
+        # Sparse tiles + the bottom row of tiles become NaN via the mask.
+        assert np.all(np.isnan(host[:64, 64:]))
+        assert np.all(np.isnan(host[64:, :]))