Refuse pack=True values the packed integer dtype cannot hold (#3272)

* Refuse pack=True values the packed integer dtype cannot hold (#3260)

The integer restore in _pack cast unguarded: finite values whose packed
form fell outside the dtype range wrapped in the astype (40000 -> -25536
in int16) and +/-Inf cast to a platform-defined integer, all silently.
A range/finiteness guard now runs after the round and before the cast,
raising ValueError. numpy and cupy buffers raise at call time; dask
backings map the guard into the graph so it raises from the write's
single compute, same shape as the #3235 NaN guard. The exclusive upper
bound is iinfo.max + 1 (a power of two, exact in float64) so 64-bit
dtypes reject exactly-2**63 instead of letting it wrap.

Also routes the eager NaN-no-sentinel scan through _pack_guard_no_nan:
the old bool(out.isnull().any()) materialised via np.asarray and raised
TypeError on cupy-backed eager input before the scan could answer.

* Record accuracy-sweep pass 27 for geotiff (#3260)

* Address review: assert on the packed output file, drop redundant attrs copy (#3260)

Author: brendancol

.claude/sweep-accuracy-state.csv

@@ -1730,6 +1730,44 @@ def _pack_fill_nan(chunk, fill):
     return out
 
 
+def _pack_guard_int_range(chunk, tgt_name, lo, hi_excl):
+    """Per-chunk range / finiteness guard for ``_pack``'s integer restore.
+
+    Runs after the round and before the ``astype``. Raises ``ValueError``
+    when the chunk holds non-finite values or values outside the target
+    integer dtype's range -- the cast would silently wrap them into
+    valid-looking integers (issue #3260). Returns the chunk unchanged
+    when clean.
+
+    ``lo`` is ``iinfo.min`` and ``hi_excl`` is ``iinfo.max + 1``, both as
+    floats. The exclusive upper bound matters for the 64-bit dtypes:
+    ``float(iinfo.max)`` rounds up to ``2**63`` / ``2**64``, so a
+    ``chunk > float(iinfo.max)`` test would pass a value that still wraps
+    in the cast. ``iinfo.max + 1`` is a power of two for every integer
+    dtype and therefore exact in float64. ``iinfo.min`` is likewise a
+    power of two (or zero) and exact.
+
+    Mapped over dask blocks so the scan fuses with the write's single
+    compute (same shape as :func:`_pack_guard_no_nan`); numpy and cupy
+    buffers are checked eagerly at ``_pack`` call time. ``np.isfinite``
+    and the comparisons dispatch on the array type, so the one function
+    covers all four backends.
+    """
+    if chunk.dtype.kind == 'f':
+        bad = ~np.isfinite(chunk)
+        bad |= chunk < lo
+        bad |= chunk >= hi_excl
+        if bool(bad.any()):
+            raise ValueError(
+                f"pack=True: data contains values that the packed "
+                f"integer dtype {tgt_name} cannot represent (valid "
+                f"range {int(lo)}..{int(hi_excl) - 1} after the "
+                f"scale/offset reversal) or that are not finite; the "
+                f"cast would silently wrap them. Clip the data to the "
+                f"packed range or rescale before writing.")
+    return chunk
+
+
 def _pack_guard_no_nan(chunk, tgt_name):
     """Per-chunk NaN guard for ``_pack``'s no-sentinel integer restore.
 
@@ -1785,12 +1823,16 @@ def _pack(data, nodata=None):
 
     Raises ``ValueError`` when ``data`` carries no unpack state to
     reverse, when an integer dtype must be restored but NaN pixels are
-    present with no declared sentinel to fill them, or when the ``nodata``
+    present with no declared sentinel to fill them, when the ``nodata``
     kwarg cannot be represented in the packed integer dtype (out of range
     or fractional) -- the cast would wrap or round the fill value away
-    from what the GDAL_NODATA tag declares.
+    from what the GDAL_NODATA tag declares -- or when any pixel's packed
+    value is non-finite or falls outside the packed integer dtype's
+    range, where the cast would silently wrap it into a valid-looking
+    integer (issue #3260).
 
-    Error timing for the NaN-with-no-sentinel case depends on the
+    Error timing for the NaN-with-no-sentinel and packed-range cases
+    depends on the
     backing: numpy-backed data raises here, at call time. dask-backed
     data defers the scan into the graph (one per-chunk check fused with
     the write's single compute) so the upstream graph is not executed a
@@ -1887,14 +1929,36 @@ def _pack(data, nodata=None):
             out = out.copy(data=out.data.map_blocks(
                 _pack_guard_no_nan, tgt.name,
                 dtype=out.dtype, meta=out.data._meta))
-        elif bool(out.isnull().any()):
-            raise ValueError(
-                f"pack=True: cannot restore integer dtype {tgt.name}: "
-                "NaN pixels are present but no nodata sentinel is "
-                "declared to fill them.")
+        else:
+            # Buffer-level scan, not ``bool(out.isnull().any())``: the
+            # xarray reduction materialises through ``np.asarray`` and
+            # raised TypeError on cupy-backed eager input, crashing
+            # every no-sentinel integer pack on the gpu backend before
+            # the scan could answer (#3260). ``_pack_guard_no_nan``
+            # dispatches ``np.isnan`` on the array type, so numpy and
+            # cupy take the same path here.
+            _pack_guard_no_nan(out.data, tgt.name)
 
     if tgt.kind in ('i', 'u'):
         out = out.round()
+        # Out-of-range and non-finite packed values would silently wrap
+        # in the astype below, exactly the corruption the NaN guard
+        # above exists to prevent: a finite value whose packed form
+        # exceeds the dtype range wraps (40000 -> -25536 in int16), and
+        # +/-Inf casts to a platform-defined integer (issue #3260).
+        # Runs after the round so a value like 32767.4 that rounds back
+        # into range is not rejected. Same eager/lazy split as the NaN
+        # guard: numpy / cupy raise here, dask raises from the write's
+        # single compute.
+        info = np.iinfo(tgt)
+        lo = float(info.min)
+        hi_excl = float(int(info.max) + 1)
+        if hasattr(out.data, 'dask'):
+            out = out.copy(data=out.data.map_blocks(
+                _pack_guard_int_range, tgt.name, lo, hi_excl,
+                dtype=out.dtype, meta=out.data._meta))
+        else:
+            _pack_guard_int_range(out.data, tgt.name, lo, hi_excl)
 
     out = out.astype(tgt)
 

xrspatial/geotiff/_attrs.py

@@ -390,7 +390,11 @@ def to_geotiff(data: xr.DataArray | np.ndarray,
         integer dtype must be restored, the ``ValueError`` is raised at
         call time for numpy-backed data; for dask-backed data the check
         runs per chunk inside the write's compute (issue #3235), so the
-        error surfaces during the write. The write itself stays atomic
+        error surfaces during the write. The same timing applies to
+        pixels whose packed value is non-finite or falls outside the
+        packed integer dtype's range: the cast would silently wrap
+        them, so the write is refused with ``ValueError`` instead
+        (issue #3260). The write itself stays atomic
         (temp file plus rename), so no partial output is left at the
         destination path.
 

xrspatial/geotiff/_writers/eager.py

@@ -0,0 +1,195 @@
+"""``to_geotiff(pack=True)`` refuses packed values the dtype cannot hold (#3260).
+
+``_pack`` reverses the unpack transform and casts to the integer dtype
+recorded on ``attrs['mask_and_scale_dtype']``. Before #3260 the cast ran
+unguarded: a finite value whose packed form fell outside the dtype range
+wrapped (40000 -> -25536 in int16) and +/-Inf cast to a platform-defined
+integer, both silently. The guard raises ``ValueError`` instead -- at
+call time for numpy / cupy buffers, from the write's single compute for
+dask backings (same timing as the #3235 NaN guard).
+"""
+import numpy as np
+import pytest
+import xarray as xr
+
+from xrspatial.geotiff import open_geotiff, to_geotiff
+from xrspatial.geotiff._attrs import _pack_guard_int_range
+
+from .._helpers.markers import requires_gpu
+
+
+def _write_scaled_int16(path, *, nodata=None):
+    """int16 source with SCALE=0.1 so unpacked values are data * 0.1."""
+    data = np.array([[100, 200], [300, 32000]], dtype=np.int16)
+    attrs = {
+        "crs": 4326,
+        "gdal_metadata": {"SCALE": "0.1", "OFFSET": "0.0"},
+    }
+    if nodata is not None:
+        attrs["nodata"] = nodata
+    da = xr.DataArray(
+        data,
+        dims=("y", "x"),
+        coords={"y": [1.5, 0.5], "x": [0.5, 1.5]},
+        attrs=attrs,
+    )
+    to_geotiff(da, str(path), nodata=nodata)
+    return str(path)
+
+
+def _unpacked(path, *, gpu=False):
+    kwargs = {"unpack": True}
+    if gpu:
+        kwargs["gpu"] = True
+    return open_geotiff(path, **kwargs).copy()
+
+
+# ---------------------------------------------------------------------------
+# Finite overflow wraps in the cast: refused on every backend
+# ---------------------------------------------------------------------------
+
+
+@pytest.mark.parametrize("chunks,gpu", [
+    pytest.param(None, False, id="numpy"),
+    pytest.param(1, False, id="dask"),
+    pytest.param(None, True, marks=requires_gpu, id="gpu"),
+    pytest.param(1, True, marks=requires_gpu, id="dask-gpu"),
+])
+def test_pack_rejects_finite_overflow(tmp_path, chunks, gpu):
+    src = _write_scaled_int16(tmp_path / "src_overflow_3260.tif")
+    mod = _unpacked(src, gpu=gpu)
+    # Packs to 40000 > int16 max 32767; pre-#3260 this wrapped to -25536.
+    mod.data[0, 0] = 4000.0
+    if chunks is not None:
+        mod = mod.chunk({"y": chunks})
+    out = str(tmp_path / "out_overflow_3260.tif")
+    with pytest.raises(ValueError, match="cannot represent"):
+        to_geotiff(mod, out, pack=True)
+
+
+@pytest.mark.parametrize("chunks", [None, 1], ids=["numpy", "dask"])
+def test_pack_rejects_underflow_unsigned(tmp_path, chunks):
+    """A negative packed value must not wrap into the top of a uint range."""
+    data = np.array([[10, 20], [30, 40]], dtype=np.uint16)
+    da = xr.DataArray(
+        data,
+        dims=("y", "x"),
+        coords={"y": [1.5, 0.5], "x": [0.5, 1.5]},
+        attrs={"crs": 4326,
+               "gdal_metadata": {"SCALE": "0.5", "OFFSET": "0.0"}},
+    )
+    src = str(tmp_path / "src_u16_underflow_3260.tif")
+    to_geotiff(da, src)
+
+    mod = _unpacked(src)
+    mod.data[0, 0] = -1.0  # packs to -2, below uint16 min
+    if chunks is not None:
+        mod = mod.chunk({"y": chunks})
+    out = str(tmp_path / "out_u16_underflow_3260.tif")
+    with pytest.raises(ValueError, match="cannot represent"):
+        to_geotiff(mod, out, pack=True)
+
+
+# ---------------------------------------------------------------------------
+# +/-Inf passes the NaN fill / NaN guard but must not reach the cast
+# ---------------------------------------------------------------------------
+
+
+@pytest.mark.parametrize("chunks", [None, 1], ids=["numpy", "dask"])
+@pytest.mark.parametrize("value", [np.inf, -np.inf], ids=["inf", "neg-inf"])
+def test_pack_rejects_inf(tmp_path, chunks, value):
+    # A declared sentinel routes Inf past the NaN fill (isnan only), so
+    # this leg pins the path where Inf used to reach the cast directly.
+    src = _write_scaled_int16(tmp_path / "src_inf_3260.tif", nodata=-32768)
+    mod = _unpacked(src)
+    mod.data[0, 0] = value
+    if chunks is not None:
+        mod = mod.chunk({"y": chunks})
+    out = str(tmp_path / "out_inf_3260.tif")
+    with pytest.raises(ValueError, match="not finite|cannot represent"):
+        to_geotiff(mod, out, pack=True)
+
+
+# ---------------------------------------------------------------------------
+# No false positives: boundary and round-back-into-range values pass
+# ---------------------------------------------------------------------------
+
+
+@pytest.mark.parametrize("chunks", [None, 1], ids=["numpy", "dask"])
+def test_pack_accepts_full_dtype_range(tmp_path, chunks):
+    """Exact iinfo.min / iinfo.max packed values are representable."""
+    data = np.array([[-32768, 0], [100, 32767]], dtype=np.int16)
+    da = xr.DataArray(
+        data,
+        dims=("y", "x"),
+        coords={"y": [1.5, 0.5], "x": [0.5, 1.5]},
+        attrs={"crs": 4326,
+               "gdal_metadata": {"SCALE": "0.1", "OFFSET": "0.0"}},
+    )
+    src = str(tmp_path / "src_bounds_3260.tif")
+    to_geotiff(da, src)
+
+    mod = _unpacked(src)
+    if chunks is not None:
+        mod = mod.chunk({"y": chunks})
+    out = str(tmp_path / "out_bounds_3260.tif")
+    to_geotiff(mod, out, pack=True)
+
+    back = open_geotiff(out)
+    assert str(back.dtype) == "int16"
+    np.testing.assert_array_equal(np.asarray(back.data), data)
+
+
+def test_pack_accepts_value_that_rounds_back_into_range(tmp_path):
+    """The guard runs after the round: 3276.74 packs to 32767.4 which
+    rounds to 32767 and fits."""
+    src = _write_scaled_int16(tmp_path / "src_round_3260.tif")
+    mod = _unpacked(src)
+    mod.data[0, 0] = 3276.74
+    out = str(tmp_path / "out_round_3260.tif")
+    to_geotiff(mod, out, pack=True)
+    back = open_geotiff(out)
+    assert np.asarray(back.data)[0, 0] == 32767
+
+
+def test_pack_float_target_not_range_guarded(tmp_path):
+    """Float packed dtypes have no wrap problem; large values pass."""
+    data = np.array([[1.0, 2.0], [3.0, -9999.0]], dtype=np.float32)
+    da = xr.DataArray(
+        data,
+        dims=("y", "x"),
+        coords={"y": [1.5, 0.5], "x": [0.5, 1.5]},
+        attrs={"crs": 4326, "nodata": -9999.0,
+               "gdal_metadata": {"SCALE": "2.0", "OFFSET": "0.0"}},
+    )
+    src = str(tmp_path / "src_float_3260.tif")
+    to_geotiff(da, src, nodata=-9999.0)
+
+    mod = _unpacked(src)
+    mod.data[0, 0] = 1e30
+    out = str(tmp_path / "out_float_3260.tif")
+    to_geotiff(mod, out, pack=True)
+    back = open_geotiff(out)
+    assert str(back.dtype) == "float32"
+    # Packed value is 1e30 / SCALE = 5e29, well inside float32 range.
+    assert np.asarray(back.data)[0, 0] == np.float32(5e29)
+
+
+# ---------------------------------------------------------------------------
+# Unit: the 64-bit exclusive upper bound is exact
+# ---------------------------------------------------------------------------
+
+
+def test_guard_rejects_two_pow_63_for_int64():
+    """float64 cannot hold INT64_MAX; float(iinfo.max) rounds up to 2**63.
+    The guard's exclusive bound must reject exactly-2**63, which an
+    inclusive ``> float(iinfo.max)`` test would let through to wrap."""
+    info = np.iinfo(np.int64)
+    chunk = np.array([float(2 ** 63)])
+    with pytest.raises(ValueError, match="cannot represent"):
+        _pack_guard_int_range(
+            chunk, "int64", float(info.min), float(int(info.max) + 1))
+    # One ULP below 2**63 is representable in int64 and passes.
+    ok = np.array([np.nextafter(float(2 ** 63), 0.0)])
+    _pack_guard_int_range(
+        ok, "int64", float(info.min), float(int(info.max) + 1))