Aim sky_view_factor rays in ground azimuth space (#3626)

.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.

xrspatial/sky_view_factor.py

@@ -139,12 +139,28 @@ 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)
     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]
@@ -153,9 +169,8 @@ def _svf_cpu(data, max_radius, n_directions, cellsize_x, cellsize_y):
 
             svf_sum = 0.0
             for d in range(n_directions):
-                angle = 2.0 * _pi * d / n_directions
-                dx = _cos(angle)
-                dy = _sin(angle)
+                dx = dxs[d]
+                dy = dys[d]
 
                 max_elev_angle = 0.0
                 for r in range(1, max_radius + 1):
@@ -201,9 +216,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 +329,13 @@ 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.  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
     ----------

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)