import numpy as np
_test_rad = np.random.randint(25, size=[6, 8])
# _test_thr = [15, 10, 5, 1, 1]
_test_thr = [5, 10, 15, 1, 1]
# this function creates a map of sunglint areas, based on the different angles set in the
# threshold file. The goal is to create an array of indices that I can use to quickly assign
# different coefficients depending on the angle interval. This will be mostly used in the
# function get_sunglint_thresholds().
# All of this is because we want to be able to process the whole array, instead of iterating
# over all pixels one by one.
def sunglint_scene(refang, sunglint_thr):
sunglint_flag = np.zeros(refang.shape)
sunglint_flag[(refang > sunglint_thr['bounds'][2]) & (refang <= sunglint_thr['bounds'][3])] = 1
sunglint_flag[(refang > sunglint_thr['bounds'][1]) & (refang <= sunglint_thr['bounds'][2])] = 2
sunglint_flag[refang <= sunglint_thr['bounds'][1]] = 3
return sunglint_flag
def get_sunglint_thresholds(refang, thresholds, band_n, sunglint_flag, thr):
band = f'band{band_n}'
sunglint_thr = np.zeros((4, len(refang)))
# if refang > thresholds['bounds'][3]:
# sunglint = sunglint
# # dosgref[2] = hicnf
# # dosgref[0] = locnf
# # dosgref[1] = mdcnf
# # sunglint[3] = doref2[3]
# if refang <= thresholds['bounds'][1]:
if sunglint_flag == 3:
sunglint_thr = np.full((4, len(refang)), thresholds[f'{band}_0deg'])
# sunglint_thr[:, :] = thresholds[f'{band}_0deg']
else:
# if (refang > thresholds['bounds'][1] and refang <= thresholds['bounds'][2]):
if sunglint_flag == 2:
lo_ang = thresholds['bounds'][1]
hi_ang = thresholds['bounds'][2]
lo_ang_val = thresholds[f'{band}_10deg'][0]
hi_ang_val = thresholds[f'{band}_10deg'][1]
power = thresholds[f'{band}_10deg'][3]
conf_range = thresholds[f'{band}_10deg'][2]
# elif (refang > thresholds['bounds'][2] and refang <= thresholds['bounds'][3]):
elif sunglint_flag == 1:
lo_ang = thresholds['bounds'][2]
hi_ang = thresholds['bounds'][3]
lo_ang_val = thresholds[f'{band}_20deg'][0]
hi_ang_val = thr[1]
power = thresholds[f'{band}_20deg'][3]
conf_range = thresholds[f'{band}_20deg'][2]
else:
raise ValueError("Wrong sunglint flag")
a = (refang - lo_ang) / (hi_ang - lo_ang)
midpt = lo_ang_val + a*(hi_ang_val - lo_ang_val)
sunglint_thr[1, :] = midpt
sunglint_thr[2, :] = midpt - conf_range
sunglint_thr[0, :] = midpt + conf_range
sunglint_thr[3, :] = power
return sunglint_thr
def conf_test(rad=_test_rad, thr=_test_thr):
'''
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 ---------->
'''
coeff = np.power(2, (thr[3] - 1))
hicut = thr[0]
beta = thr[1]
locut = thr[2]
power = thr[3]
confidence = np.zeros(rad.shape)
if hicut > locut:
gamma = thr[0]
alpha = thr[2]
flipped = False
else:
gamma = thr[2]
alpha = thr[0]
flipped = True
# Rad between alpha and beta
range_ = 2. * (beta - alpha)
s1 = (rad[rad <= beta] - alpha) / range_
if flipped is False:
confidence[rad <= beta] = coeff * np.power(s1, power)
if flipped is True:
confidence[rad <= beta] = 1. - (coeff * np.power(s1, power))
# Rad between beta and gamma
range_ = 2. * (beta - gamma)
s1 = (rad[rad > beta] - gamma) / range_
if flipped is False:
confidence[rad > beta] = 1. - (coeff * np.power(s1, power))
if flipped is True:
confidence[rad > beta] = coeff * np.power(s1, power)
# Rad outside alpha-gamma interval
if flipped is False:
confidence[rad > gamma] = 1
confidence[rad < alpha] = 0
if flipped is True:
confidence[rad > gamma] = 0
confidence[rad < alpha] = 1
confidence[confidence > 1] = 1
confidence[confidence < 0] = 0
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 ------------------------->
# '''
hicut = [coeffs[0], coeffs[1]]
locut = [coeffs[2], coeffs[3]]
midpt = [coeffs[4], coeffs[5]]
power = coeffs[6]
gamma1 = hicut[0]
gamma2 = hicut[1]
alpha1 = locut[0]
alpha2 = locut[1]
beta1 = midpt[0]
beta2 = midpt[1]
coeff = np.power(2, (power - 1))
radshape = rad.shape
rad = rad.reshape((rad.shape[0]*rad.shape[1]))
c = np.zeros(rad.shape)
# Find if interval between inner cutoffs passes or fails test
if (alpha1 - gamma1 > 0):
# Value is within range of lower set of limits
range_ = 2 * (beta1 - alpha1)
s1 = (rad[(rad <= alpha1) & (rad >= beta1)] - alpha1) / range_
c[(rad <= alpha1) & (rad >= beta1)] = coeff * np.power(s1, power)
range_ = 2 * (beta1 - gamma1)
s1 = (rad[(rad <= alpha1) & (rad < beta1)] - gamma1) / range_
c[(rad <= alpha1) & (rad < beta1)] = coeff * np.power(s1, power)
# Value is within range of upper set of limits
range_ = 2 * (beta2 - alpha2)
s1 = (rad[(rad > alpha1) & (rad <= beta2)] - alpha2) / range_
c[(rad > alpha1) & (rad <= beta2)] = coeff * np.power(s1, power)
range_ = 2 * (beta2 - gamma2)
s1 = (rad[(rad > alpha1) & (rad > beta2)] - gamma2) / range_
c[(rad > alpha1) & (rad > beta2)] = coeff * np.power(s1, power)
# Inner region fails test
# Check for value beyond function range
c[(rad > alpha1) & (rad < alpha2)] = 0
c[(rad < gamma1) | (rad > gamma2)] = 1
else:
# Value is withing range of lower set of limits
range_ = 2 * (beta1 - alpha1)
s1 = (rad[(rad <= gamma1) & (rad <= beta1)] - alpha1) / range_
c[(rad <= gamma1) & (rad <= beta1)] = coeff * np.power(s1, power)
range_ = 2 * (beta1 - gamma1)
s1 = (rad[(rad <= gamma1) & (rad > beta1)] - gamma1) / range_
c[(rad <= gamma1) & (rad > beta1)] = coeff * np.power(s1, power)
# Value is within range of upper set of limits
range_ = 2 * (beta2 - alpha2)
s1 = (rad[(rad > gamma1) & (rad >= beta2)] - alpha2) / range_
c[(rad > gamma1) & (rad >= beta2)] = coeff * np.power(s1, power)
range_ = 2 * (beta2 - gamma2)
s1 = (rad[(rad > gamma1) & (rad < beta2)] - gamma2) / range_
c[(rad > gamma1) & (rad < beta2)] = coeff * np.power(s1, power)
# Inner region passes test
# Check for value beyond function range
c[(rad > gamma1) & (rad < gamma2)] = 1
c[(rad < alpha1) | (rad > alpha2)] = 0
c[c > 1] = 1
c[c < 0] = 0
confidence = c.reshape(radshape)
return confidence
def get_b1_threshold(thresholds, ndvi_thresholds):
# this needs to be read from the thresholds file
des_ndvi = 0.25
pxin_ndvibk = 0 # placeholder until I figure out how to pass this value
pxin_sctang = 0 # placeholder until I figure out how to pass this value
delta_ndvi_bin = 0.1 # defined as is in the original C code
thr = np.zeros((4, len(pxin_ndvibk)))
thr_adj = np.zeros((4, len(pxin_ndvibk)))
coeffs = np.full((10, 3, 4, len(pxin_ndvibk)))
indvi = np.zeros((len(pxin_ndvibk), ))
x = np.zeros((len(pxin_ndvibk), ))
y1 = np.zeros((len(pxin_ndvibk), ))
y2 = np.zeros((len(pxin_ndvibk), ))
interp = np.zeros((len(pxin_ndvibk), ))
# Check for special cases and fill values ('ndvibk'=32.767)
# idx = np.nonzero((pxin_ndvibk < ndvi_thresholds['fill_ndvi'][0]) &
# (pxin_ndvibk > ndvi_thresholds['fill_ndvi'][1]) &
# (pxin_ndvibk < ndvi_thresholds['ndvi_bnd1']))
idx = np.nonzero(pxin_ndvibk < ndvi_thresholds['ndvi_bnd1'])
# all gets set to zero, so no need to do anything for this idx
# idx = np.nonzero((pxin_ndvibk < ndvi_thresholds['fill_ndvi'][0]) &
# (pxin_ndvibk > ndvi_thresholds['fill_ndvi'][1]) &
# (pxin_ndvibk >= ndvi_thresholds['ndvi_bnd2']))
idx = np.nonzero(pxin_ndvibk >= ndvi_thresholds['ndvi_bnd2'])
indvi[idx] = 9
# else
# idx = np.nonzero((pxin_ndvibk < ndvi_thresholds['fill_ndvi'][0]) &
# (pxin_ndvibk > ndvi_thresholds['fill_ndvi'][1]) &
# (pxin_ndvibk >= ndvi_thresholds['ndvi_bnd1']) &
# (pxin_ndvibk < ndvi_thresholds['ndvi_bnd2']))
idx = np.nonzero((pxin_ndvibk >= ndvi_thresholds['ndvi_bnd1']) &
(pxin_ndvibk < ndvi_thresholds['ndvi_bnd2']))
interp[idx] = 1 # this is not really needed anymore, maybe
indvi = (pxin_ndvibk[idx] / delta_ndvi_bin) - 0.5
indvi[indvi < 0] = 0
x = (pxin_ndvibk - delta_ndvi_bin*indvi + delta_ndvi_bin/2) / delta_ndvi_bin
x[x < 0] = 0
x[x > 1] = 1
for i in range(3):
# idx = np.nonzero((pxin_ndvibk < ndvi_thresholds['fill_ndvi'][0]) &
# (pxin_ndvibk > ndvi_thresholds['fill_ndvi'][1]))
y1[:, i] = (coeffs[indvi, i, 0] +
coeffs[indvi, i, 1]*pxin_sctang +
coeffs[indvi, i, 2]*pxin_sctang**2 +
coeffs[indvi, i, 3]*pxin_sctang**3)
# idx = np.nonzero((pxin_ndvibk < ndvi_thresholds['fill_ndvi'][0]) &
# (pxin_ndvibk > ndvi_thresholds['fill_ndvi'][1]) & interp == 1)
y2[idx, i] = (coeffs[indvi[idx]+1, i, 0] +
coeffs[indvi[idx]+1, i, 1]*pxin_sctang[idx] +
coeffs[indvi[idx]+1, i, 2]*pxin_sctang[idx]**2 +
coeffs[indvi[idx]+1, i, 3]*pxin_sctang[idx]**3)
thr[:, i] = ((1.0 - x) * y1) + (x * y2)
thr_adj[pxin_ndvibk < des_ndvi, i] = thr[pxin_ndvibk < des_ndvi, 0]*ndvi_thresholds['adj_fac_desert']
thr_adj[pxin_ndvibk >= des_ndvi, i] = thr[pxin_ndvibk >= des_ndvi, 0]*ndvi_thresholds['adj_fac_land']
thr[:, 3] = 2.0
return thr