Fix target_elev contaminating horizon in dask viewshed distance sweep (#1098) (#1099)

_sweep_ring used a single gradient array that included target_elev for
both the visibility test and the horizon update.  The horizon profile
should only reflect raw terrain gradients — target_elev is the height
offset of the thing you're looking for, not extra height on every cell
that might block your view.

Split the gradient into two arrays: gradients (with target_elev, for
visibility) and terrain_gradients (without, for horizon update).

Author: brendancol

xrspatial/tests/test_viewshed.py

@@ -362,3 +362,58 @@ def test_viewshed_max_distance_matches_full(backend):
                 dist_vals[r, c], full_vals[r, c],
                 atol=0.03,
                 err_msg=f"Mismatch at ({r},{c})")
+
+
+def test_viewshed_dask_distance_sweep_target_elev():
+    """Tier C distance sweep must not let target_elev contaminate the
+    horizon profile.  Regression test: previously the horizon was updated
+    with gradients that included target_elev, so intervening terrain
+    appeared artificially tall and caused false occlusion."""
+    from unittest.mock import patch
+
+    ny, nx = 15, 20
+    terrain = np.zeros((ny, nx))
+    # Place a ridge at column 12 that is tall enough to block line-of-sight
+    # to column 15 when target_elev=0, but not when target_elev=8.
+    terrain[:, 12] = 6.0
+
+    xs = np.arange(nx, dtype=float)
+    ys = np.arange(ny, dtype=float)
+
+    obs_x, obs_y = 5.0, 7.0
+    obs_elev = 5
+
+    # Numpy reference: with target_elev=8 the cells behind the ridge at
+    # col 15 should be visible (target sticks up above the ridge).
+    raster_np = xa.DataArray(terrain.copy(), coords=dict(x=xs, y=ys),
+                             dims=["y", "x"])
+    v_np = viewshed(raster_np, x=obs_x, y=obs_y,
+                    observer_elev=obs_elev, target_elev=8)
+
+    # Sanity: numpy says the cell behind the ridge IS visible.
+    assert v_np.values[7, 15] > INVISIBLE, (
+        "numpy reference should see cell (7,15) with target_elev=8")
+
+    # Dask Tier C with same parameters
+    raster_da = xa.DataArray(
+        da.from_array(terrain.copy(), chunks=(5, 5)),
+        coords=dict(x=xs, y=ys), dims=["y", "x"])
+
+    with patch('xrspatial.viewshed._available_memory_bytes',
+               return_value=10_000):
+        v_dask = viewshed(raster_da, x=obs_x, y=obs_y,
+                          observer_elev=obs_elev, target_elev=8)
+
+    result = v_dask.values
+
+    # The cell behind the ridge must be visible in Tier C too.
+    assert result[7, 15] > INVISIBLE, (
+        "Tier C distance sweep should see cell (7,15) with target_elev=8 "
+        "(horizon must not include target_elev)")
+
+    # Visible cells should have angles that roughly match numpy.
+    vis_mask = (v_np.values > INVISIBLE) & (result > INVISIBLE)
+    if vis_mask.any():
+        np.testing.assert_allclose(
+            result[vis_mask], v_np.values[vis_mask], atol=1.0,
+            err_msg="Tier C angles diverge too far from numpy reference")

xrspatial/viewshed.py

@@ -1858,8 +1858,12 @@ def _sweep_ring(rows, cols, elevs, n_cells,
     bin_width = 2.0 * PI / n_angles
     INF = 1e30
 
-    # Pre-compute per-cell values
+    # Pre-compute per-cell values.
+    # Two gradient arrays: one WITH target_elev for visibility testing, one
+    # WITHOUT for the horizon (occlusion) profile.  The terrain itself blocks
+    # the view, not imaginary targets on every cell.
     gradients = np.empty(n_cells, dtype=np.float64)
+    terrain_gradients = np.empty(n_cells, dtype=np.float64)
     center_bins = np.empty(n_cells, dtype=np.int64)
     half_bins_arr = np.empty(n_cells, dtype=np.int64)
     dist_sqs = np.empty(n_cells, dtype=np.float64)
@@ -1887,6 +1891,7 @@ def _sweep_ring(rows, cols, elevs, n_cells,
         eff_elevs[i] = eff_elev
         dist_sqs[i] = dist_sq
         gradients[i] = atan((eff_elev - obs_elev) / dist)
+        terrain_gradients[i] = atan((elev - obs_elev) / dist)
 
         angle = atan2(dy, dx)
         ang_width = cell_size / dist
@@ -1913,13 +1918,13 @@ def _sweep_ring(rows, cols, elevs, n_cells,
             visibility[r, c] = _get_vertical_ang(
                 obs_elev, dist_sqs[i], eff_elevs[i])
 
-    # Pass 2: update horizon with all ring cells
+    # Pass 2: update horizon with raw terrain gradients (no target_elev)
     for i in range(n_cells):
         if valid[i] == 0:
             continue
         cb = center_bins[i]
         hb = half_bins_arr[i]
-        grad = gradients[i]
+        grad = terrain_gradients[i]
         for b in range(-hb, hb + 1):
             idx = (cb + b) % n_angles
             if grad > horizon[idx]: