diff --git a/xrspatial/convolution.py b/xrspatial/convolution.py
index 9dd66466..3df13ed6 100644
--- a/xrspatial/convolution.py
+++ b/xrspatial/convolution.py
@@ -17,6 +17,13 @@ class cupy(object):
         ndarray = False
 
 
+def _promote_float(dtype):
+    """Return at least float32; preserve float64."""
+    if np.issubdtype(dtype, np.floating):
+        return dtype
+    return np.float32
+
+
 DEFAULT_UNIT = 'meter'
 METER = 1
 FOOT = 0.3048
@@ -286,8 +293,7 @@ def custom_kernel(kernel):
 @jit(nopython=True, nogil=True)
 def _convolve_2d_numpy(data, kernel):
     # apply kernel to data image.
-    # TODO: handle nan
-    data = data.astype(np.float32)
+    # Caller must ensure data is a float type (float32 or float64).
     nx = data.shape[0]
     ny = data.shape[1]
     nkx = kernel.shape[0]
@@ -295,7 +301,7 @@ def _convolve_2d_numpy(data, kernel):
     wkx = nkx // 2
     wky = nky // 2
 
-    out = np.zeros(data.shape, dtype=np.float32)
+    out = np.empty_like(data)
     out[:] = np.nan
     for i in prange(wkx, nx-wkx):
         iimin = max(i - wkx, 0)
@@ -315,6 +321,7 @@ def _convolve_2d_numpy(data, kernel):
 
 
 def _convolve_2d_numpy_boundary(data, kernel, boundary='nan'):
+    data = data.astype(_promote_float(data.dtype))
     if boundary == 'nan':
         return _convolve_2d_numpy(data, kernel)
     pad_h = kernel.shape[0] // 2
@@ -329,7 +336,7 @@ def _convolve_2d_numpy_boundary(data, kernel, boundary='nan'):
 
 
 def _convolve_2d_dask_numpy(data, kernel, boundary='nan'):
-    data = data.astype(np.float32)
+    data = data.astype(_promote_float(data.dtype))
     pad_h = kernel.shape[0] // 2
     pad_w = kernel.shape[1] // 2
     _func = partial(_convolve_2d_numpy, kernel=kernel)
@@ -393,8 +400,9 @@ def _convolve_2d_cupy(data, kernel, boundary='nan'):
         c1 = -pad_w if pad_w else None
         return result[r0:r1, c0:c1]
 
-    data = data.astype(cupy.float32)
-    out = cupy.empty(data.shape, dtype='f4')
+    fdtype = _promote_float(data.dtype)
+    data = data.astype(fdtype)
+    out = cupy.empty(data.shape, dtype=fdtype)
     out[:, :] = cupy.nan
     griddim, blockdim = cuda_args(data.shape)
     _convolve_2d_cuda[griddim, blockdim](data, kernel, cupy.asarray(out))
@@ -402,7 +410,7 @@ def _convolve_2d_cupy(data, kernel, boundary='nan'):
 
 
 def _convolve_2d_dask_cupy(data, kernel, boundary='nan'):
-    data = data.astype(cupy.float32)
+    data = data.astype(_promote_float(data.dtype))
     pad_h = kernel.shape[0] // 2
     pad_w = kernel.shape[1] // 2
     _func = partial(_convolve_2d_cupy, kernel=kernel)
diff --git a/xrspatial/tests/test_focal.py b/xrspatial/tests/test_focal.py
index 5517713d..47d0e3a4 100644
--- a/xrspatial/tests/test_focal.py
+++ b/xrspatial/tests/test_focal.py
@@ -840,6 +840,62 @@ def test_convolution_2d_boundary_no_nan(boundary):
         np_result.data, da_result.data.compute(), equal_nan=True, rtol=1e-5)
 
 
+# --- convolve_2d float64 preservation (issue-1096) ---
+
+
+@pytest.mark.parametrize("backend", ['numpy', 'dask+numpy', 'cupy', 'dask+cupy'])
+def test_convolve_2d_preserves_float64_1096(backend):
+    """Float64 input should produce float64 output, not float32.
+
+    Regression test for #1096: convolve_2d hardcoded .astype(float32)
+    across all backends.
+    """
+    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")
+    if 'dask' in backend and da is None:
+        pytest.skip("Requires Dask")
+
+    # Values near 1e7 where float32 loses the 0.0x differences
+    data = np.array([[1e7 + 0.01, 1e7 + 0.02, 1e7 + 0.03],
+                     [1e7 + 0.04, 1e7 + 0.05, 1e7 + 0.06],
+                     [1e7 + 0.07, 1e7 + 0.08, 1e7 + 0.09],
+                     [1e7 + 0.10, 1e7 + 0.11, 1e7 + 0.12]], dtype=np.float64)
+    kernel = np.array([[0, 1, 0], [1, 1, 1], [0, 1, 0]], dtype=np.float64)
+
+    agg = create_test_raster(data, backend=backend, chunks=(4, 3))
+    result = convolve_2d(agg.data, kernel)
+
+    if hasattr(result, 'compute'):
+        result = result.compute()
+    if hasattr(result, 'get'):
+        result = result.get()
+    result = np.asarray(result)
+
+    # Output must be float64
+    assert result.dtype == np.float64, f"got {result.dtype}"
+
+    # Interior pixel (1,1): kernel cross hits [0.02, 0.04, 0.05, 0.06, 0.08]
+    # Expected: sum = 5e7 + 0.25
+    expected_center = 5e7 + 0.25
+    assert abs(float(result[1, 1]) - expected_center) < 1e-8, (
+        f"center={result[1, 1]}, expected={expected_center}"
+    )
+
+
+def test_convolve_2d_int_promotes_to_float32_1096():
+    """Integer input should be promoted to float32 (not stay int)."""
+    data = np.array([[1, 2, 3],
+                     [4, 5, 6],
+                     [7, 8, 9]], dtype=np.int32)
+    kernel = np.array([[0, 1, 0], [1, 1, 1], [0, 1, 0]], dtype=np.float64)
+
+    result = convolve_2d(data, kernel)
+    assert np.issubdtype(result.dtype, np.floating), f"got {result.dtype}"
+    # Interior pixel (1,1): cross sum = 2+4+5+6+8 = 25
+    assert float(result[1, 1]) == 25.0
+
+
 # --- 3D (multi-band) focal tests ---