diff --git a/.claude/sweep-accuracy-state.csv b/.claude/sweep-accuracy-state.csv index 9f917411..8a7762e3 100644 --- a/.claude/sweep-accuracy-state.csv +++ b/.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." diff --git a/xrspatial/sky_view_factor.py b/xrspatial/sky_view_factor.py index 207c5c1b..acd32c34 100644 --- a/xrspatial/sky_view_factor.py +++ b/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, diff --git a/xrspatial/tests/test_sky_view_factor.py b/xrspatial/tests/test_sky_view_factor.py index f0b39b84..05f53c67 100644 --- a/xrspatial/tests/test_sky_view_factor.py +++ b/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 # ---------------------------------------------------------------------------