Add sky-view factor to surface module

[brendancol]

Author of Proposal:

Reason or Problem

Sky-view factor (SVF) measures what fraction of the sky hemisphere is visible from each cell, from 0 (fully obstructed) to 1 (flat open terrain). It comes up in LiDAR archaeology for picking out subtle terrain features under canopy (Zakek et al. 2011), urban heat island work where you need to quantify street canyon geometry, solar energy modeling as a proxy for diffuse irradiance, and terrain visualization where you want illumination without the directional bias of hillshade.

Right now there's no GDAL-free Python library that does SVF with Dask chunked processing or GPU acceleration. The options are SAGA GIS, WhiteboxTools, or GRASS r.skyview, all C-based and single-threaded.

Proposal

Design:
For each cell, cast rays at n_directions evenly spaced azimuths (default 16). Along each ray, find the maximum elevation angle to the horizon. SVF is then:

SVF(i,j) = 1 - mean(sin(max_horizon_angle(d)) for d in directions)

This is a focal operation with a circular search radius max_radius (in cells). The implementation follows the existing pattern in focal.py:

• NumPy: @ngjit kernel iterating over cells and directions
• CuPy: @cuda.jit kernel with one thread per cell, loop over directions inside the kernel
• Dask: map_overlap with depth=max_radius and boundary=np.nan
• Dask+CuPy: map_overlap dispatching to the CuPy kernel per chunk

New file: xrspatial/sky_view_factor.py

Usage:

from xrspatial import sky_view_factor

svf = sky_view_factor(dem, max_radius=100, n_directions=16)

Value: Nothing in the GDAL-free Python ecosystem does this with chunked or GPU processing. The per-cell computation is independent, so it maps well to map_overlap tiling and CUDA threads, which would make large LiDAR DEMs (100M+ cells) practical without dropping into C tools.

Stakeholders and Impacts

LiDAR analysts, urban climate researchers, solar energy modelers. Adds one new public function and source file. No changes to existing code.

Drawbacks

• The naive algorithm is O(cells * directions * radius), which gets expensive at large radii. A horizon-tracking approach (Yokoyama et al. 2002) cuts this down but is more complex to implement.
• SVF is sensitive to DEM resolution and max_radius choice. The docs will need some guidance on parameter selection.

Alternatives

• Wrap an existing C library (SAGA, WhiteboxTools). Breaks the "no C dependency" design goal.
• Approximate SVF from slope and aspect alone (Dozier & Frew 1990). Faster but much less accurate in complex terrain.
• Reuse the existing viewshed function at multiple observer directions. Would work but is inefficient since viewshed computes more than SVF needs.

Unresolved Questions

• Should the function accept a cellsize parameter for non-square pixels, or require the user to pre-scale?
• Should it support geodesic coordinates (lat/lon DEMs), or require projected CRS?
• Worth exposing the per-direction horizon angles as an optional output? They're useful for anisotropic sky radiance models.