Use ground distance for sky_view_factor horizon angle (#1407)

.claude/sweep-accuracy-state.csv

@@ -19,6 +19,7 @@ polygonize,2026-04-13T12:00:00Z,1190,,,NaN pixels not masked in numpy/dask backe
 proximity,2026-03-30T15:00:00Z,,,,Direction >= boundary fragile but works due to truncated constant. Float32 truncation is design choice. No wrong-results bugs found.
 resample,2026-04-14T12:00:00Z,1202,,,Interpolation used edge-aligned coords instead of block-centered. Fix in PR #1204.
 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
 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."
 visibility,2026-04-13T12:00:00Z,,,,"Bresenham line, LOS kernel, Fresnel zone all correct. All backends converge to numpy."

xrspatial/sky_view_factor.py

@@ -45,6 +45,7 @@ class cupy:
     _validate_raster,
     _validate_scalar,
     cuda_args,
+    get_dataarray_resolution,
     ngjit,
 )
 
@@ -131,8 +132,14 @@ def _check_gpu_memory(rows, cols):
 # ---------------------------------------------------------------------------
 
 @ngjit
-def _svf_cpu(data, max_radius, n_directions):
-    """Compute SVF over an entire 2-D array on the CPU."""
+def _svf_cpu(data, max_radius, n_directions, cellsize_x, cellsize_y):
+    """Compute SVF over an entire 2-D array on the CPU.
+
+    *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.
+    """
     rows, cols = data.shape
     out = np.empty((rows, cols), dtype=np.float64)
     out[:] = np.nan
@@ -159,7 +166,9 @@ def _svf_cpu(data, max_radius, n_directions):
                     if elev != elev:  # NaN
                         break
                     dz = elev - center
-                    dist = float(r)
+                    gx = (sx - x) * cellsize_x
+                    gy = (sy - y) * cellsize_y
+                    dist = _sqrt(gx * gx + gy * gy)
                     elev_angle = _atan2(dz, dist)
                     if elev_angle > max_elev_angle:
                         max_elev_angle = elev_angle
@@ -175,7 +184,7 @@ def _svf_cpu(data, max_radius, n_directions):
 # ---------------------------------------------------------------------------
 
 @cuda.jit
-def _svf_gpu(data, out, max_radius, n_directions):
+def _svf_gpu(data, out, max_radius, n_directions, cellsize_x, cellsize_y):
     """CUDA global kernel: one thread per cell."""
     y, x = cuda.grid(2)
     rows, cols = data.shape
@@ -205,7 +214,9 @@ def _svf_gpu(data, out, max_radius, n_directions):
             if elev != elev:  # NaN
                 break
             dz = elev - center
-            dist = float(r)
+            gx = (sx - x) * cellsize_x
+            gy = (sy - y) * cellsize_y
+            dist = _sqrt(gx * gx + gy * gy)
             elev_angle = _atan2(dz, dist)
             if elev_angle > max_elev_angle:
                 max_elev_angle = elev_angle
@@ -219,24 +230,32 @@ def _svf_gpu(data, out, max_radius, n_directions):
 # Backend wrappers
 # ---------------------------------------------------------------------------
 
-def _run_numpy(data, max_radius, n_directions):
+def _run_numpy(data, max_radius, n_directions, cellsize_x, cellsize_y):
     _check_memory(data.shape[0], data.shape[1])
     data = data.astype(np.float64)
-    return _svf_cpu(data, max_radius, n_directions)
+    return _svf_cpu(data, max_radius, n_directions, cellsize_x, cellsize_y)
 
 
-def _run_cupy(data, max_radius, n_directions):
+def _run_cupy(data, max_radius, n_directions, cellsize_x, cellsize_y):
     _check_gpu_memory(data.shape[0], data.shape[1])
     data = data.astype(cupy.float64)
     out = cupy.full(data.shape, cupy.nan, dtype=cupy.float64)
     griddim, blockdim = cuda_args(data.shape)
-    _svf_gpu[griddim, blockdim](data, out, max_radius, n_directions)
+    _svf_gpu[griddim, blockdim](
+        data, out, max_radius, n_directions, cellsize_x, cellsize_y
+    )
     return out
 
 
-def _run_dask_numpy(data, max_radius, n_directions):
+def _run_dask_numpy(data, max_radius, n_directions, cellsize_x, cellsize_y):
     data = data.astype(np.float64)
