diff --git a/docs/source/reference/focal.rst b/docs/source/reference/focal.rst index 648927b46..9cc7563f6 100644 --- a/docs/source/reference/focal.rst +++ b/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 diff --git a/docs/source/user_guide/data_types.rst b/docs/source/user_guide/data_types.rst index 77df9340c..911b3ab78 100644 --- a/docs/source/user_guide/data_types.rst +++ b/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 =================== diff --git a/xrspatial/convolution.py b/xrspatial/convolution.py index acb38d3f2..68cffb134 100644 --- a/xrspatial/convolution.py +++ b/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 """ diff --git a/xrspatial/edge_detection.py b/xrspatial/edge_detection.py index 1215d9364..80d1e11ff 100644 --- a/xrspatial/edge_detection.py +++ b/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) diff --git a/xrspatial/focal.py b/xrspatial/focal.py index 340b3936a..5f932c4ea 100644 --- a/xrspatial/focal.py +++ b/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), diff --git a/xrspatial/tests/test_convolution.py b/xrspatial/tests/test_convolution.py index db0229876..ff3851217 100644 --- a/xrspatial/tests/test_convolution.py +++ b/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): diff --git a/xrspatial/tests/test_focal.py b/xrspatial/tests/test_focal.py index 3670387b6..a886ec120 100644 --- a/xrspatial/tests/test_focal.py +++ b/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