From 42d85a9054d55b3a2f31398634c10a5169da65f0 Mon Sep 17 00:00:00 2001
From: Paolo Veglio <paolo.veglio@ssec.wisc.edu>
Date: Tue, 24 May 2022 20:53:51 +0000
Subject: [PATCH] added more tests. reworked conf_test_dble to be faster
---
conf.py | 142 +++++++++++++++++++++++++++++++++++++++++++------------
main.py | 11 ++++-
tests.py | 112 +++++++++++++++++++++++++++++++++++--------
utils.py | 17 ++++++-
4 files changed, 231 insertions(+), 51 deletions(-)
diff --git a/conf.py b/conf.py
index 3ab6030..9582fe1 100644
--- a/conf.py
+++ b/conf.py
@@ -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
diff --git a/main.py b/main.py
index a659e84..84cb602 100644
--- a/main.py
+++ b/main.py
@@ -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
diff --git a/tests.py b/tests.py
index aee86cd..04667c1 100644
--- a/tests.py
+++ b/tests.py
@@ -1,5 +1,7 @@
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):
diff --git a/utils.py b/utils.py
index 044ec89..6734164 100644
--- a/utils.py
+++ b/utils.py
@@ -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
--
GitLab