solar

solar

Compute Swift (1976) clear-sky potential solar radiation on sloped surfaces.

Port pywatershed's PRMSSolarGeometry.compute_soltab() into a standalone pure-function module for use in hydro-param's derivation step 9 (soltab). The algorithm computes potential clear-sky direct-beam solar radiation on sloped and horizontal surfaces for every day of the year (1--366) using the equivalent-slope method of Lee (1963) with Swift's (1976) extensions.

The module is vectorised over HRUs: all public functions accept arrays of shape (nhru,) for slope, aspect, and latitude, and return arrays of shape (366, nhru) -- one row per day of the year.

Units

• Radiation: Langleys (cal/cm\ :sup:2/day)
• Sunlight duration: hours
• Slope: decimal fraction (rise/run)
• Aspect: degrees clockwise from north
• Latitude: decimal degrees

Notes

This is a faithful port of the Fortran-origin algorithm preserved in pywatershed's Python codebase. Array aliasing behaviour (e.g., t3 = t7 sharing the same ndarray) is intentionally replicated to match pywatershed's results exactly.

References

Swift, L.W., 1976, Algorithm for solar radiation on mountain slopes: Water Resources Research, v. 12, no. 1, p. 108--112. Lee, R., 1963, Evaluation of solar beam irradiation as a climatic parameter of mountain watersheds: Colorado State University Hydrology Papers, no. 2.

See Also

hydro_param.derivations.pywatershed : Derivation step 9 calls compute_soltab.

NDOY module-attribute

NDOY: int = 366

Number of days used in the solar table (includes leap day).

solar_declination module-attribute

solar_declination: NDArray[floating] = (
    0.006918
    - 0.399912 * cos(_yy)
    + 0.070257 * sin(_yy)
    - 0.006758 * cos(2 * _yy)
    + 0.000907 * sin(2 * _yy)
    - 0.002697 * cos(3 * _yy)
    + 0.00148 * sin(3 * _yy)
)

Solar declination angle (radians) for each day of the year, shape (366,).

r1 module-attribute

r1: NDArray[floating] = 60.0 * 2.0 / _obliquity ** 2

Extraterrestrial irradiance corrected for Earth-Sun distance (Langleys/hr), shape (366,).

compute_soltab

compute_soltab(
    slopes: NDArray[floating],
    aspects: NDArray[floating],
    lats: NDArray[floating],
) -> tuple[
    NDArray[np.floating],
    NDArray[np.floating],
    NDArray[np.floating],
]

Compute potential clear-sky solar radiation tables for PRMS.

Implement the Swift (1976) algorithm for potential direct-beam solar radiation on sloped and horizontal surfaces for every day of the year (1--366). This is the public entry point for derivation step 9 (soltab) in the pywatershed parameterization pipeline.

The function runs the core algorithm twice: once for horizontal surfaces (zero slope, zero aspect) and once for the actual sloped surfaces. Both results are needed by PRMS -- the horizontal radiation is used as a reference for computing the slope correction factor.

PARAMETER	DESCRIPTION
slopes	Slope as decimal fraction (rise/run), non-negative. Values above 10.0 trigger a warning since they likely indicate degree input rather than fractional. TYPE: (ndarray, shape(nhru))
aspects	Aspect in degrees, 0--360 clockwise from north (0 = north, 180 = south). Flat HRUs conventionally use aspect = 0. TYPE: (ndarray, shape(nhru))
lats	Latitude in decimal degrees, range [-90, 90]. TYPE: (ndarray, shape(nhru))

RETURNS	DESCRIPTION
soltab_potsw	Potential solar radiation on the sloped surface (Langleys, cal/cm2/day). TYPE: (ndarray, shape(366, nhru))
soltab_horad_potsw	Potential solar radiation on a horizontal surface (Langleys, cal/cm2/day). TYPE: (ndarray, shape(366, nhru))
soltab_sunhrs	Hours of direct sunlight on the sloped surface. TYPE: (ndarray, shape(366, nhru))

RAISES	DESCRIPTION
ValueError	If input arrays have mismatched lengths, are empty, contain NaN, or have out-of-range values (latitude outside [-90, 90] or negative slopes).

Notes

Input validation is strict: NaN values, mismatched array lengths, and out-of-range latitudes all raise ValueError rather than producing silently incorrect output. Slopes > 10 produce a warning but do not error, in case extreme terrain is intentional.

References

