Fix cumulative_viewshed TypeError on cupy rasters (#3194)

.claude/sweep-accuracy-state.csv

@@ -35,6 +35,6 @@ sky_view_factor,2026-05-01,1407,HIGH,4,Horizon angle ignored cell size; fixed by
 terrain,2026-04-10T12:00:00Z,,,,Perlin/Worley/ridged noise correct. Dask chunk boundaries produce bit-identical results. No precision issues.
 terrain_metrics,2026-04-30,,LOW,2;5,"LOW: Inf input not rejected, propagates as Inf (consistent across backends but undocumented). LOW: dask+cupy non-nan boundary path double-pads (wasted compute, central output values still correct). No CRIT/HIGH; tests cover NaN propagation, all 4 backends, all 4 boundary modes, dtype acceptance."
 viewshed,2026-05-29,2691,HIGH,3;5,max_distance window sized from coarser axis clipped cells on anisotropic rasters (PR #2702). LOW unfixed: distance_sweep ring radius same max(res) pattern but max_distance arg always None; _calculate_event_row_col line 880 abs(x>1) precedence bug is a broken guard only. cuda+rtx paths validated.
-visibility,2026-04-13T12:00:00Z,,,,"Bresenham line, LOS kernel, Fresnel zone all correct. All backends converge to numpy."
+visibility,2026-06-10,3194,HIGH,5,"Cat5 backend bug (HIGH): cumulative_viewshed + visibility_frequency raised TypeError on cupy input. count=np.zeros accumulator added to cupy-backed viewshed() result; cupy refuses numpy+cupy mixed add. cupy path had no test coverage. Fix #3194: pull vs_data to numpy via .get() before accumulate (matches _extract_transect). Verified on GPU (RTX A6000): both funcs run, int32/float64 dtypes, valid ranges; 26 visibility tests pass. Cats 1-4 clean: Bresenham endpoints/coverage correct; LOS horizon-angle test correct (>= inclusive, observer cell forced visible); NaN terrain marks that cell invisible without poisoning max_angle downstream; distance accumulation fine. LOW (documented, not fixed): visible[] uses terrain elevations not tgt_h so target_elev does not affect the LOS mask (only los_height/fresnel); fresnel_clear=True over a NaN terrain cell (NaN<x is False); distance/Fresnel assume meter grid like rest of viewshed (geographic-degree coords give physically wrong Fresnel radii) - consistent with module-wide contract."
 worley,2026-05-01,,MEDIUM,2;5,"MEDIUM: numpy backend uses np.empty_like(data) so integer input dtype produces integer output (distances truncated to 0); cupy/dask paths always produce float32. LOW: freq=inf produces 100000 sentinel (sqrt of initial min_dist=1e10), no validation of freq/seed for non-finite values."
 zonal,2026-05-27,2528,MEDIUM,5,"Pass 2 (2026-05-27): MEDIUM fixed -- issue #2528. zonal_stats() on dask-backed inputs silently dropped 'majority' from the requested stats list. The mutable default stats_funcs included 'majority' (added in commit 7c8d5759), but the dask path filtered it out at xrspatial/zonal.py:459 (computed_stats = [s for s in stats_funcs.keys() if s in stats_dict]) because 'majority' is not in _DASK_BLOCK_STATS. Symptom: stats(zones=dask, values=dask) returned 7 columns instead of the 8 the docstring promises; stats(..., stats_funcs=['mean','majority']) returned only ['zone','mean'] with no error or warning. Both dask+numpy and dask+cupy were affected (dask+cupy delegates to dask+numpy). Fix: replaced the mutable list literal default with stats_funcs=None and resolved the default per backend inside the function -- numpy/cupy get the full 8-stat list, dask gets the 7-stat subset (no majority).  Explicit majority on dask now raises ValueError with a clear supported-stats message instead of silently filtering. 4 regression tests in test_zonal.py: explicit majority raises on dask, bare default omits majority on dask, bare default keeps majority on numpy, default list is not mutated across calls (covers the historical mutable-default pitfall). All 129 test_zonal.py tests pass (125 pre-existing + 4 new); test_dasymetric.py 61 tests still pass (dasymetric uses zonal.stats internally). Categories: Cat 5 (backend inconsistency: numpy/cupy honoured majority; dask paths silently dropped it). | Pass 1 (2026-03-30T12:00:00Z): historical entry #1090."

xrspatial/tests/test_visibility.py

@@ -2,6 +2,7 @@
 import pytest
 import xarray as xr
 
+from xrspatial.gpu_rtx import has_rtx
 from xrspatial.visibility import _bresenham_line, _extract_transect
 
 
@@ -190,6 +191,29 @@ def test_single_observer_matches_viewshed(self):
         expected = (vs.values != INVISIBLE).astype(np.int32)
         np.testing.assert_array_equal(result.values, expected)
 
+    @pytest.mark.skipif(not has_rtx(),
+                        reason="requires rtxpy for the GPU viewshed path")
+    def test_cupy_does_not_raise(self):
+        """cupy-backed input must accumulate without a numpy/cupy type error.
+
+        The GPU viewshed uses RTX ray tracing while the numpy backend uses
+        the GRASS sweep, so per-cell counts can differ (see viewshed's
+        docstring). This only asserts the accumulator runs and returns valid
+        counts; it does not require cell-for-cell parity with numpy.
+        """
+        import cupy as cp
+        data = np.random.RandomState(3).rand(15, 15).astype(float) * 100
+        raster_cp = _make_raster(data)
+        raster_cp.data = cp.asarray(data)
+        observers = [
+            {'x': 3.0, 'y': 3.0, 'observer_elev': 50},
+            {'x': 12.0, 'y': 12.0, 'observer_elev': 50},
+        ]
+        result = cumulative_viewshed(raster_cp, observers)
+        assert result.dtype == np.int32
+        assert result.values.min() >= 0
+        assert result.values.max() <= len(observers)
+
     def test_wall_blocks_one_side(self):
         """A tall wall blocks visibility from the other side."""
         data = np.zeros((5, 11), dtype=float)

xrspatial/visibility.py

@@ -282,6 +282,12 @@ def cumulative_viewshed(
                       target_elev=te, max_distance=md)
 
         vs_data = vs.data
+        # viewshed() returns a cupy-backed result for a cupy-backed raster.
+        # The numpy `count` accumulator can't be added to a cupy array, so
+        # pull it to numpy first (matching how _extract_transect handles
+        # cupy data).
+        if has_cuda_and_cupy() and is_cupy_array(vs_data):
+            vs_data = vs_data.get()
         if _is_dask and not isinstance(vs_data, da.Array):
             vs_data = da.from_array(vs_data, chunks=raster.data.chunks)
         count = count + (vs_data != INVISIBLE).astype(np.int32)