Fix flood dask backends emitting float32 output for float32 input

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)

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