diff --git a/.claude/sweep-security-state.csv b/.claude/sweep-security-state.csv
index 1ccbec891..50256deaf 100644
--- a/.claude/sweep-security-state.csv
+++ b/.claude/sweep-security-state.csv
@@ -11,7 +11,7 @@ cost_distance,2026-04-25,1262,HIGH,1,,"HIGH (fixed #1262): the recent #1252/#125
 curvature,2026-04-25,,,,,"Clean. Small (271 LOC) module computing 3x3 second-derivative stencil. Cat 1: only single output buffer matching input shape (np.empty at line 37, cupy.empty at line 101) -- bounded by caller, per audit guidance not a finding. Cat 2: _cpu numba kernel uses range(1, rows-1)/range(1, cols-1) with simple (y, x) indices; no flat indexing or queue arrays; numba range loops produce int64. Cat 3: division by cellsize*cellsize on line 44 -- cellsize comes from get_dataarray_resolution() (raster property, not user-direct); cellsize=0 is unrealistic and would produce inf consistently across backends. NaN inputs propagate correctly through float arithmetic. Cat 4: _run_gpu (line 79-86) has full bounds guard via 'i + di <= out.shape[0] - 1 and j + dj <= out.shape[1] - 1' which guarantees i < shape[0] and j < shape[1] before the out[i, j] write; no shared memory; out is pre-filled with NaN at line 102 so threads outside the guard correctly leave NaN. Cat 5: no file I/O. Cat 6: curvature() calls _validate_raster at line 253; all four backend paths explicitly cast to float32 (lines 51, 62, 97, 112) so dtype is normalized before any computation; tests cover int32/int64/uint32/uint64/float32/float64 across numpy/cupy/dask+numpy/dask+cupy."
 dasymetric,2026-04-25,1261,HIGH,1;6,,"HIGH (fixed #1261): pycnophylactic() and disaggregate(method='limiting_variable') allocated full-shape working arrays without checking available memory first. _pycnophylactic_numpy additionally stored one full-shape bool mask per zone in zone_masks, so peak memory grew with N_zones * H * W (1000 zones on a 10000x10000 raster ~ 100 GB just for masks on top of ~3.4 GB of iteration buffers). Fixed by adding _available_memory_bytes() helper and two budget functions (_check_disaggregate_memory, _check_pycnophylactic_memory) that raise MemoryError before the first allocation when projected working memory exceeds 50% of available RAM. The disaggregate guard runs only for in-RAM backends (numpy, cupy); dask paths process per-chunk and are skipped. The pycnophylactic guard scales with len(values_dict) so an exploding zone count is rejected even on a small raster. MEDIUM (unfixed, Cat 6): disaggregate() and pycnophylactic() do not call _validate_raster on zones/weight; they only check isinstance(xr.DataArray), ndim, and shape. Object-dtype or other non-numeric input would fail with confusing TypeError from inside numpy.asarray rather than a clean ValueError. Deferred to a separate PR per the security-sweep one-fix-per-PR policy."
 diffusion,2026-04-27,1267,HIGH,1;3,1281,"HIGH (fixed #1267): diffuse() had no memory guard on its core allocations and steps was unbounded. (1) The public API allocated np.full(agg.shape) for scalar diffusivity even when the dispatched backend was dask, forcing a full numpy alpha raster up front -- a 100kx100k input would OOM on an 80 GB allocation before any backend dispatch. (2) _diffuse_step_numpy and _diffuse_cupy allocated per-step buffers with no memory check. (3) steps was validated only with min_val=1, so steps=10**12 was accepted and would loop forever. Fixed by adding _check_memory/_check_gpu_memory helpers (cost_distance pattern, ~32 B/pixel budget for u + out + alpha + padded copy at 50% of available RAM/VRAM), deferring the np.full alpha allocation until after the guard runs in eager paths, teaching _diffuse_dask_cupy to handle scalar alpha lazily via cp.full per chunk (mirroring _diffuse_dask_numpy), and capping steps at _MAX_STEPS = 100_000 in _validate_scalar. GPU kernel _diffuse_step_gpu has bounds guard (if i < rows and j < cols), no shared memory, _validate_raster called on agg and on diffusivity DataArray, NaN check uses val != val correctly, no file I/O, no int32 indexing. Follow-up HIGH (fixed #1281): user-supplied dt was validated only as > 0, but explicit forward-Euler is unconditionally unstable above 0.25 * dx**2 / max(alpha); the dt=None branch already used this exact bound, so the fix hoists it into cfl_max and raises ValueError when the user-supplied dt exceeds it. Single check in the public entrypoint covers all four backends."
