Promote 32/64-bit integer rasters to float64 in _promote_float (#3684)

docs/source/reference/focal.rst

@@ -12,8 +12,9 @@ Focal
 
 .. note::
 
-   ``mean`` and ``focal_stats`` output **float32**; ``hotspots`` outputs
-   **int8** confidence bands. Both ``mean`` and ``hotspots`` exclude NaN
+   ``mean`` and ``focal_stats`` keep a float input's dtype and promote
+   integers to float32 (bool, 8/16-bit) or float64 (32/64-bit);
+   ``hotspots`` outputs **int8** confidence bands. Both ``mean`` and ``hotspots`` exclude NaN
    neighbours from their neighbourhood computation.
 
 Apply

docs/source/user_guide/data_types.rst

@@ -9,7 +9,7 @@ This guide explains how xarray-spatial handles floating-point data types (float3
 Overview
 ========
 
-xarray-spatial standardizes on **float32** (32-bit floating-point) for output data in most analytical functions. This design decision provides a balance between computational precision, memory efficiency, and performance that is well-suited for raster analysis tasks.
+xarray-spatial standardizes on **float32** (32-bit floating-point) for output data in most analytical functions (the focal and convolution functions are the exception for 32/64-bit integer input; see below). This design decision provides a balance between computational precision, memory efficiency, and performance that is well-suited for raster analysis tasks.
 
 .. warning::
 
@@ -149,11 +149,16 @@ Focal statistics and convolution operations also use float32:
    # Integer input data
    data = xr.DataArray(np.random.randint(0, 255, (100, 100), dtype=np.uint8))
 
-   # Focal operations convert to and output float32
+   # Focal operations convert 8/16-bit integers to float32
    kernel = circle_kernel(1, 1, 3)
    smoothed = convolve_2d(data, kernel)
    print(f"Convolution output dtype: {smoothed.dtype}")  # float32
 
+   # 32/64-bit integers are widened to float64 instead, because float32
+   # cannot represent integers above 2**24 exactly
+   wide = xr.DataArray(np.arange(10_000, dtype=np.int32).reshape(100, 100))
+   print(f"Convolution output dtype: {convolve_2d(wide, kernel).dtype}")  # float64
+
 
 Backend Consistency
 ===================

xrspatial/convolution.py

@@ -19,9 +19,19 @@ class cupy(object):
 
 
 def _promote_float(dtype):
-    """Return at least float32; preserve float64."""
-    if np.issubdtype(dtype, np.floating):
+    """Return a float dtype that represents ``dtype`` values exactly.
+
+    Float inputs keep their dtype. 32- and 64-bit integers go to float64
+    (exact for int32/uint32, and up to 2**53 for int64/uint64); bool and
+    8/16-bit integers go to float32, which holds them exactly. float32
+    cannot separate integers above 2**24, so casting wide ints to it
+    silently corrupted large-valued rasters (#3684).
+    """
+    dtype = np.dtype(dtype)
+    if dtype.kind == 'f':
         return dtype
+    if dtype.kind in 'iu' and dtype.itemsize > 2:
+        return np.float64
     return np.float32
 
 
@@ -563,8 +573,8 @@ def convolve_2d(data, kernel, boundary='nan'):
     Parameters
     ----------
     data : numpy.ndarray, cupy.ndarray, or dask.array.Array
-        2D array of values to convolve. Integer inputs are promoted to at
-        least float32.
+        2D array of values to convolve. Bool and 8/16-bit integer inputs
+        are promoted to float32, 32/64-bit integers to float64 (#3684).
     kernel : numpy.ndarray
         Impulse kernel with an odd number of rows and columns.
     boundary : str, default='nan'
@@ -665,7 +675,7 @@ def convolution_2d(agg, kernel, name='convolution_2d', boundary='nan'):
         array([[nan, nan, nan, nan, nan, nan],
                [nan, 35., 40., 45., 50., nan],
                [nan, 65., 70., 75., 80., nan],
-               [nan, nan, nan, nan, nan, nan]], dtype=float32)
+               [nan, nan, nan, nan, nan, nan]])
         Dimensions without coordinates: dim_0, dim_1
 
     convolution_2d() works with Dask with NumPy backed DataArray.
@@ -723,7 +733,7 @@ def convolution_2d(agg, kernel, name='convolution_2d', boundary='nan'):
         array([[ nan,  nan,  nan,  nan,  nan,  nan],
                [ nan,  50.,  58.,  66.,  74.,  nan],
                [ nan,  98., 106., 114., 122.,  nan],
-               [ nan,  nan,  nan,  nan,  nan,  nan]], dtype=float32)
+               [ nan,  nan,  nan,  nan,  nan,  nan]])
         Dimensions without coordinates: dim_0, dim_1
     """
 

xrspatial/edge_detection.py

@@ -31,14 +31,14 @@
 def _promote_wide_int(data):
     """Cast 32/64-bit integer arrays to float64 before convolution.
 
