Skip to content
Snippets Groups Projects
conf_xr.py 6.54 KiB
Newer Older
import numpy as np
import xarray as xr


def test(flipped=False):
    bt = np.arange(265, 275)
    if flipped is False:
        thr = np.array([267, 270, 273, 1, 1])
    else:
        thr = np.array([273, 270, 267, 1, 1])
    c = conf_test(bt, thr)
    print(c)


def test_dble(flipped=False):
    bt = np.arange(260, 282)
    if flipped is False:
        thr = np.array([264, 267, 270, 273, 276, 279, 1, 1])
    else:
        thr = np.array([279, 276, 273, 270, 267, 264, 1, 1])
    c = conf_test_dble(bt, thr)
    print(c)


def conf_test(data, band):
    '''
    Assuming a linear function between min and max confidence level, the plot below shows
    how the confidence (y axis) is computed as function of radiance (x axis).
    This case illustrates alpha < gamma, obviously in case alpha > gamma, the plot would be
    flipped.
                       gamma
    c  1                 ________
    o  |                /
    n  |               /
    f  |              /
    i  |     beta    /
    d 1/2    |....../
    e  |           /
    n  |          /
    c  |         /
    e  0________/
       |      alpha
    --------- radiance ---------->
    '''

    hicut = data.threshold[:, :, 2]
    beta = data.threshold[:, :, 1]
    locut = data.threshold[:, :, 0]
    power = data.threshold[:, :, 3]
    coeff = np.power(2, (power - 1))

    gamma = data.threshold.where(hicut > locut, data.threshold[:, :, 0])[:, :, 2]
    alpha = data.threshold.where(hicut > locut, data.threshold[:, :, 2])[:, :, 0]
    flipped = xr.zeros_like(data[band]).where(hicut > locut, 1)

    # Rad between alpha and beta
    range_ = 2. * (beta - alpha)
    s1 = (data[band].values - alpha)/range_
    conf_tmp1 = (coeff * np.power(s1, power)).where((data[band] <= beta) & (flipped == 0))
    conf_tmp2 = (1.0 - coeff * np.power(s1, power)).where((data[band] <= beta) & (flipped == 1))
    conf_tmp12 = conf_tmp1.where(flipped == 0, conf_tmp2)

    # Rad between beta and gamma
    range_ = 2. * (beta - gamma)
    s1 = (data[band].values - gamma)/range_
    conf_tmp3 = (1.0 - coeff * np.power(s1, power)).where((data[band] <= beta) & (flipped == 0))
    conf_tmp4 = (coeff * np.power(s1, power)).where((data[band] <= beta) & (flipped == 1))
    conf_tmp34 = conf_tmp3.where(flipped == 0, conf_tmp4)

    confidence = conf_tmp12.where(data[band] <= beta, conf_tmp34)

    confidence = confidence.where(confidence > 0, 0)
    confidence = confidence.where(confidence < 1, 1)

    return confidence


def conf_test_dble(rad, coeffs):
    # '''
    #            gamma1                         gamma2
    #    c  1_______                               ________
    #    o  |       \                             /
    #    n  |        \                           /
    #    f  |         \                         /
    #    i  |          \   beta1       beta2   /
    #    d 1/2          \....|          |...../
    #    e  |            \                   /
    #    n  |             \                 /
    #    c  |              \               /
    #    e  0               \_____________/
    #       |             alpha1       alpha2
    #    --------------------- radiance ------------------------->
    # '''

    coeffs = np.array(coeffs)
    radshape = rad.shape
    rad = rad.reshape(np.prod(radshape))
    confidence = np.zeros(rad.shape)

    alpha1, gamma1 = np.empty(rad.shape), np.empty(rad.shape)
    alpha2, gamma2 = np.empty(rad.shape), np.empty(rad.shape)

    if coeffs.ndim == 1:
        coeffs = np.full((rad.shape[0], 7), coeffs[:7]).T

    gamma1 = coeffs[0, :]
    beta1 = coeffs[1, :]
    alpha1 = coeffs[2, :]
    alpha2 = coeffs[3, :]
    beta2 = coeffs[4, :]
    gamma2 = coeffs[5, :]
    power = coeffs[6, :]

    coeff = np.power(2, (power - 1))
    # radshape = rad.shape
    # rad = rad.reshape((rad.shape[0]*rad.shape[1]))

    # ## Find if interval between inner cutoffs passes or fails test

    # Inner region fails test

    # Value is within range of lower set of limits
    range_ = 2. * (beta1 - alpha1)
    s1 = (rad - alpha1) / range_
    idx = np.nonzero((rad <= alpha1) & (rad >= beta1) & (alpha1 - gamma1 > 0))
    confidence[idx] = coeff[idx] * np.power(s1[idx], power[idx])

    range_ = 2. * (beta1 - gamma1)
    s1 = (rad - gamma1) / range_
    idx = np.nonzero((rad >= gamma1) & (rad < beta1) & (alpha1 - gamma1 > 0))
    confidence[idx] = 1.0 - coeff[idx] * np.power(s1[idx], power[idx])

    # Value is within range of upper set of limits
    range_ = 2. * (beta2 - alpha2)
    s1 = (rad - alpha2) / range_
    idx = np.nonzero((rad > alpha1) & (rad <= beta2) & (alpha1 - gamma1 > 0))
    confidence[idx] = coeff[idx] * np.power(s1[idx], power[idx])

    range_ = 2. * (beta2 - gamma2)
    s1 = (rad - gamma2) / range_
    idx = np.nonzero((rad > alpha1) & (rad > beta2) & (alpha1 - gamma1 > 0))
    confidence[idx] = 1.0 - coeff[idx] * np.power(s1[idx], power[idx])

    # Check for value beyond function range
    confidence[(alpha1 - gamma1 > 0) & (rad > alpha1) & (rad < alpha2)] = 0
    confidence[(alpha1 - gamma1 > 0) & ((rad < gamma1) | (rad > gamma2))] = 1

    ###

    # Inner region passes test
    print("I NEED TO REVIEW THIS TO WRITE IT MORE CLEARLY")
    # FOR NOW ALPHA AND GAMMA ARE SWITCHED BECAUSE OF HOW THE ARRAYS ARE DEFINED.
    # THINK ON HOW THIS COULD BE WRITTEN SO THAT IT'S EASIER TO UNDERSTAND (AND DEBUG)
    # Value is within range of lower set of limits
    range_ = 2 * (beta1 - alpha1)
    s1 = (rad - alpha1) / range_
    idx = np.nonzero((rad > alpha1) & (rad <= gamma1) & (rad <= beta1) & (alpha1 - gamma1 <= 0))
    confidence[idx] = 1.0 - coeff[idx] * np.power(s1[idx], power[idx])

    range_ = 2 * (beta1 - gamma1)
    s1 = (rad - gamma1) / range_
    idx = np.nonzero((rad > alpha1) & (rad <= gamma1) & (rad > beta1) & (alpha1 - gamma1 <= 0))
    confidence[idx] = coeff[idx] * np.power(s1[idx], power[idx])

    # Values is within range of upper set of limits
    range_ = 2 * (beta2 - alpha2)
    s1 = (rad - alpha2) / range_
    idx = np.nonzero((rad > gamma2) & (rad < alpha2) & (rad >= beta2) & (alpha1 - gamma1 <= 0))
    confidence[idx] = 1.0 - coeff[idx] * np.power(s1[idx], power[idx])

    range_ = 2 * (beta2 - gamma2)
    s1 = (rad - gamma2) / range_
    idx = np.nonzero((rad > gamma2) & (rad < alpha2) & (rad < beta2) & (alpha1 - gamma1 <= 0))
    confidence[idx] = coeff[idx] * np.power(s1[idx], power[idx])

    confidence[(alpha1 - gamma1 <= 0) & ((rad > gamma1) | (rad < gamma2))] = 0
    confidence[(alpha1 - gamma1 <= 0) & (rad <= alpha1) & (rad >= alpha2)] = 1

    confidence[confidence > 1] = 1
    confidence[confidence < 0] = 0

    return confidence


if __name__ == "__main__":
    test_dble()