-edge_detection,2026-04-25,1271,MEDIUM,6,,"MEDIUM (fixed #1271): the five public functions sobel_x, sobel_y, laplacian, prewitt_x, prewitt_y did not call _validate_raster on agg. Non-DataArray inputs raised AttributeError from agg.data and wrong-ndim DataArrays failed inside numba/cupy with confusing errors instead of clean TypeError/ValueError. Numerical correctness was unaffected because convolve_2d._promote_float casts integer dtypes to float32 before the kernel runs. Fixed by adding _validate_raster(agg, func_name=..., name='agg') at the top of each function. No CRITICAL/HIGH findings: convolve_2d enforces 3x3 odd kernels and 2D agg.data, allocations match input shape, no CUDA kernels owned by this module, no file I/O."
+edge_detection,2026-07-18,3680,HIGH,6,3684,"HIGH (fixed #3680, PR #3686): all five wrappers pass agg.data to convolve_2d, whose _promote_float casts every int dtype to float32; integers above 2**24 lose low bits so unit-step gradients on large int32/int64 rasters returned exactly 0.0 on all four backends (repro: sobel_x on int32 ramp at offset 1e8 -> 0.0 vs float64 ref 8.0; CUDA available, cupy and dask+cupy executed for real). Fixed by pre-promoting 32/64-bit ints to float64 in edge_detection; 8/16-bit ints stay float32. Root cause lives in convolution._promote_float and also hits convolution_2d and focal.mean -> follow-up #3684 (contract ratified by #1096/#2769/#3214, needs maintainer call). Other cats clean: no allocations, index math, file IO, or CUDA kernels owned by the module; convolve_2d validates raster, boundary, and kernel; NaN propagates under all boundary modes; adversarial shapes (1x1, 1xN, Nx1, 2x2, primes) OK; invalid boundary strings and bool/datetime/complex dtypes rejected; 119 tests pass incl new wide-int coverage."
 emerging_hotspots,2026-04-25,1274,HIGH,1,,"HIGH (fixed #1274): emerging_hotspots() public API only validated ndim and shape[0] >= 2. The numpy and cupy backends each materialised three full (T, H, W) cubes (a float32 input copy, gi_zscore float32, gi_bin int8) plus H*W temporaries with no memory check; a (100, 20000, 20000) input projected to ~480 GB. Fixed by adding _available_memory_bytes()/_check_memory(n_times, ny, nx) (12 bytes per cube cell budget) and calling it from the public API for non-dask inputs. Dask paths skip the guard because their map_blocks/map_overlap chunk functions do not materialise the full cube. MEDIUM (unfixed, Cat 6): public API does not call _validate_raster() so non-numeric dtypes fail later with a confusing error rather than a clean TypeError. No GPU kernels in this module (uses convolve_2d). No file I/O. Cat 3 statistical paths are robust: _mann_kendall_statistic_numpy guards var_s <= 0 before sqrt, both numpy and cupy backends raise ZeroDivisionError on global_std == 0, and _mk_pvalue handles z==0 explicitly."
 erosion,2026-04-25,1275,HIGH,1;3;6,,"HIGH (fixed #1275): erode() accepted three user-controlled parameters with no upper bound. (1) iterations sized rng.random((iterations, 2)) on the host (16 B/particle) and was copied to the GPU via cupy.asarray, so iterations=10**12 attempted ~16 TB on each side. (2) params['radius'] drove _build_brush which iterates (2r+1)**2 cells and stores three arrays of the same length, so radius=10**6 allocated ~12 TB of brush data. (3) params['max_lifetime'] is the inner per-particle JIT loop in both _erode_cpu and _erode_gpu_kernel, so max_lifetime=10**12 with the default iterations=50000 ran 5e16 step iterations. The existing _check_erosion_memory helper only fired on dask paths and ignored the random_pos and brush working sets. Fixed by capping all three parameters at the public erode() entry via _validate_scalar(max_val=...) (_MAX_ITERATIONS=1e8, _MAX_RADIUS=1024, _MAX_LIFETIME=1e5), rewriting _check_erosion_memory to include the random_pos buffer and brush bytes in its budget, and wiring the guard into _erode_numpy and _erode_cupy so every backend benefits (the dask paths inherit it via their _erode_numpy/_erode_cupy calls). Mirrors diffuse #1268 pattern. Deferred follow-ups (separate PRs): Cat 3 HIGH NaN input is not guarded in _erode_cpu / _erode_gpu_kernel -- a NaN cell propagates through bilinear interpolation into dir_x/dir_y, NaN bounds checks fall through, and particles can deposit NaN into arbitrary cells via cuda.atomic.add. Cat 6 MEDIUM erode() does not call _validate_raster() on agg -- non-numeric or wrong-ndim input fails inside numba/cupy with a confusing error. No Cat 2 (no int32 flat-index math), no Cat 4 (GPU kernel has bounds guard at line 184 plus per-step bounds checks before every read/write, brush writes are explicitly bounds-checked, no shared memory), no Cat 5 (no file I/O)."
 fire,2026-04-25,,,,,"Clean. Despite the module's size hint, fire.py is purely per-cell raster ops -- not cellular-automaton or front-tracking. Seven public APIs: dnbr, rdnbr, burn_severity_class, fireline_intensity, flame_length, rate_of_spread, kbdi. No iteration, no queues, no multi-channel state, no random numbers, no file paths. Cat 1: every output allocation matches input shape (single buffer, bounded by caller). Anderson-13 fuel table is a fixed 13x8 constant. _rothermel_fuel_constants returns 12 scalars before dispatch (no per-pixel state). Cat 2: no flat-index math, all indexing is 2-D (y, x); no height*width multiplication. Cat 3: rdnbr guards denom < 1e-10; burn_severity_class is threshold-only; flame_length guards v <= 0.0 before fractional power; rate_of_spread guards M_x>0/beta>0/denom>0 and clamps eta_M, U_mmin, R; kbdi clamps Q to [0, 800] and net_P to >= 0. Adversarial wind=inf or T=inf would push exp/power to inf in rate_of_spread/kbdi but inputs are user-controlled rasters, fire model is research-quality (LOW only). Cat 4: all 7 CUDA kernels (_dnbr_gpu L157, _rdnbr_gpu L246, _bsc_gpu L362, _fli_gpu L455, _fl_gpu L552, _ros_gpu L681, _kbdi_gpu L870) have 'y < out.shape[0] and x < out.shape[1]' bounds guard; every kernel is point-wise (no neighbour stencil) so the simple guard is sufficient; no shared memory, no syncthreads needed. Cat 5: no file I/O. Cat 6: every public function calls _validate_raster on each input raster (dnbr/rdnbr/fireline_intensity/rate_of_spread/kbdi pass 2-3 rasters each, all validated), validate_arrays enforces equal shape, _validate_scalar gates heat_content/fuel_model (1-13)/annual_precip, and every input is .astype('f4') before reaching any kernel so dtype is normalized."