-    ``convolve_2d`` promotes integer inputs to float32, whose 24-bit
-    mantissa cannot represent integers above 2**24: unit steps between
-    large values vanish in the cast and gradients silently collapse to
-    zero (#3680). float64 is exact for every int32/uint32 value and for
-    int64/uint64 up to 2**53. 8- and 16-bit integers are exactly
-    representable in float32, so they keep ``convolve_2d``'s promotion.
-    Works on numpy, cupy, and dask arrays alike (``astype`` is lazy on
-    dask).
+    ``convolve_2d`` used to promote every integer input to float32, whose
+    24-bit mantissa cannot represent integers above 2**24: unit steps
+    between large values vanished in the cast and gradients silently
+    collapsed to zero (#3680). Since #3684 ``convolve_2d``'s own
+    ``_promote_float`` applies the same split (8/16-bit ints to float32,
+    wider ints to float64), so this pre-cast is redundant and kept only
+    so the edge detectors do not depend on that helper's rule. Works on
+    numpy, cupy, and dask arrays alike (``astype`` is lazy on dask).
     """
     if data.dtype.kind in 'iu' and data.dtype.itemsize > 2:
         return data.astype(np.float64)

xrspatial/focal.py

@@ -101,9 +101,10 @@ def _check_kernel_vs_raster_memory(kernel, rows, cols, func_name, chunks=None,
 
     ``itemsize`` is the byte width of the dtype the focal internals
     compute in.  Since #2805 ``apply`` and ``focal_stats`` preserve
-    float64 input (``_promote_float``), so callers must pass 8 for
-    float64 rasters; ``hotspots`` always computes in float32 and uses
-    the default of 4 (issue #3223).
+    float64 input, and since #3684 32/64-bit integer input is widened to
+    float64 too (``_promote_float``), so callers must pass 8 for those
+    rasters; ``hotspots`` always computes in float32 and uses the
+    default of 4 (issue #3223).
     """
     krows, kcols = kernel.shape
     pad_h = krows // 2
@@ -372,7 +373,8 @@ def mean(agg, passes=1, excludes=None, name='mean', boundary='nan'):
         for each data variable.
         2D aggregate array of filtered values.
         The input float dtype is preserved (float64 in, float64 out);
-        integer inputs are promoted to float32, matching ``apply`` and
+        bool and 8/16-bit integer inputs are promoted to float32 and
+        32/64-bit integers to float64, matching ``apply`` and
         ``focal_stats``.
 
     Examples
@@ -466,8 +468,9 @@ def mean(agg, passes=1, excludes=None, name='mean', boundary='nan'):
                          coords=agg.coords,
                          attrs=agg.attrs)
 
-    # Preserve the input float dtype, promoting ints to float32 -- the same
-    # contract as apply() / focal_stats() (#2769) and convolve_2d() (#1096).
+    # Preserve the input float dtype, promoting narrow ints to float32 and
+    # wide ints to float64 (#3684) -- the same contract as apply() /
+    # focal_stats() (#2769) and convolve_2d() (#1096).
     # Cast excludes to the working dtype so value matching behaves the same
     # on every backend (issue #3214).
     fdtype = _promote_float(agg.data.dtype)
@@ -780,7 +783,8 @@ def apply(agg=None, kernel=None, func=None, name='focal_apply',
 
     rows, cols = agg.shape[-2], agg.shape[-1]
     # Budget for the dtype the internals will actually compute in:
-    # 8 bytes/cell for float64 input, 4 otherwise (issue #3223).
+    # 8 bytes/cell for float64 and 32/64-bit int input, 4 otherwise
+    # (issues #3223, #3684).
     itemsize = np.dtype(_promote_float(agg.dtype)).itemsize
     _check_kernel_vs_raster_memory(kernel, rows, cols, func_name='apply',
                                    chunks=getattr(agg.data, 'chunks', None),
@@ -1433,7 +1437,8 @@ def focal_stats(agg,
 
     rows, cols = agg.shape[-2], agg.shape[-1]
     # Budget for the dtype the internals will actually compute in:
-    # 8 bytes/cell for float64 input, 4 otherwise (issue #3223).
+    # 8 bytes/cell for float64 and 32/64-bit int input, 4 otherwise
+    # (issues #3223, #3684).
     itemsize = np.dtype(_promote_float(agg.dtype)).itemsize
     _check_kernel_vs_raster_memory(kernel, rows, cols, func_name='focal_stats',
                                    chunks=getattr(agg.data, 'chunks', None),

xrspatial/tests/test_convolution.py

@@ -10,7 +10,8 @@
 from xrspatial.tests.general_checks import (assert_boundary_mode_correctness,
                                             assert_numpy_equals_cupy, assert_numpy_equals_dask_cupy,
                                             assert_numpy_equals_dask_numpy, create_test_raster,
-                                            cuda_and_cupy_available, dask_array_available)
+                                            cuda_and_cupy_available, dask_array_available,
+                                            has_cuda_and_cupy)
 from xrspatial.utils import VALID_BOUNDARY_MODES
 
 KERNEL = circle_kernel(1, 1, 1)
@@ -30,6 +31,50 @@ def test_convolve_2d_dask_dtype_matches_numpy(dtype):
     assert lazy.compute().dtype == eager.dtype
 
 
+@pytest.mark.parametrize("dtype, expected", [
+    (np.uint8, np.float32),
+    (np.int16, np.float32),
+    (np.int32, np.float64),
+    (np.uint32, np.float64),
+    (np.int64, np.float64),
+    (np.float32, np.float32),
+    (np.float64, np.float64),
+])
+def test_convolve_2d_promotion_rule_3684(dtype, expected):
+    # Floats keep their dtype; narrow ints go to float32, which holds them
+    # exactly; 32/64-bit ints go to float64 because float32 cannot separate
+    # integers above 2**24 (#3684).
+    data = np.arange(64).reshape(8, 8).astype(dtype)
+    assert convolve_2d(data, KERNEL).dtype == expected
+
+
+@pytest.mark.parametrize("backend", ['numpy', 'dask+numpy', 'cupy', 'dask+cupy'])
+@pytest.mark.parametrize("dtype", [np.int32, np.int64])
+def test_convolution_2d_wide_int_exact_above_2_24_3684(backend, dtype):
+    # Regression for #3684: a horizontal-gradient kernel on an int32 raster
+    # offset by 1e8 returned 0 because the float32 cast erased the unit
+    # steps. Every backend must match the float64 reference exactly.
+    if 'cupy' in backend and not has_cuda_and_cupy():
+        pytest.skip("Requires CUDA and CuPy")
+
+    base = 100_000_000
+    data = (np.arange(30).reshape(5, 6) % 6 + base).astype(dtype)
+    kernel = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]], dtype=np.float64)
+    expected = convolution_2d(create_test_raster(data.astype(np.float64)), kernel)
+    assert expected.data[2, 2] == 8.0
+
+    # chunks only applies to the dask backends; (2, 3) splits the 5x6 raster
+    agg = create_test_raster(data, backend=backend, chunks=(2, 3))
+    out = convolution_2d(agg, kernel).data
+    if hasattr(out, 'compute'):
+        out = out.compute()
+    if hasattr(out, 'get'):
+        out = out.get()
+    assert out.dtype == np.float64
+    # exact, not allclose: at 1e8 a relative tolerance would hide the bug
+    np.testing.assert_array_equal(out, expected.data)
+
+
 @cuda_and_cupy_available
 @pytest.mark.parametrize("dtype", [np.int32, np.float32, np.float64])
 def test_convolve_2d_cupy_dtype_matches_numpy(dtype):

xrspatial/tests/test_focal.py

@@ -10,8 +10,8 @@
 import xarray as xr
 
 from xrspatial import mean
-from xrspatial.convolution import (annulus_kernel, calc_cellsize, circle_kernel, convolution_2d,
-                                   convolve_2d, custom_kernel)
+from xrspatial.convolution import (_promote_float, annulus_kernel, calc_cellsize, circle_kernel,
+                                   convolution_2d, convolve_2d, custom_kernel)
 from xrspatial.focal import apply, focal_stats, hotspots
 from xrspatial.tests.general_checks import (assert_boundary_mode_correctness,
                                             create_test_raster,
@@ -956,7 +956,8 @@ def test_apply_keeps_float32(backend):
 # and to float32 on cupy/dask+cupy, so the output dtype depended on the
 # backend and float64 rasters lost precision on GPU. It now follows the same
 # _promote_float contract as apply()/focal_stats() (#2769) and convolve_2d()
-# (#1096): float dtypes are preserved, ints promote to float32.
+# (#1096): float dtypes are preserved, narrow ints promote to float32 and
+# wide ints to float64 (#3684).
 
 
 @pytest.mark.parametrize("backend", ['numpy', 'cupy', 'dask+numpy', 'dask+cupy'])
@@ -988,15 +989,50 @@ def test_mean_keeps_float32_3214(backend):
 
 
 @pytest.mark.parametrize("backend", ['numpy', 'cupy'])
-def test_mean_int_promotes_to_float32_3214(backend):
+@pytest.mark.parametrize("dtype, expected", [
+    (np.uint8, np.float32),
+    (np.int16, np.float32),
+    (np.int32, np.float64),
+    (np.uint32, np.float64),
+    (np.int64, np.float64),
+])
+def test_mean_int_promotion_3214(backend, dtype, expected):
+    # #3214 fixed ints promoting to float32 on every backend; #3684 widened
+    # 32/64-bit ints to float64 because float32 cannot hold them exactly.
     from xrspatial.tests.general_checks import has_cuda_and_cupy
     if 'cupy' in backend and not has_cuda_and_cupy():
         pytest.skip("Requires CUDA and CuPy")
 
-    data = np.arange(20, dtype=np.int32).reshape(4, 5)
+    data = np.arange(20).reshape(4, 5).astype(dtype)
     agg = create_test_raster(data, backend=backend)
     result = mean(agg)
-    assert _compute_dtype(result) == np.float32
+    assert _compute_dtype(result) == expected
+
+
+@pytest.mark.parametrize("backend", ['numpy', 'cupy', 'dask+numpy', 'dask+cupy'])
+@pytest.mark.parametrize("dtype", [np.int32, np.int64])
+def test_mean_wide_int_exact_above_2_24_3684(backend, dtype):
+    # Regression for #3684: int32/int64 rasters with values above 2**24
+    # were cast to float32 before the kernel ran, so the unit steps between
+    # neighbours vanished and mean() returned the (rounded) base value.
+    from xrspatial.tests.general_checks import has_cuda_and_cupy
+    if 'cupy' in backend and not has_cuda_and_cupy():
+        pytest.skip("Requires CUDA and CuPy")
+
+    base = 100_000_000
+    data = (np.arange(30).reshape(5, 6) % 6 + base).astype(dtype)
+    # chunks only applies to the dask backends; (2, 3) splits the 5x6 raster
+    agg = create_test_raster(data, backend=backend, chunks=(2, 3))
+    result = mean(agg)
+    expected = mean(create_test_raster(data.astype(np.float64)))
+    out = result.data
+    if hasattr(out, 'compute'):
+        out = out.compute()
+    if hasattr(out, 'get'):
+        out = out.get()
+    assert out.dtype == np.float64
+    np.testing.assert_array_equal(out, expected.data)
+    assert out[2, 2] == base + 2
 
 
 @cuda_and_cupy_available
@@ -2270,7 +2306,7 @@ def test_mean_dtype_consistent_across_backends_3217(backend, in_dtype):
     result = mean(agg)
     # mean() promotes via _promote_float before dispatch (#3214); every
     # backend must keep that dtype, lazily and computed.
-    expected = np.float64 if in_dtype == np.float64 else np.float32
+    expected = _promote_float(in_dtype)
     assert result.dtype == expected
     assert _computed_dtype(result) == expected
 
@@ -2303,7 +2339,7 @@ def test_apply_dask_advertised_dtype_matches_computed_3217(backend, in_dtype):
         result = apply(agg, kernel, func=_focal_mean_cuda)
     else:
         result = apply(agg, kernel)
-    expected = np.float64 if in_dtype == np.float64 else np.float32
+    expected = _promote_float(in_dtype)
     assert result.dtype == expected
     assert _computed_dtype(result) == expected
 
@@ -2321,7 +2357,7 @@ def test_focal_stats_dask_advertised_dtype_matches_computed_3217(backend, in_dty
     kernel = custom_kernel(np.array([[0, 1, 0], [1, 1, 1], [0, 1, 0]]))
     agg = create_test_raster(data, backend=backend, chunks=(2, 2))
     result = focal_stats(agg, kernel, stats_funcs=['mean', 'std'])
-    expected = np.float64 if in_dtype == np.float64 else np.float32
+    expected = _promote_float(in_dtype)
     assert result.dtype == expected
     assert _computed_dtype(result) == expected