Skip to content

sky_view_factor: ray azimuths uniform in cell-index space, biasing SVF on rectangular cells #3626

Description

@brendancol

Describe the bug

sky_view_factor() spaces its n_directions rays evenly in cell-index angle: angle = 2 * pi * d / n_directions, then dx = cos(angle), dy = sin(angle) in units of cells (CPU kernel around line 156, GPU kernel around line 204). Since #1410 the horizon angle along each ray uses true ground distance, but the ray directions are still cell-index directions. With square cells the two agree. With rectangular cells (cellsize_x != cellsize_y) a direction set that is uniform in cell space is not uniform on the ground, so the equal-weight mean over directions oversamples some ground azimuths and undersamples others. The SVF comes out biased.

The docstring says rays are cast at "evenly spaced azimuths", which a user will read as ground azimuths. Anisotropic cell sizes are a supported input; test_anisotropic_cellsize exercises them.

Reproduction

Same physical scene, different grids: a planar ramp rising east at 45 degrees (analytic SVF = 0.75 away from edges), sampled at several cell aspect ratios.

import numpy as np
import xarray as xr
from xrspatial import sky_view_factor

def svf_center(cx, cy):
    ny = nx = 121
    cols = np.arange(nx, dtype=np.float64)
    z = np.broadcast_to(cols * cx, (ny, nx)).copy()  # 45-deg ramp on the ground
    agg = xr.DataArray(z, dims=['y', 'x'])
    out = sky_view_factor(agg, max_radius=40, n_directions=64,
                          cellsize_x=cx, cellsize_y=cy)
    return float(out.data[60, 60])

for cx, cy in [(1.0, 1.0), (1.0, 2.0), (2.0, 1.0), (1.0, 4.0)]:
    print((cx, cy), round(svf_center(cx, cy), 4))

Output on main (b9ba8a8):

(1.0, 1.0) 0.7284
(1.0, 2.0) 0.7731
(2.0, 1.0) 0.6927
(1.0, 4.0) 0.8147

Cell aspect alone moves the result by up to 0.12 (16% relative) on an identical surface. The square-cell value sits below 0.75 too, but that offset comes from nearest-cell ray sampling (the max along a rounded ray picks up azimuth jitter) and applies equally at every aspect ratio. The spread between the rows is the azimuth-weighting bias. The weighting effect can also be computed without any raster, from exact horizon angles: equal weights at ground azimuths atan2(cy*sin(t), cx*cos(t)) give +0.052 error at 2:1 and +0.105 at 4:1.

Expected behavior

The same terrain gives roughly the same SVF whatever the grid's cell aspect ratio. Rays should be spaced evenly in ground azimuth: pick phi = 2 * pi * d / n_directions, convert to cell space as (cos(phi) / cellsize_x, sin(phi) / cellsize_y), and renormalize to a unit cell-space step. For square cells this reduces to the current directions, so results on square grids do not change.

Environment

  • xarray-spatial main (b9ba8a8), Linux, Python 3.14, numba CPU + CUDA

Additional context

Found by the accuracy sweep. All four backends share the direction code, so all four are affected the same way. Follow-up to #1407, which fixed the ground-distance half of anisotropy handling but not the direction spacing.

Metadata

Metadata

Assignees

No one assigned

    Labels

    bugSomething isn't working

    Type

    No type

    Fields

    No fields configured for issues without a type.

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions