From 98b023f4bf989d4f497f478d30ba495557db050e Mon Sep 17 00:00:00 2001
From: Brendan Collins <brendancol@gmail.com>
Date: Tue, 23 Jun 2026 05:16:49 -0700
Subject: [PATCH 1/2] Preserve input dtype on perlin() dask backends

The dask+numpy and dask+cupy paths hardcoded float32 in their
map_blocks calls, so a float64 input was silently downcast to float32
while the numpy and cupy backends kept float64. Derive the output
dtype from data.dtype in both dask functions and add cross-backend
dtype regression tests covering all four backends.
---
 xrspatial/perlin.py            | 15 +++++++++------
 xrspatial/tests/test_perlin.py | 33 +++++++++++++++++++++++++++++++++
 2 files changed, 42 insertions(+), 6 deletions(-)

diff --git a/xrspatial/perlin.py b/xrspatial/perlin.py
index fe08cff8c..e4593be0f 100644
--- a/xrspatial/perlin.py
+++ b/xrspatial/perlin.py
@@ -2,7 +2,6 @@
 
 # std lib
 import math
-from functools import partial
 
 # 3rd-party
 import numpy as np
@@ -112,6 +111,7 @@ def _perlin_dask_numpy(data: da.Array,
                        freq: tuple,
                        seed: int) -> da.Array:
     p = _make_perm_table(seed)
+    out_dtype = data.dtype
 
     height, width = data.shape
     linx = da.linspace(0, freq[0], width, endpoint=False, dtype=np.float32,
@@ -120,8 +120,10 @@ def _perlin_dask_numpy(data: da.Array,
                        chunks=data.chunks[0][0])
     x, y = da.meshgrid(linx, liny)
 
-    _func = partial(_perlin, p)
-    data = da.map_blocks(_func, x, y, meta=np.array((), dtype=np.float32),
+    def _func(x_blk, y_blk):
+        return _perlin(p, x_blk, y_blk).astype(out_dtype)
+
+    data = da.map_blocks(_func, x, y, meta=np.array((), dtype=out_dtype),
                          **_dask_task_name_kwargs('xrspatial.perlin'))
 
     # persist so min/ptp don't recompute the noise from scratch
@@ -254,6 +256,7 @@ def _perlin_dask_cupy(data: da.Array,
                       freq: tuple,
                       seed: int) -> da.Array:
     p = cupy.asarray(_make_perm_table(seed))
+    out_dtype = data.dtype
 
     height, width = data.shape
 
@@ -266,13 +269,13 @@ def _chunk_perlin(block, block_info=None):
         y0 = freq[1] * y_start / height
         y1 = freq[1] * y_end / height
 
-        out = cupy.empty(block.shape, dtype=cupy.float32)
+        out = cupy.empty(block.shape, dtype=out_dtype)
         griddim, blockdim = cuda_args(block.shape)
         _perlin_gpu[griddim, blockdim](p, x0, x1, y0, y1, 1.0, out)
         return out
 
-    data = da.map_blocks(_chunk_perlin, data, dtype=cupy.float32,
-                         meta=cupy.array((), dtype=cupy.float32),
+    data = da.map_blocks(_chunk_perlin, data, dtype=out_dtype,
+                         meta=cupy.array((), dtype=out_dtype),
                          **_dask_task_name_kwargs('xrspatial.perlin'))
 
     # persist so min/max don't recompute the noise from scratch
diff --git a/xrspatial/tests/test_perlin.py b/xrspatial/tests/test_perlin.py
index 8942e6456..4517641fa 100644
--- a/xrspatial/tests/test_perlin.py
+++ b/xrspatial/tests/test_perlin.py
@@ -118,3 +118,36 @@ def test_perlin_float64_input():
     # Normalized to [0, 1]
     assert result.data.min() >= 0.0
     assert result.data.max() <= 1.0
+
+
+@pytest.mark.parametrize("dtype", [np.float32, np.float64])
+def test_perlin_dask_cpu_preserves_dtype(dtype):
+    # Regression: the dask backends used to hardcode float32, silently
+    # downcasting float64 input while numpy/cupy preserved it.
+    import dask.array as da
+    data = da.from_array(np.zeros((20, 20), dtype=dtype), chunks=(10, 10))
+    raster = xr.DataArray(data, dims=['y', 'x'])
+    result = perlin(raster)
+    assert result.dtype == dtype
+
+
+@cuda_and_cupy_available
+@pytest.mark.parametrize("dtype", [np.float32, np.float64])
+def test_perlin_gpu_preserves_dtype(dtype):
+    import cupy
+    data = cupy.zeros((20, 20), dtype=dtype)
+    raster = xr.DataArray(data, dims=['y', 'x'])
+    result = perlin(raster)
+    assert result.dtype == dtype
+
+
+@cuda_and_cupy_available
+@dask_array_available
+@pytest.mark.parametrize("dtype", [np.float32, np.float64])
+def test_perlin_dask_gpu_preserves_dtype(dtype):
+    import cupy
+    import dask.array as da
+    data = da.from_array(cupy.zeros((20, 20), dtype=dtype), chunks=(10, 10))
+    raster = xr.DataArray(data, dims=['y', 'x'])
+    result = perlin(raster)
+    assert result.dtype == dtype

From 1f2254f247f2366c78ea446dcc49348f176229e8 Mon Sep 17 00:00:00 2001
From: Brendan Collins <brendancol@gmail.com>
Date: Tue, 23 Jun 2026 05:20:26 -0700
Subject: [PATCH 2/2] Address review nit: assert dask float64 perlin is finite
 and normalized

---
 xrspatial/tests/test_perlin.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/xrspatial/tests/test_perlin.py b/xrspatial/tests/test_perlin.py
index 4517641fa..81b959c59 100644
--- a/xrspatial/tests/test_perlin.py
+++ b/xrspatial/tests/test_perlin.py
@@ -129,6 +129,10 @@ def test_perlin_dask_cpu_preserves_dtype(dtype):
     raster = xr.DataArray(data, dims=['y', 'x'])
     result = perlin(raster)
     assert result.dtype == dtype
+    computed = result.data.compute()
+    assert np.isfinite(computed).all()
+    assert computed.min() >= 0.0
+    assert computed.max() <= 1.0
 
 
 @cuda_and_cupy_available