Sky view factor calculation from Sirguey et al. 2009

sky_view.py

def skyview(horz_data, slp, asp):
    """ Compute the sky view factor.
    Compute the sky view factor for every pixel, i.e. the ratio between diffuse
    irradiance received to that on unobstructed surface according to its slope,
    aspect and horizon. Algorithm by Dozier et al. 1981, 1990.
    Args:
        horz_data (ndarray): horizon elevation data computed from the DEM
        slp (ndarray): slope calculated from the DEM
        asp (ndarray): aspect calculated from the DEM
    Returns:
        vd (ndarray): skyview factor
        ct (ndarray): terrain configuration factor (counterpart of vt)
    """

    Xsize = horz_data.shape[1]  # get the horizon raster size
    Ysize = horz_data.shape[2]
    N = horz_data.shape[0]  # get the number of viewing azimuths from horizon

    # Preallocate variables
    dphi = 2 * np.pi / N
    cosS = np.cos(np.deg2rad(slp))
    sinS = np.sin(np.deg2rad(slp))

    horz_data = np.deg2rad(90 - horz_data)
    cosH = np.cos(horz_data)
    sinH = np.sin(horz_data)

    # Compute vd (skyview factor)
    vd = np.zeros((Xsize, Ysize))  # Initialise vd

    phi = np.degrees(
        np.multiply.outer(
            np.pi - np.array(range(N)) * dphi, np.ones((Xsize, Ysize))
        )
    )

    cos_diff_phi_az = np.cos(np.deg2rad(phi - asp))

    prod = cos_diff_phi_az * (horz_data - sinH * cosH)
    prod[np.isnan(prod)] = 0

    vd_tmp = (dphi / (2 * np.pi)) * (cosS * sinH ** 2 + sinS * prod)
    vd = np.sum(vd_tmp, axis=0)

    ct = 1 - vd

    return vd, ct