Swift, L.W., 1976, Algorithm for solar radiation on mountain slopes: Water Resources Research, v. 12, no. 1, p. 108--112.

See Also

hydro_param.derivations.pywatershed : Step 9 calls this function.

Source code in src/hydro_param/solar.py

def compute_soltab(
    slopes: NDArray[np.floating],
    aspects: NDArray[np.floating],
    lats: NDArray[np.floating],
) -> tuple[NDArray[np.floating], NDArray[np.floating], NDArray[np.floating]]:
    """Compute potential clear-sky solar radiation tables for PRMS.

    Implement the Swift (1976) algorithm for potential direct-beam solar
    radiation on sloped and horizontal surfaces for every day of the
    year (1--366).  This is the public entry point for derivation
    step 9 (soltab) in the pywatershed parameterization pipeline.

    The function runs the core algorithm twice: once for horizontal
    surfaces (zero slope, zero aspect) and once for the actual sloped
    surfaces.  Both results are needed by PRMS -- the horizontal
    radiation is used as a reference for computing the slope correction
    factor.

    Parameters
    ----------
    slopes : ndarray, shape (nhru,)
        Slope as decimal fraction (rise/run), non-negative.  Values
        above 10.0 trigger a warning since they likely indicate
        degree input rather than fractional.
    aspects : ndarray, shape (nhru,)
        Aspect in degrees, 0--360 clockwise from north (0 = north,
        180 = south).  Flat HRUs conventionally use aspect = 0.
    lats : ndarray, shape (nhru,)
        Latitude in decimal degrees, range [-90, 90].

    Returns
    -------
    soltab_potsw : ndarray, shape (366, nhru)
        Potential solar radiation on the sloped surface
        (Langleys, cal/cm2/day).
    soltab_horad_potsw : ndarray, shape (366, nhru)
        Potential solar radiation on a horizontal surface
        (Langleys, cal/cm2/day).
    soltab_sunhrs : ndarray, shape (366, nhru)
        Hours of direct sunlight on the sloped surface.

    Raises
    ------
    ValueError
        If input arrays have mismatched lengths, are empty, contain NaN,
        or have out-of-range values (latitude outside [-90, 90] or
        negative slopes).

    Notes
    -----
    Input validation is strict: NaN values, mismatched array lengths,
    and out-of-range latitudes all raise ``ValueError`` rather than
    producing silently incorrect output.  Slopes > 10 produce a warning
    but do not error, in case extreme terrain is intentional.

    References
    ----------
    Swift, L.W., 1976, Algorithm for solar radiation on mountain slopes:
        Water Resources Research, v. 12, no. 1, p. 108--112.

    See Also
    --------
    hydro_param.derivations.pywatershed : Step 9 calls this function.
    """
    # --- Input validation ---
    if len(slopes) != len(aspects) or len(slopes) != len(lats):
        raise ValueError(
            f"Input arrays must have equal length: "
            f"slopes={len(slopes)}, aspects={len(aspects)}, lats={len(lats)}"
        )
    if len(slopes) == 0:
        raise ValueError("Input arrays must not be empty")
    if np.any(np.isnan(slopes)) or np.any(np.isnan(aspects)) or np.any(np.isnan(lats)):
        nan_counts = (
            np.count_nonzero(np.isnan(slopes)),
            np.count_nonzero(np.isnan(aspects)),
            np.count_nonzero(np.isnan(lats)),
        )
        raise ValueError(
            f"NaN values in input arrays: slopes={nan_counts[0]}, "
            f"aspects={nan_counts[1]}, lats={nan_counts[2]}"
        )
    if np.any(np.abs(lats) > 90):
        raise ValueError(
            f"Latitude values must be in [-90, 90], got range [{lats.min():.2f}, {lats.max():.2f}]"
        )
    if np.any(slopes < 0):
        raise ValueError("Slope values must be non-negative (decimal fraction rise/run)")
    if np.any(slopes > 10):
        logger.warning(
            "Slopes > 10 detected (max=%.1f). Slopes must be decimal fraction "
            "(rise/run), not degrees. Verify input units.",
            float(slopes.max()),
        )

    nhru = len(slopes)

    # Horizontal surface (zero slope, zero aspect)
    horad_potsw, _horad_sunh = _compute_soltab_core(np.zeros(nhru), np.zeros(nhru), lats)

    # Sloped surface
    sloped_potsw, sloped_sunh = _compute_soltab_core(slopes, aspects, lats)

    return sloped_potsw, horad_potsw, sloped_sunh