From 74879a1ab3884711ffe725e22619a5e2eb32c2a5 Mon Sep 17 00:00:00 2001
From: Brendan Collins <brendancol@gmail.com>
Date: Thu, 25 Jun 2026 11:30:50 -0700
Subject: [PATCH] Fix flood dask backends emitting float32 output for float32
 input

flood_depth() and curve_number_runoff() document a float64 output and
force float64 on the numpy/cupy backends. The dask and dask+cupy
backends ran directly on the input dtype, so a float32 input produced a
float32 result and the output dtype depended on the active backend.

Cast the input to float64 in _flood_depth_dask and _cn_runoff_dask to
match the numpy/cupy backends and the documented contract. Adds float32
input dtype tests across all four backends.

Fixes #3496
---
 xrspatial/flood.py            |  3 ++
 xrspatial/tests/test_flood.py | 69 +++++++++++++++++++++++++++++++++++
 2 files changed, 72 insertions(+)

diff --git a/xrspatial/flood.py b/xrspatial/flood.py
index 27e6f6f3a..94ea3530f 100644
--- a/xrspatial/flood.py
+++ b/xrspatial/flood.py
@@ -176,6 +176,7 @@ def _flood_depth_cupy(hand, wl):
 
 def _flood_depth_dask(hand, wl):
     import dask.array as _da
+    hand = hand.astype(np.float64)
     depth = _da.where(hand <= wl, wl - hand, np.nan)
     depth = _da.where(_da.isnan(hand), np.nan, depth)
     return depth
@@ -366,6 +367,8 @@ def _cn_runoff_cupy(p, cn):
 
 def _cn_runoff_dask(p, cn):
     import dask.array as _da
+    p = p.astype(np.float64)
+    cn = cn.astype(np.float64) if hasattr(cn, 'astype') else float(cn)
     s = (25400.0 / cn) - 254.0
     ia = 0.2 * s
     q = _da.where(p > ia, (p - ia) ** 2 / (p + 0.8 * s), 0.0)
diff --git a/xrspatial/tests/test_flood.py b/xrspatial/tests/test_flood.py
index fa6973979..2ac5ba643 100644
--- a/xrspatial/tests/test_flood.py
+++ b/xrspatial/tests/test_flood.py
@@ -1049,3 +1049,72 @@ def test_flood_depth_vegetation_rejects_negative_mannings_n_dataarray(self):
         bad_n = self._good_raster(-0.1)
         with pytest.raises(ValueError, match="mannings_n"):
             flood_depth_vegetation(hand, slope, bad_n, unit_discharge=1.0)
+
+
+# ===================================================================
+# Output dtype is float64 across all backends (issue #3496)
+# ===================================================================
+
+class TestFloat64OutputAcrossBackends:
+    """float32 input must still produce float64 output on every backend.
+
+    The numpy/cupy backends force float64; the dask backends previously
+    inherited the float32 input dtype, so output dtype depended on the
+    backend.
+    """
+
+    def test_flood_depth_numpy_float32_input(self):
+        hand_data = np.array([[0.0, 2.0], [5.0, 10.0]], dtype=np.float32)
+        hand = create_test_raster(hand_data, backend='numpy', name='hand')
+        assert flood_depth(hand, water_level=5.0).dtype == np.float64
+
+    @dask_array_available
+    def test_flood_depth_dask_float32_input(self):
+        hand_data = np.array([[0.0, 2.0], [5.0, 10.0]], dtype=np.float32)
+        hand = create_test_raster(hand_data, backend='dask', name='hand')
+        assert flood_depth(hand, water_level=5.0).dtype == np.float64
+
+    @cuda_and_cupy_available
+    def test_flood_depth_cupy_float32_input(self):
+        hand_data = np.array([[0.0, 2.0], [5.0, 10.0]], dtype=np.float32)
+        hand = create_test_raster(hand_data, backend='cupy', name='hand')
+        assert flood_depth(hand, water_level=5.0).dtype == np.float64
+
+    @cuda_and_cupy_available
+    @dask_array_available
+    def test_flood_depth_dask_cupy_float32_input(self):
+        hand_data = np.array([[0.0, 2.0], [5.0, 10.0]], dtype=np.float32)
+        hand = create_test_raster(hand_data, backend='dask+cupy', name='hand')
+        assert flood_depth(hand, water_level=5.0).dtype == np.float64
+
+    def test_curve_number_runoff_numpy_float32_input(self):
+        p_data = np.array([[10.0, 100.0]], dtype=np.float32)
+        rainfall = create_test_raster(p_data, backend='numpy', name='rain')
+        assert curve_number_runoff(rainfall, 80.0).dtype == np.float64
+
+    @dask_array_available
+    def test_curve_number_runoff_dask_float32_input(self):
+        p_data = np.array([[10.0, 100.0]], dtype=np.float32)
+        rainfall = create_test_raster(p_data, backend='dask', name='rain')
+        assert curve_number_runoff(rainfall, 80.0).dtype == np.float64
+
+    @dask_array_available
+    def test_curve_number_runoff_dask_float32_cn_dataarray(self):
+        p_data = np.array([[10.0, 100.0]], dtype=np.float32)
+        cn_data = np.array([[80.0, 90.0]], dtype=np.float32)
+        rainfall = create_test_raster(p_data, backend='dask', name='rain')
+        cn = create_test_raster(cn_data, backend='dask', name='cn')
+        assert curve_number_runoff(rainfall, cn).dtype == np.float64
+
+    @cuda_and_cupy_available
+    def test_curve_number_runoff_cupy_float32_input(self):
+        p_data = np.array([[10.0, 100.0]], dtype=np.float32)
+        rainfall = create_test_raster(p_data, backend='cupy', name='rain')
+        assert curve_number_runoff(rainfall, 80.0).dtype == np.float64
+
+    @cuda_and_cupy_available
+    @dask_array_available
+    def test_curve_number_runoff_dask_cupy_float32_input(self):
+        p_data = np.array([[10.0, 100.0]], dtype=np.float32)
+        rainfall = create_test_raster(p_data, backend='dask+cupy', name='rain')
+        assert curve_number_runoff(rainfall, 80.0).dtype == np.float64