From 1bf097f7e38c1d28dd1426df58006024d53232ec Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Fri, 3 Jul 2026 09:16:20 -0400 Subject: [PATCH 1/2] Aim sky_view_factor rays in ground azimuth space (#3626) Ray directions were uniform in cell-index angle, so rectangular cells (cellsize_x != cellsize_y) sampled ground azimuths non-uniformly and biased the SVF average. Convert each ground azimuth to a unit step in cell-index space instead. Square-cell results are unchanged. Also records the 2026-07-03 accuracy sweep state for the module. --- .claude/sweep-accuracy-state.csv | 2 +- xrspatial/sky_view_factor.py | 29 +++++++++++--- xrspatial/tests/test_sky_view_factor.py | 50 +++++++++++++++++++++++++ 3 files changed, 74 insertions(+), 7 deletions(-) diff --git a/.claude/sweep-accuracy-state.csv b/.claude/sweep-accuracy-state.csv index cc43f3153..595d0f14a 100644 --- a/.claude/sweep-accuracy-state.csv +++ b/.claude/sweep-accuracy-state.csv @@ -32,7 +32,7 @@ rasterize,2026-06-18,3384,HIGH,1;5,dask all_touched polygon-boundary supercover reproject,2026-06-12,3274,HIGH,1;4,"3 confirmed bugs, all kernel-vs-PROJ parity: #3274 HIGH LAEA inverse spurious /rq (2.6 km err for 3035) + _authalic_apa inverse series wrong (4.8 m in AEA/CEA inverses; PROJ 3-term = 1.6 mm), CPU+CUDA kernels both; #3275 HIGH _is_wgs84_compatible_ellipsoid passes R-defined spheres (MODIS sinusoidal 18.9 km err) and _aea_params/_cea_params lack the guard entirely (23.8 km on spherical aea/cea); #3276 MEDIUM itrf helmert scale 1e-9 but PROJ +s is ppm (1e-6), ~23 mm err. Verified clean: merc/emerc/UTM/tmerc/LCC/polar stere (incl lat_ts akm1) forward+inverse <=1e-5 m vs pyproj; resampling kernels NaN handling and GDAL renorm match across numpy/cupy/dask (CUDA run, gpu-vs-cpu 1.3e-7); dask footprint chunk-skip bbox is a superset in all probed cases (no holes). LOW (documented only): _source_footprint_in_target probe array typo uses x-midpoint mx as a latitude in last 3 ys entries (bbox superset, correctness unaffected)." resample,2026-05-29,2610,HIGH,3;5,"dask interp (nearest/bilinear) overlap depth=1 too small on downsample; block-centered source coord landed past chunk, map_coordinates clamped to edge -> wrong seam rows. Fixed PR #2627 via per-axis _downsample_radius. cupy+dask+cupy verified." sieve,2026-04-13T12:00:00Z,,,,Union-find CCL correct. NaN excluded from labeling. All backends funnel through _sieve_numpy. -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 +sky_view_factor,2026-07-03,3626,MEDIUM,4,"Re-audit post-#1407 fix (ground distance OK). MEDIUM #3626: ray azimuths uniform in cell-INDEX angle not ground azimuth, so anisotropic cellsizes bias the azimuthal quadrature (45-deg ramp, true SVF 0.75: 0.728/0.773/0.693 at 1:1/1:2/2:1 aspect; +0.105 pure weighting error at 4:1); fixed by aiming rays in ground space (cos(phi)/cx, sin(phi)/cy normalized), square-cell results bit-identical. Cats 1-3,5 clean: float64 throughout, NaN center->NaN out, flat=1.0 exactly all 4 backends, flip/rot90/translate invariants at machine eps, numpy==dask==cupy==dask+cupy (max 2.2e-16) incl NaN holes straddling chunk boundaries, map_overlap depth==max_radius==stencil radius, depth>chunksize works (dask 2025.7). Cat6 skipped: richdem-unavailable rvt-unavailable grass-unavailable. LOWs (not fixed): NaN neighbor silently truncates ray (horizon beyond nodata ignored, undocumented but keeps numpy/dask edge semantics consistent); nearest-cell ray rounding gives ~-2.9% jitter bias on 45-deg slopes even square cells (inherent to integer-cell marching, reference tools interpolate along ray). For test-coverage sweep: cross-backend tests had no NaN-input and no anisotropic-cellsize cases (added aniso numpy/dask/cupy parity in #3626 PR); for style sweep: unused import _boundary_to_dask (F401)." 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. diff --git a/xrspatial/sky_view_factor.py b/xrspatial/sky_view_factor.py index b48db9b1f..c00cc53eb 100644 --- a/xrspatial/sky_view_factor.py +++ b/xrspatial/sky_view_factor.py @@ -139,7 +139,10 @@ def _svf_cpu(data, max_radius, n_directions, cellsize_x, cellsize_y): *cellsize_x* and *cellsize_y* are the ground distance per cell along the x and y axes, in the same units as the elevation values. They are used to convert the integer ray step *r* into a true ground - distance for the horizon angle calculation. + distance for the horizon angle calculation, and to aim the rays so + the *n_directions* azimuths are evenly spaced on the ground. For + rectangular cells a direction set that is uniform in cell-index + space is not uniform in ground azimuth (issue #3626). """ rows, cols = data.shape out = np.empty((rows, cols), dtype=np.float64) @@ -153,9 +156,15 @@ def _svf_cpu(data, max_radius, n_directions, cellsize_x, cellsize_y): svf_sum = 0.0 for d in range(n_directions): + # Ground azimuth, converted to a unit step in cell-index + # space so the sampled azimuths stay evenly spaced on the + # ground even when cells are rectangular. angle = 2.0 * _pi * d / n_directions - dx = _cos(angle) - dy = _sin(angle) + dx = _cos(angle) / cellsize_x + dy = _sin(angle) / cellsize_y + dnorm = _sqrt(dx * dx + dy * dy) + dx = dx / dnorm + dy = dy / dnorm max_elev_angle = 0.0 for r in range(1, max_radius + 1): @@ -201,9 +210,14 @@ def _svf_gpu(data, out, max_radius, n_directions, cellsize_x, cellsize_y): pi2 = 2.0 * 3.141592653589793 for d in range(n_directions): + # Ground azimuth, converted to a unit step in cell-index space + # (see _svf_cpu). angle = pi2 * d / n_directions - dx = _cos(angle) - dy = _sin(angle) + dx = _cos(angle) / cellsize_x + dy = _sin(angle) / cellsize_y + dnorm = _sqrt(dx * dx + dy * dy) + dx = dx / dnorm + dy = dy / dnorm max_elev_angle = 0.0 for r in range(1, max_radius + 1): @@ -309,7 +323,10 @@ def sky_view_factor( The horizon angle along each ray uses true ground distance (``cell_step * cellsize``), so cell size must be in the same unit - as the elevation values (e.g. meters for both). + as the elevation values (e.g. meters for both). Ray azimuths are + evenly spaced on the ground: with rectangular cells + (``cellsize_x != cellsize_y``) each ray's cell-index direction is + adjusted so the ground azimuth spacing stays uniform. Parameters ---------- diff --git a/xrspatial/tests/test_sky_view_factor.py b/xrspatial/tests/test_sky_view_factor.py index 05f53c676..7ca0a3988 100644 --- a/xrspatial/tests/test_sky_view_factor.py +++ b/xrspatial/tests/test_sky_view_factor.py @@ -346,6 +346,56 @@ def test_anisotropic_cellsize(): ) +def test_anisotropic_cellsize_ground_truth_consistency(): + """The same physical surface sampled on square and rectangular cells + gives about the same SVF (#3626). + + A planar 45-degree ramp has analytic SVF = 0.75 away from edges, so + the value cannot depend on how the grid happens to be cellated. + Before the fix, ray azimuths were uniform in cell-index space, and + the ramp gave 0.728 / 0.773 / 0.693 at cell aspect 1:1 / 1:2 / 2:1. + """ + def svf_center(cx, cy): + ny = nx = 121 + cols = np.arange(nx, dtype=np.float64) + # z rises 1 elevation unit per ground unit along x: 45-deg ramp + z = np.broadcast_to(cols * cx, (ny, nx)).copy() + agg = xr.DataArray(z, dims=['y', 'x']) + out = sky_view_factor(agg, max_radius=40, n_directions=32, + cellsize_x=cx, cellsize_y=cy) + return float(out.data[60, 60]) + + square = svf_center(1.0, 1.0) + tall = svf_center(1.0, 2.0) + wide = svf_center(2.0, 1.0) + assert abs(tall - square) < 0.02, ( + f"1:2 cells diverge from square: {tall} vs {square}") + assert abs(wide - square) < 0.02, ( + f"2:1 cells diverge from square: {wide} vs {square}") + + +@dask_array_available +def test_anisotropic_numpy_equals_dask(): + """Backend parity for rectangular cells (#3626).""" + data = np.random.default_rng(3626).random((30, 30)) * 100 + numpy_agg = create_test_raster(data, backend='numpy') + dask_agg = create_test_raster(data, backend='dask', chunks=(10, 10)) + func = lambda agg: sky_view_factor(agg, max_radius=5, n_directions=16, + cellsize_x=1.0, cellsize_y=3.0) + assert_numpy_equals_dask_numpy(numpy_agg, dask_agg, func, nan_edges=False) + + +@cuda_and_cupy_available +def test_anisotropic_numpy_equals_cupy(): + """Backend parity for rectangular cells (#3626).""" + data = np.random.default_rng(3626).random((30, 30)) * 100 + numpy_agg = create_test_raster(data, backend='numpy') + cupy_agg = create_test_raster(data, backend='cupy') + func = lambda agg: sky_view_factor(agg, max_radius=5, n_directions=16, + cellsize_x=1.0, cellsize_y=3.0) + assert_numpy_equals_cupy(numpy_agg, cupy_agg, func, nan_edges=False) + + def test_flat_surface_unaffected_by_cellsize(): """SVF on flat terrain is 1 regardless of cell size.""" data = np.full((30, 30), 100.0, dtype=np.float64) From 66d805db01f81ee9b6050c3d8183e61b6ebd2551 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Fri, 3 Jul 2026 09:21:30 -0400 Subject: [PATCH 2/2] Address review: hoist CPU direction table, note azimuth quantization limit (#3626) --- xrspatial/sky_view_factor.py | 29 +++++++++++++++++++---------- 1 file changed, 19 insertions(+), 10 deletions(-) diff --git a/xrspatial/sky_view_factor.py b/xrspatial/sky_view_factor.py index c00cc53eb..9dea6e66a 100644 --- a/xrspatial/sky_view_factor.py +++ b/xrspatial/sky_view_factor.py @@ -148,6 +148,19 @@ def _svf_cpu(data, max_radius, n_directions, cellsize_x, cellsize_y): out = np.empty((rows, cols), dtype=np.float64) out[:] = np.nan + # Ground azimuths, converted to unit steps in cell-index space so + # the sampled azimuths stay evenly spaced on the ground even when + # cells are rectangular. Computed once, reused for every pixel. + dxs = np.empty(n_directions, dtype=np.float64) + dys = np.empty(n_directions, dtype=np.float64) + for d in range(n_directions): + angle = 2.0 * _pi * d / n_directions + dx = _cos(angle) / cellsize_x + dy = _sin(angle) / cellsize_y + dnorm = _sqrt(dx * dx + dy * dy) + dxs[d] = dx / dnorm + dys[d] = dy / dnorm + for y in range(rows): for x in range(cols): center = data[y, x] @@ -156,15 +169,8 @@ def _svf_cpu(data, max_radius, n_directions, cellsize_x, cellsize_y): svf_sum = 0.0 for d in range(n_directions): - # Ground azimuth, converted to a unit step in cell-index - # space so the sampled azimuths stay evenly spaced on the - # ground even when cells are rectangular. - angle = 2.0 * _pi * d / n_directions - dx = _cos(angle) / cellsize_x - dy = _sin(angle) / cellsize_y - dnorm = _sqrt(dx * dx + dy * dy) - dx = dx / dnorm - dy = dy / dnorm + dx = dxs[d] + dy = dys[d] max_elev_angle = 0.0 for r in range(1, max_radius + 1): @@ -326,7 +332,10 @@ def sky_view_factor( as the elevation values (e.g. meters for both). Ray azimuths are evenly spaced on the ground: with rectangular cells (``cellsize_x != cellsize_y``) each ray's cell-index direction is - adjusted so the ground azimuth spacing stays uniform. + adjusted so the ground azimuth spacing stays uniform. The spacing + is nominal to within one cell of rounding, since rays sample whole + cells; at extreme aspect ratios (roughly beyond 4:1) the sampled + azimuths still quantize toward the grid axes. Parameters ----------