Parallelize convolve_2d numpy kernel (#3615) (#3616)

* Parallelize convolve_2d numpy kernel (#3615)

_convolve_2d_numpy iterated with numba.prange but was decorated
@jit(nopython=True, nogil=True) with no parallel=True, so prange
degraded to a serial range. On a 20-core host a 2000x2000 float64
raster with a 15x15 kernel dropped from ~356 ms to ~53 ms end-to-end
once parallel=True was enabled, with identical results.

Add parallel=True and serialize the kernel launch behind a
module-level threading.Lock, since dask calls it per chunk under a
threaded scheduler and numba's workqueue layer is not threadsafe
across host threads (SIGABRT on macOS, same hazard as #3141 fixed in
terrain.py). A single numpy call takes the lock uncontended and still
runs across all cores.

Add benchmarks/benchmarks/convolution.py (Convolve2d, CircleKernel)
so the regression cannot silently return.

* bench(convolution): force dask compute in convolve_2d benchmark

The dask case of time_convolve_2d returned a lazy array, so asv timed
graph construction instead of the kernel this PR parallelizes. Force the
compute for backends that produce a lazy array, matching the convention
in the twi/interpolate/flood benchmarks. numpy and cupy results have no
.compute and are unaffected. (#3615)

Author: brendancol

.claude/sweep-performance-state.csv

@@ -5,7 +5,7 @@ bilateral,2026-03-31T18:00:00Z,SAFE,compute-bound,0,,
 bump,2026-07-02,SAFE,graph-bound,0,3612,"Perf re-audit 2026-07-02: MEDIUM (#3612/PR) _partition_bumps rescanned all locs per chunk -> O(n_chunks*count) dask graph build (7.8s@25.6k chunks, count=100k); fixed with searchsorted bucketing, byte-identical partitions, 85x@25.6k chunks. Memory SAFE (per-chunk lazy build + guard). LOW: cupy path runs CPU kernel then single host->device transfer (by design, no GPU kernel); np.random.choice(range) slower than randint; unused rows,cols in _finish_bump. cuda-available, cupy+dask+cupy parity verified."
 classify,2026-06-20,RISKY,graph-bound,1,3412,"Re-audit 2026-06-20 (CUDA host). 1 HIGH: _generate_sample_indices >10M branch used RandomState.choice(replace=False) which builds a full arange(num_data) permutation -> O(num_data) host alloc (160MB for 20M pop, OOM at 30TB) despite docstring claiming O(num_sample). Backed dask/dask+cupy natural_breaks/maximum_breaks/quantile/percentiles/box_plot. Fixed via np.random.default_rng().choice (Floyd, O(num_sample), still deterministic); peak 160MB->0.4MB. Other paths SAFE: head_tail_breaks already persists+fuses; box_plot samples; cupy kernels low-register; no .values/np.asarray-on-dask/.compute-in-loop. 93 classify tests pass incl GPU."
 contour,2026-03-31T18:00:00Z,SAFE,compute-bound,0,,
-convolution,2026-03-31T18:00:00Z,SAFE,compute-bound,0,,
+convolution,2026-07-02,SAFE,compute-bound,1,3615,_convolve_2d_numpy used prange w/o parallel=True -> ran serial (~7-10x slow); fixed via parallel=True + threading.Lock (macOS SIGABRT hazard #3141); cuda kernel 40 regs OK; dask ~20 tasks/chunk
 corridor,2026-03-31T18:00:00Z,SAFE,compute-bound,0,,
 cost_distance,2026-06-15,RISKY,memory-bound,1,3342,"Perf sweep 2026-06-15. HIGH: bounded map_overlap branch in _cost_distance_dask gated on full dims (pad>=height/width) not chunk size; pad>chunk collapses to single chunk (#880-class OOM, verified npartitions=1 at chunks=10/pad=96). Fixed: compare pad vs max chunk dim, route to iterative when pad>=chunk (matches GPU path L484). dask+cupy path already correct. Register count 37 (no pressure). nanmin().compute() L478/L1149 intentional scalar. iterative tile_cache full-dataset materialization is documented MemoryError-guarded design (#1118). All 56 tests pass incl GPU."
 curvature,2026-03-31T18:00:00Z,SAFE,compute-bound,0,,

benchmarks/benchmarks/convolution.py

@@ -0,0 +1,32 @@
+import numpy as np
+
+from xrspatial.convolution import circle_kernel, convolve_2d
+
+from .common import get_xr_dataarray
+
+
+class Convolve2d:
+    params = ([300, 1000, 3000], [(5, 5), (25, 25)], ["numpy", "cupy", "dask"])
+    param_names = ("nx", "kernelsize", "type")
+
+    def setup(self, nx, kernelsize, type):
+        ny = nx // 2
+        self.agg = get_xr_dataarray((ny, nx), type)
+        kernel_h, kernel_w = kernelsize
+        self.kernel = np.ones((kernel_h, kernel_w), dtype=np.float64)
+
+    def time_convolve_2d(self, nx, kernelsize, type):
+        # convolve_2d takes the backing array, not the DataArray wrapper.
+        result = convolve_2d(self.agg.data, self.kernel)
+        # dask returns a lazy array; force the compute so the benchmark
+        # times the kernel, not just graph construction.
+        if hasattr(result, 'compute'):
+            result.compute()
+
+
+class CircleKernel:
+    params = ([3, 25, 100],)
+    param_names = ("radius",)
+
+    def time_circle_kernel(self, radius):
+        circle_kernel(1, 1, radius)

xrspatial/convolution.py

@@ -1,4 +1,5 @@
 import re
+import threading
 from functools import partial
 
 import numpy as np
@@ -343,7 +344,16 @@ def custom_kernel(kernel):
     return kernel
 
 
-@jit(nopython=True, nogil=True)
+# Numba parallel=True kernels must not be launched concurrently from multiple
+# Python threads: the default 'workqueue' threading layer is not threadsafe and
+# aborts the process (SIGABRT on macOS) when two host threads enter a parallel
+# region at once.  _convolve_2d_dask_numpy calls the kernel per chunk under
+# dask's threaded scheduler, so the kernel launch is serialized behind this
+# lock.  Same hazard and fix as the terrain and reproject kernels (#3141).
+_PARALLEL_KERNEL_LOCK = threading.Lock()
+
+
+@jit(nopython=True, nogil=True, parallel=True)
 def _convolve_2d_numpy(data, kernel):
     # apply kernel to data image.
     # Caller must ensure data is a float type (float32 or float64).
@@ -373,14 +383,23 @@ def _convolve_2d_numpy(data, kernel):
     return out
 
 
+def _convolve_2d_numpy_locked(data, kernel):
+    # Serialize the parallel=True kernel launch across host threads; see the
+    # comment on _PARALLEL_KERNEL_LOCK. A single numpy call takes the lock
+    # uncontended and still runs across all cores; concurrent dask chunk calls
+    # run one at a time, each internally parallel.
+    with _PARALLEL_KERNEL_LOCK:
+        return _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)
+        return _convolve_2d_numpy_locked(data, kernel)
     pad_h = kernel.shape[0] // 2
     pad_w = kernel.shape[1] // 2
     padded = _pad_array(data, (pad_h, pad_w), boundary)
-    result = _convolve_2d_numpy(padded, kernel)
+    result = _convolve_2d_numpy_locked(padded, kernel)
     r0 = pad_h if pad_h else None
     r1 = -pad_h if pad_h else None
     c0 = pad_w if pad_w else None
@@ -392,7 +411,7 @@ def _convolve_2d_dask_numpy(data, kernel, boundary='nan'):
     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)
+    _func = partial(_convolve_2d_numpy_locked, kernel=kernel)
     out = data.map_overlap(_func,
                            depth=(pad_h, pad_w),
                            boundary=_boundary_to_dask(boundary),