diff --git a/.claude/sweep-error-handling-state.csv b/.claude/sweep-error-handling-state.csv
index 38d2bfd25..63471c748 100644
--- a/.claude/sweep-error-handling-state.csv
+++ b/.claude/sweep-error-handling-state.csv
@@ -1,5 +1,6 @@
 module,last_inspected,issue,severity_max,categories_found,notes
 bump,2026-07-02,,HIGH,1;2;3;4,"bump() agg template unvalidated: plain ndarray -> ArrayTypeFunctionMapping 'Unsupported Array Type'; 3D/1D DataArray -> 'too many values to unpack'. count/spread unvalidated (internal numpy errors / silent). Added _validate_raster(agg,ndim=2) + _validate_scalar for count,spread. all 4 backends verified (CUDA present)."
 convolution,2026-07-02,,HIGH,1;2;3;4,"convolve_2d/convolution_2d skipped kernel + DataArray validation: None/1D/3D/list kernel -> numba TypingError, even kernel silently off-center (custom_kernel rejects it), numpy agg -> memoryview astype error. Fixed via _validate_kernel + _validate_raster; branch deep-sweep-error-handling-convolution-2026-07-02 pushed to fork; issue/PR create blocked by auto-mode, open from parent. MEDIUM(unfixed): annulus_kernel inner>outer -> cryptic np.pad 'index cant contain negative values'. LOW: circle_kernel cellsize=0 ZeroDivisionError, cellsize<0 cryptic linspace; calc_cellsize non-DataArray -> AttributeError attrs. cupy verified."
+edge_detection,2026-07-18,3676,MEDIUM,3,"All 5 funcs validate via _validate_raster; convolve_2d validates boundary/kernel; observed failure messages name func/param/offending value, consistent across siblings. MEDIUM(fixed #3676): NaN input propagates to 3x3 footprint, undocumented in docstrings -> added Notes sections + exact-footprint test. LOW(unfixed): (0,0) raster diverges by backend (numpy empty result; dask ValueError 'overlapping depth 1 is larger than your array 0'; cupy CudaAPIError INVALID_VALUE) - engine-level in convolve_2d, adversarial input. LOW(unfixed): int->float32 promotion documented only in convolve_2d docstring. Battery executed on numpy/dask+numpy/cupy/dask+cupy (CUDA present)."
 geotiff,2026-07-02,3604,MEDIUM,2;4,"to_geotiff 0D/1D DataArray raised opaque IndexError from _coords.py coords_to_transform (dims[-2]) instead of clean 'Expected 2D or 3D' ValueError; numpy path + 4D DataArray already clean. Fixed via early ndim guard before dispatch (eager/vrt/gpu) + 3 tests; PR #3604. Read-side param validation + typed-error hierarchy + allow_rotated/allow_invalid_nodata VRT+chunked opt-in threading verified clean (CUDA available, GPU paths run). gh issue create blocked by auto-mode; PR opened. Cat 2+4."
 pathfinding,2026-07-08,3649,CRITICAL,1;2;3;4,"CRITICAL: string barriers (['0']) silently ignored -> path crosses wall, no error/warning (numba float==unicode always False); scalar barriers=0 -> numba TypingError. HIGH: search_radius unvalidated: -1 silently all-NaN 'no path' (numpy+cupy), other geometries 'negative dimensions are not allowed'; float radius -> slice TypeError; via multi_stop -> misleading 'no path between waypoints'. MEDIUM: multi_stop_search missing dims check -> bare KeyError 'y' vs a_star ValueError; scalar start/goal -> 'float object is not subscriptable' in _get_pixel_id. All fixed: _validate_barriers/_validate_search_radius/_validate_point/_validate_surface_dims + tests. LOW (unfixed): (0,0)-size raster -> 'zero-size array to reduction' from calc_res. All 4 backends executed (CUDA present). Battery: numpy/dask/cupy/dask+cupy."
diff --git a/xrspatial/edge_detection.py b/xrspatial/edge_detection.py
index 8e7e766dd..d6fe0d98b 100644
--- a/xrspatial/edge_detection.py
+++ b/xrspatial/edge_detection.py
@@ -71,6 +71,13 @@ def sobel_x(agg, name='sobel_x', boundary='nan'):
         integers as float32, 32/64-bit integers as float64 so that
         large values keep unit precision (exact up to 2**53 for
         64-bit integers).