diff --git a/xrspatial/edge_detection.py b/xrspatial/edge_detection.py
index e908800f7..8e7e766dd 100644
--- a/xrspatial/edge_detection.py
+++ b/xrspatial/edge_detection.py
@@ -28,6 +28,23 @@
                              [0,  1, 0]], dtype=np.float64)
 
 
+def _promote_wide_int(data):
+    """Cast 32/64-bit integer arrays to float64 before convolution.
+
+    ``convolve_2d`` promotes integer inputs to float32, whose 24-bit
+    mantissa cannot represent integers above 2**24: unit steps between
+    large values vanish in the cast and gradients silently collapse to
+    zero (#3680). float64 is exact for every int32/uint32 value and for
+    int64/uint64 up to 2**53. 8- and 16-bit integers are exactly
+    representable in float32, so they keep ``convolve_2d``'s promotion.
+    Works on numpy, cupy, and dask arrays alike (``astype`` is lazy on
+    dask).
+    """
+    if data.dtype.kind in 'iu' and data.dtype.itemsize > 2:
+        return data.astype(np.float64)
+    return data
+
+
 def sobel_x(agg, name='sobel_x', boundary='nan'):
     """Compute the horizontal gradient of a raster using the Sobel operator.
 
@@ -50,9 +67,13 @@ def sobel_x(agg, name='sobel_x', boundary='nan'):
     -------
     xarray.DataArray
         Horizontal gradient with the same shape and backend as the input.
+        Integer inputs are computed in floating point: 8/16-bit
+        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).
     """
     _validate_raster(agg, func_name='sobel_x', name='agg')
