Fixes #875: add out-of-core dask CPU viewshed

Add a dask+numpy backend for viewshed() using a three-tier strategy:

A. max_distance specified → extract spatial window, compute, run exact
   R2 angular sweep on the small numpy window.
B. Full R2 fits in memory (280 bytes/cell < 50% available RAM) →
   compute the full dask array, run R2 directly.
C. Otherwise → horizon-profile distance sweep algorithm that processes
   cells in Chebyshev-ring order with an LRU chunk cache, keeping
   memory bounded regardless of grid size.

Also adds a max_distance parameter to viewshed() for all backends,
a memory guard for Tier C, and four new dask-specific tests.

Author: brendancol

xrspatial/tests/test_viewshed.py

@@ -1,3 +1,4 @@
+import dask.array as da
 import datashader as ds
 import numpy as np
 import pandas as pd
@@ -6,6 +7,7 @@
 
 from xrspatial import viewshed
 from xrspatial.tests.general_checks import general_output_checks
+from xrspatial.viewshed import INVISIBLE
 
 from ..gpu_rtx import has_rtx
 
@@ -159,3 +161,103 @@ def test_viewshed_flat(backend, observer_elev, target_elev):
         # Should do better with viewshed gpu output angles.
         mask = (v.data < 90)
         np.testing.assert_allclose(v.data[mask], angle[mask], atol=0.03)
+
+
+# -------------------------------------------------------------------
+# Dask backend tests
+# -------------------------------------------------------------------
+
+@pytest.mark.parametrize("observer_elev", [5, 2])
+@pytest.mark.parametrize("target_elev", [0, 1])
+def test_viewshed_dask_flat(observer_elev, target_elev):
+    """Flat terrain: dask should match analytical formula."""
+    x, y = 0, 0
+    ny, nx = 5, 4
+    arr = da.from_array(np.full((ny, nx), 1.3), chunks=(3, 2))
+    xs = np.arange(nx) * 0.5
+    ys = np.arange(ny) * 1.5
+    xarr = xa.DataArray(arr, coords=dict(x=xs, y=ys), dims=["y", "x"])
+    v = viewshed(xarr, x=x, y=y,
+                 observer_elev=observer_elev, target_elev=target_elev)
+    result = v.values
+
+    xs2, ys2 = np.meshgrid(xs, ys)
+    d_vert = observer_elev - target_elev
+    d_horz = np.sqrt((xs2 - x) ** 2 + (ys2 - y) ** 2)
+    expected = np.rad2deg(np.arctan2(d_horz, d_vert))
+    expected[0, 0] = result[0, 0]  # observer cell
+    np.testing.assert_allclose(result, expected)
+
+
+def test_viewshed_dask_matches_numpy():
+    """Dask should closely match numpy R2 sweep on small varied terrain."""
+    np.random.seed(42)
+    ny, nx = 12, 10
+    terrain = np.random.uniform(0, 5, (ny, nx))
+    xs = np.arange(nx, dtype=float)
+    ys = np.arange(ny, dtype=float)
+
+    # numpy reference
+    raster_np = xa.DataArray(terrain.copy(), coords=dict(x=xs, y=ys),
+                             dims=["y", "x"])
+    v_np = viewshed(raster_np, x=5.0, y=6.0, observer_elev=10)
+
+    # dask (will use Tier B — full compute + R2 — on this small grid)
+    raster_da = xa.DataArray(
+        da.from_array(terrain.copy(), chunks=(4, 5)),
+        coords=dict(x=xs, y=ys), dims=["y", "x"])
+    v_da = viewshed(raster_da, x=5.0, y=6.0, observer_elev=10)
+
+    np.testing.assert_allclose(v_da.values, v_np.values)
+
+
+def test_viewshed_dask_max_distance():
+    """max_distance should produce partial viewshed within radius."""
+    ny, nx = 20, 20
+    arr_np = np.zeros((ny, nx))
+    xs = np.arange(nx, dtype=float)
+    ys = np.arange(ny, dtype=float)
+
+    raster_da = xa.DataArray(
+        da.from_array(arr_np, chunks=(10, 10)),
+        coords=dict(x=xs, y=ys), dims=["y", "x"])
+    v = viewshed(raster_da, x=10.0, y=10.0,
+                 observer_elev=5, max_distance=5.0)
+    result = v.values
+
+    # Observer cell is visible
+    assert result[10, 10] == 180.0
+
+    # Cells far beyond max_distance should be INVISIBLE
+    assert result[0, 0] == INVISIBLE
+    assert result[19, 19] == INVISIBLE
+
+    # Cells within max_distance should be visible (flat terrain, observer up)
+    assert result[10, 12] > INVISIBLE  # 2 cells away
+    assert result[8, 10] > INVISIBLE   # 2 cells away
+
+
+def test_viewshed_dask_distance_sweep():
+    """Force the distance-sweep path (Tier C) and verify flat terrain."""
+    from unittest.mock import patch
+
+    ny, nx = 10, 10
+    arr_np = np.full((ny, nx), 0.0)
+    xs = np.arange(nx, dtype=float)
+    ys = np.arange(ny, dtype=float)
+
+    raster_da = xa.DataArray(
+        da.from_array(arr_np, chunks=(5, 5)),
+        coords=dict(x=xs, y=ys), dims=["y", "x"])
+
+    # R2 needs 280*10*10=28000 bytes; Tier B requires <50% of avail.
+    # Output grid is 10*10*8=800 bytes; memory guard requires <80% of avail.
+    # 10000 bytes: skips Tier B (28000 > 5000) and passes guard (800 < 8000).
+    with patch('xrspatial.viewshed._available_memory_bytes',
+               return_value=10_000):
+        v = viewshed(raster_da, x=5.0, y=5.0, observer_elev=5)
+
+    result = v.values
+    assert result[5, 5] == 180.0
+    # All cells on flat terrain should be visible
+    assert (result > INVISIBLE).all()