+
+    Notes
+    -----
+    NaN cells in the input propagate: every output cell whose 3x3
+    neighborhood (as extended by the boundary mode) contains a NaN
+    becomes NaN. With the default ``boundary='nan'``, the outer
+    one-cell ring of the output is also NaN.
     """
     _validate_raster(agg, func_name='sobel_x', name='agg')
     out = convolve_2d(_promote_wide_int(agg.data), SOBEL_X, boundary)
@@ -104,6 +111,13 @@ def sobel_y(agg, name='sobel_y', boundary='nan'):
         integers as float32, 32/64-bit integers as float64 so that
         large values keep unit precision (exact up to 2**53 for
         64-bit integers).
+
+    Notes
+    -----
+    NaN cells in the input propagate: every output cell whose 3x3
+    neighborhood (as extended by the boundary mode) contains a NaN
+    becomes NaN. With the default ``boundary='nan'``, the outer
+    one-cell ring of the output is also NaN.
     """
     _validate_raster(agg, func_name='sobel_y', name='agg')
     out = convolve_2d(_promote_wide_int(agg.data), SOBEL_Y, boundary)
@@ -137,6 +151,13 @@ def laplacian(agg, name='laplacian', boundary='nan'):
         integers as float32, 32/64-bit integers as float64 so that
         large values keep unit precision (exact up to 2**53 for
         64-bit integers).
+
+    Notes
+    -----
+    NaN cells in the input propagate: every output cell whose 3x3
+    neighborhood (as extended by the boundary mode) contains a NaN
+    becomes NaN. With the default ``boundary='nan'``, the outer
+    one-cell ring of the output is also NaN.
     """
     _validate_raster(agg, func_name='laplacian', name='agg')
     out = convolve_2d(_promote_wide_int(agg.data), LAPLACIAN_KERNEL, boundary)
@@ -170,6 +191,13 @@ def prewitt_x(agg, name='prewitt_x', boundary='nan'):
         integers as float32, 32/64-bit integers as float64 so that
         large values keep unit precision (exact up to 2**53 for
         64-bit integers).
+
+    Notes
+    -----
+    NaN cells in the input propagate: every output cell whose 3x3
+    neighborhood (as extended by the boundary mode) contains a NaN
+    becomes NaN. With the default ``boundary='nan'``, the outer
+    one-cell ring of the output is also NaN.
     """
     _validate_raster(agg, func_name='prewitt_x', name='agg')
     out = convolve_2d(_promote_wide_int(agg.data), PREWITT_X, boundary)
@@ -203,6 +231,13 @@ def prewitt_y(agg, name='prewitt_y', boundary='nan'):
         integers as float32, 32/64-bit integers as float64 so that
         large values keep unit precision (exact up to 2**53 for
         64-bit integers).
+
+    Notes
+    -----
+    NaN cells in the input propagate: every output cell whose 3x3
+    neighborhood (as extended by the boundary mode) contains a NaN
+    becomes NaN. With the default ``boundary='nan'``, the outer
+    one-cell ring of the output is also NaN.
     """
     _validate_raster(agg, func_name='prewitt_y', name='agg')
     out = convolve_2d(_promote_wide_int(agg.data), PREWITT_Y, boundary)
diff --git a/xrspatial/tests/test_edge_detection.py b/xrspatial/tests/test_edge_detection.py
index e7701e47c..6c459b6a8 100644
--- a/xrspatial/tests/test_edge_detection.py
+++ b/xrspatial/tests/test_edge_detection.py
@@ -162,8 +162,11 @@ def test_nan_in_input_propagates(self, func):
         data[2, 3] = np.nan
         agg = create_test_raster(data, backend='numpy')
         result = func(agg, boundary='reflect')
-        # NaN should propagate to neighbors
-        assert np.isnan(result.data[2, 3])
+        # per the docstring: every output cell whose 3x3 neighborhood
+        # contains the NaN becomes NaN, and no other cell does
+        expected_nan = np.zeros((5, 6), dtype=bool)
+        expected_nan[1:4, 2:5] = True
+        np.testing.assert_array_equal(np.isnan(result.data), expected_nan)
 
     @pytest.mark.parametrize('func', [sobel_x, sobel_y, laplacian, prewitt_x, prewitt_y])
     def test_all_nan_input(self, func):