-    out = convolve_2d(agg.data, SOBEL_X, boundary)
+    out = convolve_2d(_promote_wide_int(agg.data), SOBEL_X, boundary)
     return xr.DataArray(out, name=name, coords=agg.coords,
                         dims=agg.dims, attrs=agg.attrs)
 
@@ -79,9 +100,13 @@ def sobel_y(agg, name='sobel_y', boundary='nan'):
     -------
     xarray.DataArray
         Vertical gradient with the same shape and backend as the input.
+        Integer inputs are computed in floating point: 8/16-bit
+        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).
     """
     _validate_raster(agg, func_name='sobel_y', name='agg')
-    out = convolve_2d(agg.data, SOBEL_Y, boundary)
+    out = convolve_2d(_promote_wide_int(agg.data), SOBEL_Y, boundary)
     return xr.DataArray(out, name=name, coords=agg.coords,
                         dims=agg.dims, attrs=agg.attrs)
 
@@ -108,9 +133,13 @@ def laplacian(agg, name='laplacian', boundary='nan'):
     -------
     xarray.DataArray
         Laplacian response with the same shape and backend as the input.
+        Integer inputs are computed in floating point: 8/16-bit
+        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).
     """
     _validate_raster(agg, func_name='laplacian', name='agg')
-    out = convolve_2d(agg.data, LAPLACIAN_KERNEL, boundary)
+    out = convolve_2d(_promote_wide_int(agg.data), LAPLACIAN_KERNEL, boundary)
     return xr.DataArray(out, name=name, coords=agg.coords,
                         dims=agg.dims, attrs=agg.attrs)
 
@@ -137,9 +166,13 @@ def prewitt_x(agg, name='prewitt_x', boundary='nan'):
     -------
     xarray.DataArray
         Horizontal gradient with the same shape and backend as the input.
+        Integer inputs are computed in floating point: 8/16-bit
+        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).
     """
     _validate_raster(agg, func_name='prewitt_x', name='agg')
-    out = convolve_2d(agg.data, PREWITT_X, boundary)
+    out = convolve_2d(_promote_wide_int(agg.data), PREWITT_X, boundary)
     return xr.DataArray(out, name=name, coords=agg.coords,
                         dims=agg.dims, attrs=agg.attrs)
 
@@ -166,8 +199,12 @@ def prewitt_y(agg, name='prewitt_y', boundary='nan'):
     -------
     xarray.DataArray
         Vertical gradient with the same shape and backend as the input.
+        Integer inputs are computed in floating point: 8/16-bit
+        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).
     """
     _validate_raster(agg, func_name='prewitt_y', name='agg')
-    out = convolve_2d(agg.data, PREWITT_Y, boundary)
+    out = convolve_2d(_promote_wide_int(agg.data), PREWITT_Y, boundary)
     return xr.DataArray(out, name=name, coords=agg.coords,
                         dims=agg.dims, attrs=agg.attrs)
