Skip to content
Snippets Groups Projects
Commit 42d85a90 authored by Paolo Veglio's avatar Paolo Veglio
Browse files

added more tests. reworked conf_test_dble to be faster

parent f3de871b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
conf.py
+ 112
− 30
View file @ 42d85a90
......@@ -29,56 +29,36 @@ def conf_test(rad, thr):
thr = np.full((rad.shape[0], 4), thr[:4]).T
coeff = np.power(2, (thr[3] - 1))
hicut = thr[0]
beta = thr[1]
locut = thr[2]
power = thr[3]
hicut = thr[0, :]
beta = thr[1, :]
locut = thr[2, :]
power = thr[3, :]
confidence = np.zeros(rad.shape)
alpha, gamma = np.empty(rad.shape), np.empty(rad.shape)
flipped = np.zeros(rad.shape)
# if thr.ndim == 2:
gamma[hicut > locut] = thr[0, hicut > locut]
alpha[hicut > locut] = thr[2, hicut > locut]
flipped[hicut > locut] = False
flipped[hicut > locut] = 0
gamma[hicut < locut] = thr[2, hicut < locut]
alpha[hicut < locut] = thr[0, hicut < locut]
flipped[hicut < locut] = True
# else:
# gamma[hicut > locut] = thr[0]
# alpha[hicut > locut] = thr[2]
# flipped[hicut > locut] = False
# gamma[hicut < locut] = thr[2]
# alpha[hicut < locut] = thr[0]
# flipped[hicut < locut] = True
flipped[hicut < locut] = 1
# 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))
range_ = 2. * (beta - alpha)
s1 = (rad - alpha) / range_
idx = np.nonzero((rad <= beta) & (flipped is False))
idx = np.nonzero((rad <= beta) & (flipped == 0))
confidence[idx] = coeff[idx] * np.power(s1[idx], power[idx])
idx = np.nonzero((rad <= beta) & (flipped is True))
idx = np.nonzero((rad <= beta) & (flipped == 1))
confidence[idx] = coeff[idx] * np.power(s1[idx], power[idx])
# 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)
range_ = 2. * (beta - gamma)
s1 = (rad - alpha) / range_
idx = np.nonzero((rad > beta) & (flipped is False))
idx = np.nonzero((rad > beta) & (flipped == 0))
confidence[idx] = coeff[idx] * np.power(s1[idx], power[idx])
idx = np.nonzero((rad > beta) & (flipped is True))
idx = np.nonzero((rad > beta) & (flipped == 1))
confidence[idx] = coeff[idx] * np.power(s1[idx], power[idx])
# Rad outside alpha-gamma interval
......@@ -93,4 +73,106 @@ def conf_test(rad, thr):
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 <= alpha1) & (rad < beta1) & (alpha1 - gamma1 > 0))
confidence[idx] = 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] = coeff[idx] * np.power(s1[idx], power[idx])
# Check for value beyond function range
confidence[(rad > alpha1) & (rad < alpha2)] = 0
confidence[(rad < gamma1) | (rad > gamma2)] = 1
###
# Inner region passes test
# Value is within range of lower set of limits
range_ = 2 * (beta1 - alpha1)
s1 = (rad - alpha1) / range_
idx = np.nonzero((rad <= gamma1) & (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] = 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 > gamma1) & (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 > gamma1) & (rad < beta2) & (alpha1 - gamma1 > 0))
confidence[idx] = coeff[idx] * np.power(s1[idx], power[idx])
confidence[(rad > gamma1) & (rad < gamma2)] = 1
confidence[(rad < alpha1) & (rad > alpha2)] = 0
confidence[confidence > 1] = 1
confidence[confidence < 0] = 0
return confidence
main.py
+ 10
− 1
View file @ 42d85a90
......@@ -17,7 +17,7 @@ def main():
text = f.read()
thresholds = yml.safe_load(text)
confidence = np.zeros((3, viirs_data['M01'].shape[0], viirs_data['M01'].shape[1]))
confidence = np.zeros((5, viirs_data['M01'].shape[0], viirs_data['M01'].shape[1]))
confidence[0, :, :] = tests.test_11um(viirs_data.M15.values, thresholds['Daytime_Ocean'])
confidence[1, :, :] = tests.test_11_4diff(viirs_data.M15.values, viirs_data.M13.values,
......@@ -27,6 +27,15 @@ def main():
confidence[2, :, :] = tests.nir_refl_test(viirs_data.M07.values, thresholds['Daytime_Ocean'],
thresholds['Sun_Glint'], viirs_data)
# Note that here I'm using M05/M07 but the corresponding hi-res channelsa re I1/I2
# IMPORTANT: conf_test_dble() needs to be verified. I don't think it's working as intended at the moment
confidence[3, :, :] = tests.vis_nir_ratio_test(viirs_data.M05.values, viirs_data.M07.values,
thresholds['Daytime_Ocean'], thresholds['Sun_Glint'])
# This test needs to be verified, for the granule I'm running everything is zero
confidence[4, :, :] = tests.test_11um_var(viirs_data.M15.values, thresholds['Nighttime_Ocean'],
thresholds['Daytime_Ocean_Spatial_Variability'])
return confidence
......
tests.py
+ 93
− 19
View file @ 42d85a90
import numpy as np
from numpy.lib.stride_tricks import sliding_window_view
import utils
import conf
......@@ -13,19 +15,44 @@ def test_11um(rad, threshold):
confidence = np.zeros(rad.shape)
if thr[4] == 1:
print("11um test running")
# the C code has the line below that I don't quite understand the purpose of.
# It seems to be setting the bit to 0 if the BT value is greater than the midpoint
#
# if (m31 >= dobt11[1]) (void) set_bit(13, pxout.testbits);
#confidence = utils.conf_test(rad, thr)
# confidence = utils.conf_test(rad, thr)
confidence = conf.conf_test(rad, thr)
return confidence.reshape(radshape)
def test_11um_var(rad, threshold, var_threshold):
print("11um variability test running")
thr = np.array(threshold['11um_var'])
radshape = rad.shape
var = np.zeros((radshape[0], radshape[1], 9))
# chk_spatial2() need to figure out what this is
# np = rg_var.num_small_diffs * 1.0
test = sliding_window_view(np.pad(rad, [1, 1], mode='constant'), (3, 3)) - np.expand_dims(rad, (2, 3))
var[np.abs(test).reshape(radshape[0], radshape[1], 9) < var_threshold['dovar11']] = 1
var = var.sum(axis=2).reshape(np.prod(radshape))
rad = rad.reshape(np.prod(radshape))
confidence = np.zeros(rad.shape)
confidence[var == 9] = conf.conf_test(rad[var == 9], thr)
return confidence.reshape(radshape)
def test_11_4diff(rad1, rad2, threshold, viirs_data, sg_thresh):
print("11um - 4um difference test running")
radshape = rad1.shape
raddiff = (rad1 - rad2).reshape(np.prod(radshape))
......@@ -35,16 +62,48 @@ def test_11_4diff(rad1, rad2, threshold, viirs_data, sg_thresh):
sunglint = np.zeros(rad1.shape)
sunglint[viirs_data.sunglint_angle <= sg_thresh] = 1
sunglint = sunglint.reshape(raddiff.shape)
thr = threshold['test11_4lo']
thr = np.array(threshold['test11_4lo'])
confidence = np.zeros(raddiff.shape)
confidence[(day == 1) & (sunglint == 0)] = utils.conf_test(raddiff[(day == 1) & (sunglint == 0)], thr)
# confidence[(day == 1) & (sunglint == 0)] = utils.conf_test(raddiff[(day == 1) & (sunglint == 0)], thr)
confidence[(day == 1) & (sunglint == 0)] = conf.conf_test(raddiff[(day == 1) & (sunglint == 0)], thr)
return confidence.reshape(radshape)
def vir_refl_test(rad, threshold, viirs_data):
print('Visible reflectance test running')
thr = threshold['vis_refl_test']
radshape = rad.shape()
rad = rad.reshape(np.prod(radshape))
confidence = np.zeros(radshape)
vzcpow = 0.75 # THIS NEEDS TO BE READ FROM THE THRESHOLDS FILE
vza = viirs_data.sensor_zenith.values
dtr = np.pi/180
cosvza = np.cos(vza*dtr)
coeffs = utils.get_b1_thresholds()
# this quantity is the return of get_b1_thresholds() in the C code
# it's defined here to keep a consistent logic with the original source
irtn = 1
if irtn != 0:
coeffs = thr
coeffs[:3] = coeffs[:3] * 1/np.power(cosvza, vzcpow)
confidence = conf.conf_test(rad, coeffs)
return confidence.reshape(radshape)
def nir_refl_test(rad, threshold, sunglint_thresholds, viirs_data):
print("NIR reflectance test running")
sza = viirs_data.solar_zenith.values
refang = viirs_data.sunglint_angle.values
vza = viirs_data.sensor_zenith.values
......@@ -77,34 +136,49 @@ def nir_refl_test(rad, threshold, sunglint_thresholds, viirs_data):
midpt0 = hicut0 + (threshold['b2mid'] * threshold['b2bias_adj'])
locut0 = midpt0 + (threshold['b2lo'] * threshold['b2bias_adj'])
thr = np.array([locut0, midpt0, hicut0, threshold['ref2'][3]*np.ones(rad.shape)])
print(thr.shape)
# corr_thr = np.zeros((4, 4))
corr_thr = np.zeros((4, rad.shape[0]))
corr_thr[:3, sunglint_flag == 0] = thr[:3, sunglint_flag == 0] * (1./np.power(cosvza[sunglint_flag == 0], vzcpow))
corr_thr[3, sunglint_flag == 0] = thr[3, sunglint_flag == 0]
# for flag in range(1, 4):
# sunglint_thr = utils.get_sunglint_thresholds(refang, sunglint_thresholds, band_n, flag, thr)
# corr_thr[:3, sunglint_flag == flag] = sunglint_thr[:3, :] * (1./np.power(cosvza, vzcpow[0]))
# corr_thr[3, sunglint_flag == flag] = sunglint_thr[3, :]
# corr_thr[flag, :3] = sunglint_thr[:3] * (1./np.power(cosvza, vzcpow[0]))
for flag in range(1, 4):
if len(refang[sunglint_flag == flag]) > 0:
sunglint_thr = utils.get_sunglint_thresholds(refang, sunglint_thresholds, band_n, flag, thr)
corr_thr[:3, sunglint_flag == flag] = sunglint_thr[:3, sunglint_flag == flag] * (1./np.power(cosvza[sunglint_flag == flag], vzcpow))
corr_thr[3, sunglint_flag == flag] = sunglint_thr[3, sunglint_flag == flag]
confidence = conf.conf_test(rad, corr_thr)
#confidence[sunglint_flag == 1] = utils.conf_test(rad[sunglint_flag == 1], corr_thr[1, :])
#confidence[sunglint_flag == 2] = utils.conf_test(rad[sunglint_flag == 2], corr_thr[2, :])
#confidence[sunglint_flag == 3] = utils.conf_test(rad[sunglint_flag == 3], corr_thr[3, :])
#confidence[sunglint_flag == 0] = utils.conf_test(rad[sunglint_flag == 0], corr_thr[0, :])
#for flag in range(1, 4):
# sunglint_thr = utils.get_sunglint_thresholds(refang, sunglint_thresholds, band_n, flag, thr)
# pass
return confidence.reshape(radshape)
def vis_nir_ratio_test(rad1, rad2, threshold):
pass
def vis_nir_ratio_test(rad1, rad2, threshold, sg_threshold):
print("NIR-Visible ratio test running")
if threshold['vis_nir_ratio'][6] == 1:
radshape = rad1.shape
rad1 = rad1.reshape(np.prod(radshape))
rad2 = rad2.reshape(np.prod(radshape))
vrat = rad2/rad1
thresh = np.zeros((7,))
# temp value to avoid linter bitching at me
# eventually we would have the test run in two blocks as:
# confidence[sunglint == 1] = conf.conf_test_dble(vrat[sunglint == 1], sg_threshold['snglnt'])
# confidence[sunglint == 0] = conf.conf_test_dble(vrat[sunglint == 0], threshold['vis_nir_ratio'])
# sunglint needs to be defined somewhere
sunglint = 0
if sunglint:
thresh = threshold['snglnt']
else:
thresh = threshold['vis_nir_ratio']
confidence = conf.conf_test_dble(vrat, thresh)
return confidence.reshape(radshape)
class CloudMaskTests(object):
......
utils.py
+ 16
− 1
View file @ 42d85a90
......@@ -32,7 +32,7 @@ def get_sunglint_thresholds(refang, thresholds, band_n, sunglint_flag, thr):
# if refang <= thresholds['bounds'][1]:
if sunglint_flag == 3:
sunglint_thr = np.full((4, len(refang)), thresholds[f'{band}_0deg'].T)
sunglint_thr = np.full((4, len(refang)), thresholds[f'{band}_0deg'])
# sunglint_thr[:, :] = thresholds[f'{band}_0deg']
else:
......@@ -225,3 +225,18 @@ def conf_test_dble(rad, coeffs):
confidence = c.reshape(radshape)
return confidence
def get_b1_threshold(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
if pxin_ndvibk < des_ndvi:
coeffs = thresholds['Coeffs_Band8_land_thresh']
else:
coeffs = thresholds['Coeffs_Band1_land_thresh']
# to be continued ...
return coeffs
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment