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()