-    _func = partial(_svf_cpu, max_radius=max_radius, n_directions=n_directions)
+    _func = partial(
+        _svf_cpu,
+        max_radius=max_radius,
+        n_directions=n_directions,
+        cellsize_x=cellsize_x,
+        cellsize_y=cellsize_y,
+    )
     out = data.map_overlap(
         _func,
         depth=(max_radius, max_radius),
@@ -246,9 +265,15 @@ def _run_dask_numpy(data, max_radius, n_directions):
     return out
 
 
-def _run_dask_cupy(data, max_radius, n_directions):
+def _run_dask_cupy(data, max_radius, n_directions, cellsize_x, cellsize_y):
     data = data.astype(cupy.float64)
-    _func = partial(_run_cupy, max_radius=max_radius, n_directions=n_directions)
+    _func = partial(
+        _run_cupy,
+        max_radius=max_radius,
+        n_directions=n_directions,
+        cellsize_x=cellsize_x,
+        cellsize_y=cellsize_y,
+    )
     out = data.map_overlap(
         _func,
         depth=(max_radius, max_radius),
@@ -267,6 +292,8 @@ def sky_view_factor(
     agg: xr.DataArray,
     max_radius: int = 10,
     n_directions: int = 16,
+    cellsize_x: Optional[float] = None,
+    cellsize_y: Optional[float] = None,
     name: Optional[str] = 'sky_view_factor',
 ) -> xr.DataArray:
     """Compute the sky-view factor for each cell of a DEM.
@@ -277,6 +304,10 @@ def sky_view_factor(
     *max_radius* cells, and the maximum elevation angle along each
     ray determines the horizon obstruction.
 
+    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).
+
     Parameters
     ----------
     agg : xarray.DataArray or xr.Dataset
@@ -286,10 +317,17 @@ def sky_view_factor(
         data variable independently.
     max_radius : int, default=10
         Maximum search distance in cells along each ray direction.
-        Cells within *max_radius* of the raster edge will be NaN.
     n_directions : int, default=16
         Number of azimuth directions to sample, evenly spaced
         around 360 degrees.
+    cellsize_x : float, optional
+        Ground distance per cell along the x axis, in the same unit
+        as the elevation values.  If not provided, it is read from
+        ``agg.attrs['res']`` or computed from the x coordinate.
+    cellsize_y : float, optional
+        Ground distance per cell along the y axis, in the same unit
+        as the elevation values.  If not provided, it is read from
+        ``agg.attrs['res']`` or computed from the y coordinate.
     name : str, default='sky_view_factor'
         Name of the output DataArray.
 
@@ -316,13 +354,33 @@ def sky_view_factor(
     _validate_scalar(n_directions, func_name='sky_view_factor',
                      name='n_directions', dtype=int, min_val=1)
 
+    if cellsize_x is None or cellsize_y is None:
+        try:
+            res_x, res_y = get_dataarray_resolution(agg)
+        except Exception:
+            res_x, res_y = 1.0, 1.0
+        if cellsize_x is None:
+            cellsize_x = float(abs(res_x)) if res_x else 1.0
+        if cellsize_y is None:
+            cellsize_y = float(abs(res_y)) if res_y else 1.0
+    else:
+        cellsize_x = float(abs(cellsize_x))
+        cellsize_y = float(abs(cellsize_y))
+
+    if cellsize_x <= 0:
+        cellsize_x = 1.0
+    if cellsize_y <= 0:
+        cellsize_y = 1.0
+
     mapper = ArrayTypeFunctionMapping(
         numpy_func=_run_numpy,
         cupy_func=_run_cupy,
         dask_func=_run_dask_numpy,
         dask_cupy_func=_run_dask_cupy,
     )
-    out = mapper(agg)(agg.data, max_radius, n_directions)
+    out = mapper(agg)(
+        agg.data, max_radius, n_directions, cellsize_x, cellsize_y
+    )
     return xr.DataArray(
         out,
         name=name,

xrspatial/tests/test_sky_view_factor.py

@@ -230,14 +230,135 @@ def test_known_svf_simple_ramp():
     rows, cols = 30, 30
     x = np.arange(cols, dtype=np.float64)
     data = np.broadcast_to(x * 1.0, (rows, cols)).copy()
-    agg = create_test_raster(data)
-    result = sky_view_factor(agg, max_radius=5, n_directions=16)
+    agg = create_test_raster(data, attrs={'res': (1.0, 1.0)})
+    result = sky_view_factor(agg, max_radius=5, n_directions=16,
+                             cellsize_x=1.0, cellsize_y=1.0)
     center = result.data[15, 15]
     assert np.isfinite(center)
     assert center < 1.0, f"Expected SVF < 1 on tilted plane, got {center}"
     assert center > 0.5, f"Expected SVF > 0.5 on gentle tilt, got {center}"
 
 
+# ---------------------------------------------------------------------------
+# Cell size correctness
+# ---------------------------------------------------------------------------
+
+def test_cellsize_changes_horizon_angle():
+    """Halving the cell size should double the horizon angle.
+
+    For a uniform z = x ramp with rise-per-cell = 1:
+      cell_size=1.0 -> horizon angle = atan(1/1)   = 45 deg, sin = 0.7071
+      cell_size=0.5 -> horizon angle = atan(1/0.5) = 63.43 deg, sin = 0.8944
+
+    The buggy code used dist=r (cells) and produced 45 deg in both
+    cases, so SVF was the same regardless of cell size.
+    """
+    rows, cols = 30, 30
+    x = np.arange(cols, dtype=np.float64)
+    data = np.broadcast_to(x, (rows, cols)).copy()
+
+    agg_unit = create_test_raster(data, attrs={'res': (1.0, 1.0)})
+    agg_half = create_test_raster(data, attrs={'res': (0.5, 0.5)})
+
+    res_unit = sky_view_factor(agg_unit, max_radius=5, n_directions=16)
+    res_half = sky_view_factor(agg_half, max_radius=5, n_directions=16)
+
+    # Smaller cell size -> steeper horizon -> lower SVF
+    assert res_half.data[15, 15] < res_unit.data[15, 15] - 0.05, (
+        "Expected SVF to drop noticeably when cell size shrinks; "
+        f"got unit={res_unit.data[15, 15]} half={res_half.data[15, 15]}"
+    )
+
+
+def test_cellsize_explicit_override():
+    """Explicit cellsize_x / cellsize_y override the DataArray attrs."""
+    rows, cols = 30, 30
+    x = np.arange(cols, dtype=np.float64)
+    data = np.broadcast_to(x, (rows, cols)).copy()
+
+    # DataArray says res=1.0 but we override to 0.5
+    agg = create_test_raster(data, attrs={'res': (1.0, 1.0)})
+    res_default = sky_view_factor(agg, max_radius=5, n_directions=16)
+    res_override = sky_view_factor(
+        agg, max_radius=5, n_directions=16, cellsize_x=0.5, cellsize_y=0.5,
+    )
+
+    assert res_override.data[15, 15] < res_default.data[15, 15] - 0.05, (
+        "Explicit cellsize override should change the result; "
+        f"got default={res_default.data[15, 15]} "
+        f"override={res_override.data[15, 15]}"
+    )
+
+
+def test_cellsize_45_degree_ramp_value():
+    """A z=x ramp with cellsize=1 should produce a horizon angle of
+    exactly 45 degrees in the uphill direction.
+
+    SVF = 1 - mean_d sin(max_horizon_angle_d).  For 16 directions on
+    a uniform planar ramp, the contribution along the uphill ray
+    dominates and hits sin(45) = sqrt(2)/2 ~= 0.7071.  A coarse
+    sanity check: SVF should sit near (1 - 0.7071/16) = 0.956 to
+    well below 1, but never above the no-correction value when cell
+    size is unit, and below it when cell size is sub-unit.
+    """
+    rows, cols = 30, 30
+    x = np.arange(cols, dtype=np.float64)
+    data = np.broadcast_to(x, (rows, cols)).copy()
+    agg = create_test_raster(data, attrs={'res': (1.0, 1.0)})
+    result = sky_view_factor(agg, max_radius=5, n_directions=16,
+                             cellsize_x=1.0, cellsize_y=1.0)
+    center = result.data[15, 15]
+    # 45 deg horizon means sin = 0.7071; with 16 directions and
+    # symmetric horizons, SVF stays in a plausible band.
+    assert 0.55 < center < 0.95, (
+        f"Expected SVF in (0.55, 0.95) for 45-deg ramp, got {center}"
+    )
+
+
+def test_anisotropic_cellsize():
+    """Anisotropic cell sizes (cellsize_x != cellsize_y).
+
+    Build a uniform z = x ramp (gradient along x only) and probe an
+    interior cell.  Scaling cellsize_x changes the ground distance to
+    the uphill horizon, which must change the horizon angle and SVF.
+    cellsize_y does not affect a ramp that is constant in y, so we
+    isolate the x-axis scaling.
+    """
+    rows, cols = 30, 30
+    x = np.arange(cols, dtype=np.float64)
+    data = np.broadcast_to(x, (rows, cols)).copy()
+
+    agg = create_test_raster(data, attrs={'res': (1.0, 1.0)})
+    res_wide_x = sky_view_factor(
+        agg, max_radius=5, n_directions=16,
+        cellsize_x=10.0, cellsize_y=1.0,
+    )
+    res_narrow_x = sky_view_factor(
+        agg, max_radius=5, n_directions=16,
+        cellsize_x=0.1, cellsize_y=1.0,
+    )
+    # Larger cellsize_x -> smaller horizon angle -> larger SVF
+    assert (res_wide_x.data[15, 15]
+            > res_narrow_x.data[15, 15] + 0.05), (
+        "Wider cellsize_x should give larger SVF; "
+        f"got wide={res_wide_x.data[15, 15]} "
+        f"narrow={res_narrow_x.data[15, 15]}"
+    )
+
+
+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)
+    for csx, csy in [(0.5, 0.5), (1.0, 1.0), (30.0, 30.0)]:
+        agg = create_test_raster(data, attrs={'res': (csx, csy)})
+        result = sky_view_factor(agg, max_radius=5, n_directions=16)
+        interior = result.data[5:-5, 5:-5]
+        np.testing.assert_allclose(
+            interior, 1.0, atol=1e-10,
+            err_msg=f"flat terrain at cellsize=({csx}, {csy})",
+        )
+
+
 # ---------------------------------------------------------------------------
 # Memory guard
 # ---------------------------------------------------------------------------