diff --git a/xrspatial/tests/test_edge_detection.py b/xrspatial/tests/test_edge_detection.py
index a73e1285f..e7701e47c 100644
--- a/xrspatial/tests/test_edge_detection.py
+++ b/xrspatial/tests/test_edge_detection.py
@@ -271,3 +271,73 @@ def test_numpy_equals_dask_cupy(self, ramp_data, func):
         np_agg = create_test_raster(ramp_data, backend='numpy')
         dcu_agg = create_test_raster(ramp_data, backend='dask+cupy', chunks=(5, 6))
         assert_numpy_equals_dask_cupy(np_agg, dcu_agg, func)
+
+
+# ---------------------------------------------------------------------------
+# Wide-integer precision (issue-3680)
+# ---------------------------------------------------------------------------
+# convolve_2d promotes integer inputs to float32, whose 24-bit mantissa
+# cannot separate integers above 2**24, so unit-step gradients on large
+# int32/int64 values silently collapsed to zero. edge_detection now
+# pre-promotes 32/64-bit integers to float64.
+
+
+def _to_numpy(data):
+    if hasattr(data, 'compute'):
+        data = data.compute()
+    if hasattr(data, 'get'):
+        data = data.get()
+    return data
+
+
+class TestWideIntegerPrecision:
+    @pytest.fixture
+    def wide_int32_data(self):
+        # unit-step ramp on an offset beyond float32's 24-bit mantissa
+        return (np.arange(30).reshape(5, 6) % 6 + 100_000_000).astype(np.int32)
+
+    @pytest.mark.parametrize('func', [sobel_x, sobel_y, laplacian, prewitt_x, prewitt_y])
+    @pytest.mark.parametrize('backend', ['numpy', 'dask+numpy', 'cupy', 'dask+cupy'])
+    def test_large_int32_matches_float64_reference(self, wide_int32_data, func, backend):
+        from xrspatial.tests.general_checks import has_cuda_and_cupy, has_dask_array
+        if 'cupy' in backend and not has_cuda_and_cupy():
+            pytest.skip("Requires CUDA and CuPy")
+        if 'dask' in backend and not has_dask_array():
+            pytest.skip("Requires Dask")
+        ref = func(create_test_raster(
+            wide_int32_data.astype(np.float64), backend='numpy'))
+        agg = create_test_raster(wide_int32_data, backend=backend, chunks=(5, 6))
+        result = _to_numpy(func(agg).data)
+        np.testing.assert_allclose(result, ref.data, equal_nan=True)
+
+    def test_large_int32_gradient_nonzero(self, wide_int32_data):
+        # regression guard for the exact silent failure from issue-3680:
+        # the old float32 path returned 0.0 at every interior cell
+        result = sobel_x(create_test_raster(wide_int32_data, backend='numpy'))
+        interior = result.data[1:-1, 1:-1]
+        assert np.all(interior != 0)
+
+    def test_large_int64_matches_float64_reference(self):
+        data = (np.arange(30).reshape(5, 6) % 6
+                + 10_000_000_000).astype(np.int64)
+        ref = sobel_x(create_test_raster(data.astype(np.float64),
+                                         backend='numpy'))
+        result = sobel_x(create_test_raster(data, backend='numpy'))
+        np.testing.assert_allclose(result.data, ref.data, equal_nan=True)
+
+    @pytest.mark.parametrize('dtype,expected', [
+        (np.int8, np.float32),
+        (np.uint8, np.float32),
+        (np.int16, np.float32),
+        (np.uint16, np.float32),
+        (np.int32, np.float64),
+        (np.uint32, np.float64),
+        (np.int64, np.float64),
+        (np.uint64, np.float64),
+        (np.float32, np.float32),
+        (np.float64, np.float64),
+    ])
+    def test_output_dtype_by_input_dtype(self, dtype, expected):
+        data = np.ones((5, 6), dtype=dtype)
+        result = laplacian(create_test_raster(data, backend='numpy'))
+        assert result.dtype == expected