Add solar radiation modeling over DEMs

Start by locating the existing viewshed infrastructure and its xarray DataArray interfaces, then assess whether horizon tracing can be reused. Review the dask map_overlap and GPU requirements before defining the design. Done means direct, diffuse, and terrain-reflected radiation are separately available, with hourly, daily, and annual integration over terrain and configurable atmospheric transmittance.

Sky view factor exists, but there's no actual solar irradiance model. Solar radiation modeling over terrain is a common need in renewable energy siting, ecology, and agriculture.

Scope

Radiation components. Calculate direct beam, diffuse sky, and terrain-reflected radiation separately. The total is the sum of all three.

Terrain effects. Account for slope, aspect, and horizon shadowing derived from the DEM. Cells behind ridgelines should receive reduced or zero direct beam depending on the sun's position.

Atmospheric transmittance. Configurable transmittance parameter (or a simple atmospheric model) to control how much radiation reaches the surface.

Integration periods. Support hourly snapshots, daily totals, or annual sums. Daily/annual modes integrate over sun positions throughout the period.

Sun position calculation. Compute solar declination and hour angle from date/time and latitude. No external ephemeris library should be required for basic operation.

Prior art

Similar in scope to GRASS r.sun or ArcGIS Area Solar Radiation. Both are widely used but not available as lightweight Python functions that work on xarray DataArrays.

Implementation notes

The horizon calculation is the expensive part. For each cell, you need to trace rays in multiple azimuth directions and find the maximum elevation angle. This could reuse or extend the existing viewshed infrastructure if the interface fits.

GPU acceleration would help a lot here since the horizon sweep is embarrassingly parallel across cells. A CUDA kernel that walks each azimuth ray per thread block is a natural fit.

For dask, the horizon calculation requires overlap (the search radius can be large), so map_overlap with sufficient depth is the right approach.