viewf.py

import numpy as np

from topocalc.gradient import gradient_d8
from topocalc.horizon import horizon


def d2r(a):
    """Angle to radians

    Arguments:
        a {float} -- angle in degrees

    Returns:
        v {float} -- angle in radians
    """
    v = a * np.pi / 180
    v = round(v, 6)  # just for testing at the moment
    return v


def viewf(dem, spacing, nangles=72, sin_slope=None, aspect=None):
    """
    Calculate the sky view factor of a dem.

    The sky view factor from equation 7b from Dozier and Frew 1990

    .. math::
        V_d \approx \frac{1}{2\pi} \int_{0}^{2\pi}\left [ cos(S) sin^2{H_\phi} 
        + sin(S)cos(\phi-A) \times \left ( H_\phi - sin(H_\phi) cos(H_\phi)
        \right )\right ] d\phi

    terrain configuration factor (tvf) is defined as:
        (1 + cos(slope))/2 - sky view factor

    Based on the paper Dozier and Frew, 1990 and modified from
    the Image Processing Workbench code base (Frew, 1990). The
    Python version of sky view factor will be an almost exact
    replication of the IPW command `viewf` minus rounding errors
    from type and linear quantization.

    Args:
        dem: numpy array for the DEM
        spacing: grid spacing of the DEM
        nangles: number of angles to estimate the horizon, defaults
                to 72 angles
        sin_slope: optional, will calculate if not provided
                    sin(slope) with range from 0 to 1
        aspect: optional, will calculate if not provided
                Aspect as radians from south (aspect 0 is toward
                the south) with range from -pi to pi, with negative
                values to the west and positive values to the east.

    Returns:
        svf: sky view factor
        tcf: terrain configuration factor

    """  # noqa

    if dem.ndim != 2:
        raise ValueError('viewf input of dem is not a 2D array')

    if nangles < 16:
        raise ValueError('viewf number of angles should be 16 or greater')

    if sin_slope is not None:
        if np.max(sin_slope) > 1:
            raise ValueError('slope must be sin(slope) with range from 0 to 1')

    # calculate the gradient if not provided
    # The slope is returned as radians so convert to sin(S)
    if sin_slope is None:
        slope, aspect = gradient_d8(
            dem, dx=spacing, dy=spacing, aspect_rad=True)
        sin_slope = np.sin(slope)

    # -180 is North
    angles = np.linspace(-180, 180, num=nangles, endpoint=False)

    # perform the integral
    cos_slope = np.sqrt((1 - sin_slope) * (1 + sin_slope))
    svf = np.zeros_like(sin_slope)
    for angle in angles:

        # horizon angles
        hcos = horizon(angle, dem, spacing)
        azimuth = d2r(angle)

        # sin^2(H)
        sin_squared = (1 - hcos) * (1 + hcos)

        # H - sin(H)cos(H)
        h_mult = np.arccos(hcos) - np.sqrt(sin_squared) * hcos

        # cosines of difference between horizon aspect and slope aspect
        cos_aspect = np.cos(azimuth - aspect)

        # integral in equation 7b
        intgrnd = cos_slope * sin_squared + \
            sin_slope * cos_aspect * h_mult

        ind = intgrnd > 0
        svf[ind] = svf[ind] + intgrnd[ind]

    svf = svf / len(angles)

    tcf = (1 + cos_slope)/2 - svf

    return svf, tcf