From 1da317a5f5301242412a6b23f08df7a2c90fab36 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Fri, 29 May 2026 16:42:18 -0700 Subject: [PATCH 1/3] viewshed: size max_distance window per axis (#2691) _viewshed_windowed() computed a single cell radius from the coarser of ew_res/ns_res and applied it to both rows and columns. On anisotropic rasters this under-sized the window along the finer axis and clipped cells within max_distance to INVISIBLE. Use radius_rows from ns_res and radius_cols from ew_res. Adds a cross-backend regression test (numpy, cupy, dask). --- xrspatial/tests/test_viewshed.py | 60 ++++++++++++++++++++++++++++++++ xrspatial/viewshed.py | 18 ++++++---- 2 files changed, 71 insertions(+), 7 deletions(-) diff --git a/xrspatial/tests/test_viewshed.py b/xrspatial/tests/test_viewshed.py index c821a90e5..29658e374 100644 --- a/xrspatial/tests/test_viewshed.py +++ b/xrspatial/tests/test_viewshed.py @@ -459,3 +459,63 @@ def test_viewshed_cpu_memory_guard_passes_with_max_distance(): v = viewshed(raster, x=50.0, y=50.0, observer_elev=5, max_distance=3.0) assert v.values[50, 50] == 180.0 + + +@pytest.mark.parametrize("backend", ["numpy", "cupy", "dask"]) +def test_viewshed_max_distance_anisotropic(backend): + """max_distance must not clip cells on anisotropic-resolution rasters. + + Regression test: the windowed path sized its analysis window from the + coarser of ew_res / ns_res and used that single radius for both axes, + so cells within max_distance along the finer axis were dropped from + the window and returned INVISIBLE. + """ + if backend == "cupy": + if not has_rtx(): + pytest.skip("rtxpy not available") + import cupy as cp + + # ew_res = 1 (x), ns_res = 10 (y): strongly anisotropic. + # Use a uniform nonzero elevation: flat enough that all in-range cells + # stay visible, but the RTX mesh builder needs positive max elevation. + ny, nx = 21, 21 + terrain = np.full((ny, nx), 1.0) + xs = np.arange(nx, dtype=float) * 1.0 + ys = np.arange(ny, dtype=float) * 10.0 + + obs_x, obs_y = xs[10], ys[10] + obs_elev = 50 + + base = xa.DataArray(terrain, coords=dict(x=xs, y=ys), dims=["y", "x"]) + full = viewshed(base, x=obs_x, y=obs_y, observer_elev=obs_elev) + full_vals = full.values + + # max_distance=8 reaches 8 columns east (ew_res=1) but only 0.8 rows. + arr = terrain.copy() + if backend == "cupy": + arr = cp.asarray(arr) + raster = xa.DataArray(arr, coords=dict(x=xs, y=ys), dims=["y", "x"]) + if backend == "dask": + raster = xa.DataArray( + da.from_array(terrain.copy(), chunks=(7, 7)), + coords=dict(x=xs, y=ys), dims=["y", "x"]) + + v = viewshed(raster, x=obs_x, y=obs_y, observer_elev=obs_elev, + max_distance=8.0) + if backend == "cupy": + result = v.data.get() + else: + # numpy and dask both resolve through .values + result = v.values + + # Cells within max_distance along the finer (x) axis must be evaluated + # and match the full viewshed, not be clipped to INVISIBLE. + for c in (12, 15, 17): # 2, 5, 7 units east — all < 8 + dist = abs(xs[c] - obs_x) + assert dist < 8.0 + assert result[10, c] > INVISIBLE, ( + f"cell (10,{c}) at distance {dist} wrongly clipped") + np.testing.assert_allclose(result[10, c], full_vals[10, c], atol=0.03) + + # Cells outside max_distance stay INVISIBLE. + assert result[10, 19] == INVISIBLE # 9 units east, > 8 diff --git a/xrspatial/viewshed.py b/xrspatial/viewshed.py index b5d210192..5288eaeb9 100644 --- a/xrspatial/viewshed.py +++ b/xrspatial/viewshed.py @@ -2077,13 +2077,17 @@ def _viewshed_windowed(raster, x, y, observer_elev, target_elev, x_range = (x_coords[0], x_coords[-1]) ew_res = (x_range[1] - x_range[0]) / (width - 1) if width > 1 else 1.0 ns_res = (y_range[1] - y_range[0]) / (height - 1) if height > 1 else 1.0 - cell_size = max(abs(ew_res), abs(ns_res)) - radius_cells = int(np.ceil(max_distance / cell_size)) - - r_lo = max(0, obs_r - radius_cells) - r_hi = min(height, obs_r + radius_cells + 1) - c_lo = max(0, obs_c - radius_cells) - c_hi = min(width, obs_c + radius_cells + 1) + # Size the window per axis: rows are spaced by ns_res, columns by ew_res. + # Using a single radius from the coarser resolution under-sizes the + # window along the finer axis and clips cells that are within + # max_distance there. + radius_rows = int(np.ceil(max_distance / abs(ns_res))) + radius_cols = int(np.ceil(max_distance / abs(ew_res))) + + r_lo = max(0, obs_r - radius_rows) + r_hi = min(height, obs_r + radius_rows + 1) + c_lo = max(0, obs_c - radius_cols) + c_hi = min(width, obs_c + radius_cols + 1) window = raster.isel(y=slice(r_lo, r_hi), x=slice(c_lo, c_hi)) From d1013328864cbc56eaf53303da060a58f4c43417 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Fri, 29 May 2026 16:44:38 -0700 Subject: [PATCH 2/3] Address review nit: cover both anisotropy orientations (#2691) --- xrspatial/tests/test_viewshed.py | 45 ++++++++++++++++++++------------ 1 file changed, 28 insertions(+), 17 deletions(-) diff --git a/xrspatial/tests/test_viewshed.py b/xrspatial/tests/test_viewshed.py index 29658e374..bec9d19bf 100644 --- a/xrspatial/tests/test_viewshed.py +++ b/xrspatial/tests/test_viewshed.py @@ -462,26 +462,32 @@ def test_viewshed_cpu_memory_guard_passes_with_max_distance(): @pytest.mark.parametrize("backend", ["numpy", "cupy", "dask"]) -def test_viewshed_max_distance_anisotropic(backend): +@pytest.mark.parametrize("fine_axis", ["x", "y"]) +def test_viewshed_max_distance_anisotropic(backend, fine_axis): """max_distance must not clip cells on anisotropic-resolution rasters. Regression test: the windowed path sized its analysis window from the coarser of ew_res / ns_res and used that single radius for both axes, so cells within max_distance along the finer axis were dropped from - the window and returned INVISIBLE. + the window and returned INVISIBLE. Both anisotropy orientations are + checked so a mix-up between radius_rows / radius_cols is caught. """ if backend == "cupy": if not has_rtx(): pytest.skip("rtxpy not available") import cupy as cp - # ew_res = 1 (x), ns_res = 10 (y): strongly anisotropic. + # Strongly anisotropic: one axis spaced by 1, the other by 10. # Use a uniform nonzero elevation: flat enough that all in-range cells # stay visible, but the RTX mesh builder needs positive max elevation. ny, nx = 21, 21 terrain = np.full((ny, nx), 1.0) - xs = np.arange(nx, dtype=float) * 1.0 - ys = np.arange(ny, dtype=float) * 10.0 + if fine_axis == "x": + xs = np.arange(nx, dtype=float) * 1.0 # fine + ys = np.arange(ny, dtype=float) * 10.0 # coarse + else: + xs = np.arange(nx, dtype=float) * 10.0 # coarse + ys = np.arange(ny, dtype=float) * 1.0 # fine obs_x, obs_y = xs[10], ys[10] obs_elev = 50 @@ -490,7 +496,8 @@ def test_viewshed_max_distance_anisotropic(backend): full = viewshed(base, x=obs_x, y=obs_y, observer_elev=obs_elev) full_vals = full.values - # max_distance=8 reaches 8 columns east (ew_res=1) but only 0.8 rows. + # max_distance=8 reaches 8 cells along the fine axis but < 1 along the + # coarse axis, so the buggy single-radius window clipped the fine axis. arr = terrain.copy() if backend == "cupy": arr = cp.asarray(arr) @@ -508,14 +515,18 @@ def test_viewshed_max_distance_anisotropic(backend): # numpy and dask both resolve through .values result = v.values - # Cells within max_distance along the finer (x) axis must be evaluated - # and match the full viewshed, not be clipped to INVISIBLE. - for c in (12, 15, 17): # 2, 5, 7 units east — all < 8 - dist = abs(xs[c] - obs_x) - assert dist < 8.0 - assert result[10, c] > INVISIBLE, ( - f"cell (10,{c}) at distance {dist} wrongly clipped") - np.testing.assert_allclose(result[10, c], full_vals[10, c], atol=0.03) - - # Cells outside max_distance stay INVISIBLE. - assert result[10, 19] == INVISIBLE # 9 units east, > 8 + # Cells within max_distance along the fine axis must be evaluated and + # match the full viewshed, not be clipped to INVISIBLE. obs is at + # index 10 on both axes; step along the fine axis from the observer. + for offset in (2, 5, 7): # all < 8 cells along the fine axis + if fine_axis == "x": + rc = (10, 10 + offset) + else: + rc = (10 + offset, 10) + assert result[rc] > INVISIBLE, ( + f"cell {rc} ({offset} fine-axis cells away) wrongly clipped") + np.testing.assert_allclose(result[rc], full_vals[rc], atol=0.03) + + # A cell 9 fine-axis cells away (> max_distance) stays INVISIBLE. + far = (10, 19) if fine_axis == "x" else (19, 10) + assert result[far] == INVISIBLE From 76f41ae852718f6c6625e39c977dfc36c4ef3b30 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Fri, 29 May 2026 17:07:56 -0700 Subject: [PATCH 3/3] sweep-accuracy: record viewshed inspection (issue #2691, PR #2702) --- .claude/sweep-accuracy-state.csv | 1 + 1 file changed, 1 insertion(+) diff --git a/.claude/sweep-accuracy-state.csv b/.claude/sweep-accuracy-state.csv index 4888d6885..d0bc02d83 100644 --- a/.claude/sweep-accuracy-state.csv +++ b/.claude/sweep-accuracy-state.csv @@ -31,6 +31,7 @@ sieve,2026-04-13T12:00:00Z,,,,Union-find CCL correct. NaN excluded from labeling sky_view_factor,2026-05-01,1407,HIGH,4,Horizon angle ignored cell size; fixed by passing cellsize_x/cellsize_y into CPU+GPU kernels and using ground distance 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." 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."