diff --git a/xrspatial/slope.py b/xrspatial/slope.py
index 56b41df3f..90a6aac13 100644
--- a/xrspatial/slope.py
+++ b/xrspatial/slope.py
@@ -2,7 +2,7 @@
 
 # std lib
 from functools import partial
-from math import atan
+from math import atan, isnan, nan
 from typing import Optional, Union
 
 # 3rd-party
@@ -52,6 +52,8 @@ def _cpu(data, cellsize_x, cellsize_y):
     rows, cols = data.shape
     for y in range(1, rows - 1):
         for x in range(1, cols - 1):
+            if np.isnan(data[y, x]):
+                continue
             a = data[y + 1, x - 1]
             b = data[y + 1, x]
             c = data[y + 1, x + 1]
@@ -112,6 +114,8 @@ def _run_dask_cupy(data: da.Array,
 
 @cuda.jit(device=True)
 def _gpu(arr, cellsize_x, cellsize_y):
+    if isnan(arr[1, 1]):
+        return nan
     a = arr[2, 0]
     b = arr[2, 1]
     c = arr[2, 2]
diff --git a/xrspatial/tests/test_slope.py b/xrspatial/tests/test_slope.py
index ed90af00d..3cf1a1d32 100644
--- a/xrspatial/tests/test_slope.py
+++ b/xrspatial/tests/test_slope.py
@@ -274,3 +274,48 @@ def test_degenerate_shape_geodesic(shape):
     general_output_checks(raster, result)
     assert result.shape == shape
     assert np.all(np.isnan(result.data))
+
+
+# A NoData (NaN) center cell must stay NaN in the slope output, even when all 8
+# of its neighbours are valid. The planar kernels read the center cell, so a
+# hole in the DEM can't masquerade as valid flat terrain. Regression for #2761.
+def _center_nan_data():
+    data = np.zeros((5, 5), dtype=np.float64)
+    data[2, 2] = np.nan
+    return data
+
+
+def test_center_nan_propagates_numpy():
+    agg = create_test_raster(_center_nan_data(), backend='numpy', attrs={'res': (1, 1)})
+    result = slope(agg)
+    general_output_checks(agg, result)
+    assert np.isnan(result.data[2, 2])
+
+
+@dask_array_available
+def test_center_nan_propagates_dask_numpy():
+    data = _center_nan_data()
+    numpy_agg = create_test_raster(data, backend='numpy', attrs={'res': (1, 1)})
+    dask_agg = create_test_raster(data, backend='dask+numpy',
+                                  attrs={'res': (1, 1)}, chunks=(3, 3))
+    assert_numpy_equals_dask_numpy(numpy_agg, dask_agg, slope, nan_edges=False)
+    assert np.isnan(slope(dask_agg).data.compute()[2, 2])
+
+
+@cuda_and_cupy_available
+def test_center_nan_propagates_cupy():
+    data = _center_nan_data()
+    numpy_agg = create_test_raster(data, backend='numpy', attrs={'res': (1, 1)})
+    cupy_agg = create_test_raster(data, backend='cupy', attrs={'res': (1, 1)})
+    assert_numpy_equals_cupy(numpy_agg, cupy_agg, slope, nan_edges=False)
+    assert np.isnan(slope(cupy_agg).data.get()[2, 2])
+
+
+@dask_array_available
+@cuda_and_cupy_available
+def test_center_nan_propagates_dask_cupy():
+    data = _center_nan_data()
+    numpy_agg = create_test_raster(data, backend='numpy', attrs={'res': (1, 1)})
+    dask_cupy_agg = create_test_raster(data, backend='dask+cupy',
+                                       attrs={'res': (1, 1)}, chunks=(3, 3))
+    assert_numpy_equals_dask_cupy(numpy_agg, dask_cupy_agg, slope, nan_edges=False)