diff --git a/conf.py b/conf.py
new file mode 100644
index 0000000000000000000000000000000000000000..3ab6030c56f2b10fa33cff076dce390af4a0a03f
--- /dev/null
+++ b/conf.py
@@ -0,0 +1,96 @@
+import numpy as np
+
+
+def conf_test(rad, 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 ---------->
+ '''
+
+ radshape = rad.shape
+ rad = rad.reshape(np.prod(radshape))
+
+ if thr.ndim == 1:
+ 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]
+ 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
+ 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
+
+ # 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))
+ confidence[idx] = coeff[idx] * np.power(s1[idx], power[idx])
+ idx = np.nonzero((rad <= beta) & (flipped is True))
+ 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))
+ confidence[idx] = coeff[idx] * np.power(s1[idx], power[idx])
+ idx = np.nonzero((rad > beta) & (flipped is True))
+ confidence[idx] = coeff[idx] * np.power(s1[idx], power[idx])
+
+ # Rad outside alpha-gamma interval
+ confidence[(rad > gamma) & (flipped is False)] = 1
+ confidence[(rad < alpha) & (flipped is False)] = 0
+ confidence[(rad > gamma) & (flipped is True)] = 0
+ confidence[(rad < alpha) & (flipped is True)] = 1
+
+ confidence[confidence > 1] = 1
+ confidence[confidence < 0] = 0
+
+ return confidence
+
+
+
diff --git a/main.py b/main.py
index 29acf57277e0e8608d5b4e0593fe54316bb4f297..a659e84e730e0d483ff164df68154a4d3afdbb65 100644
--- a/main.py
+++ b/main.py
@@ -11,17 +11,21 @@ def main():
fname_geo = 'VNP03MOD.A2014213.1548.001.2017301015705.uwssec.nc'
thresh_file = '/home/pveglio/mvcm_leo/thresholds/new_thresholds.mvcm.snpp.v1.0.0.yaml'
- viirs_data = rd.read_data(f'{datapath}/{fname_l1b}', f'{datapath}/{fname_geo}')
+ viirs_data = rd.read_data('viirs', f'{datapath}/{fname_l1b}', f'{datapath}/{fname_geo}')
with open(thresh_file) as f:
text = f.read()
thresholds = yml.safe_load(text)
- confidence = np.zeros(2, viirs_data['M01'].shape[0], viirs_data['M01'].shape[1])
+ confidence = np.zeros((3, 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,
- thresholds['Daytime_Ocean'])
+ thresholds['Daytime_Ocean'], viirs_data,
+ thresholds['Sun_Glint']['bounds'][3])
+
+ confidence[2, :, :] = tests.nir_refl_test(viirs_data.M07.values, thresholds['Daytime_Ocean'],
+ thresholds['Sun_Glint'], viirs_data)
return confidence
diff --git a/read_data.py b/read_data.py
index 036847406c3b61b314104e19b377806265d99c5f..06ac0306ddbf96de760d1b4fb4ac1b4ca0069e27 100644
--- a/read_data.py
+++ b/read_data.py
@@ -31,11 +31,13 @@ def read_data(sensor: str, l1b_filename: str, geo_filename: str):
in_data = in_data.merge(data)
- relazi = relative_azimuth_angle(data['sensor_azimuth'].values, data['solar_azimuth'].values)
- sunglint = sun_glint_angle(data['sensor_zenith'].values, data['solar_zenith'].values, relazi)
+ relazi = relative_azimuth_angle(data.sensor_azimuth.values, data.solar_azimuth.values)
+ sunglint = sun_glint_angle(data.sensor_zenith.values, data.solar_zenith.values, relazi)
+ scatt_angle = scattering_angle(data.solar_zenith.values, data.sensor_zenith.values, relazi)
in_data['relative_azimuth'] = (('number_of_lines', 'number_of_pixels'), relazi)
in_data['sunglint_angle'] = (('number_of_lines', 'number_of_pixels'), sunglint)
+ in_data['scattering_angle'] = (('number_of_lines', 'number_of_pixels'), scatt_angle)
return in_data
@@ -48,9 +50,18 @@ def relative_azimuth_angle(sensor_azimuth: float, solar_azimuth: float) -> float
def sun_glint_angle(sensor_zenith: float, solar_zenith: float, rel_azimuth: float) -> float:
cossna = (np.sin(sensor_zenith*_DTR) * np.sin(solar_zenith*_DTR) * np.cos(rel_azimuth*_DTR) +
np.cos(sensor_zenith*_DTR) * np.cos(solar_zenith*_DTR))
+ cossna[cossna > 1] = 1
sunglint_angle = np.arccos(cossna) * _RTD
return sunglint_angle
+def scattering_angle(solar_zenith, sensor_zenith, relative_azimuth):
+ cos_scatt_angle = -1. * (np.cos(solar_zenith*_DTR) * np.cos(sensor_zenith*_DTR) -
+ np.sin(solar_zenith*_DTR) * np.sin(sensor_zenith*_DTR) *
+ np.cos(relative_azimuth*_DTR))
+ scatt_angle = np.arccos(cos_scatt_angle) * _RTD
+ return scatt_angle
+
+
def correct_reflectances():
pass
diff --git a/tests.py b/tests.py
index 5564ea7e334850130b4972ef5172cf65b9f78e3b..aee86cd9f0ea62b107646ac77f913f4bda8eb8fd 100644
--- a/tests.py
+++ b/tests.py
@@ -1,46 +1,67 @@
import numpy as np
import utils
+import conf
def test_11um(rad, threshold):
- thr = threshold['bt11']
+ radshape = rad.shape
+ rad = rad.reshape(np.prod(radshape))
+
+ thr = np.array(threshold['bt11'])
+ confidence = np.zeros(rad.shape)
+
if thr[4] == 1:
# 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)
- return confidence
+def test_11_4diff(rad1, rad2, threshold, viirs_data, sg_thresh):
-def test_11_4diff(rad1, rad2, threshold):
+ radshape = rad1.shape
+ raddiff = (rad1 - rad2).reshape(np.prod(radshape))
- day = True
- sunglint = False
+ day = np.zeros(radshape)
+ day[viirs_data.solar_zenith <= 85] = 1
+ day = day.reshape(raddiff.shape)
+ sunglint = np.zeros(rad1.shape)
+ sunglint[viirs_data.sunglint_angle <= sg_thresh] = 1
+ sunglint = sunglint.reshape(raddiff.shape)
thr = threshold['test11_4lo']
- raddiff = rad1 - rad2
+ confidence = np.zeros(raddiff.shape)
- if day is True and sunglint is False:
- confidence = utils.conf_test(raddiff, thr)
+ confidence[(day == 1) & (sunglint == 0)] = utils.conf_test(raddiff[(day == 1) & (sunglint == 0)], thr)
- return confidence
+ return confidence.reshape(radshape)
-def nir_refl_test(rad, threshold):
+def nir_refl_test(rad, threshold, sunglint_thresholds, viirs_data):
- sunglint = False
- sza = 0
- refang = 0
- vza = 0
- dtr = 0
+ sza = viirs_data.solar_zenith.values
+ refang = viirs_data.sunglint_angle.values
+ vza = viirs_data.sensor_zenith.values
+ dtr = np.pi/180
+ # Keep in mind that band_n uses MODIS band numbers (i.e. 2=0.86um and 7=2.1um)
+ # For VIIRS would be 2=M07 (0.865um) and 7=M11 (2.25um)
band_n = 2
- vzcpow = 0
+ vzcpow = 0.75 # THIS NEEDS TO BE READ FROM THE THRESHOLDS FILE
+
+ radshape = rad.shape
+ rad = rad.reshape(np.prod(radshape))
+ confidence = np.zeros(rad.shape)
+ sza = sza.reshape(rad.shape)
+ vza = vza.reshape(rad.shape)
+ refang = refang.reshape(rad.shape)
+ sunglint_flag = utils.sunglint_scene(refang, sunglint_thresholds).reshape(rad.shape)
- sunglint_thr = np.zeros((4,))
# ref2 [5]
# b2coeffs [4]
# b2mid [1]
@@ -48,35 +69,38 @@ def nir_refl_test(rad, threshold):
# b2lo [1]
# vzcpow [3] (in different place)
+ cosvza = np.cos(vza*dtr)
coeffs = threshold['b2coeffs']
- hicut0 = coeffs[0] + coeffs[1]*sza + coeffs[2]*np.power(sza, 2) + coeffs[3]*np.power(sza, 3)
+ hicut0 = np.array(coeffs[0] + coeffs[1]*sza + coeffs[2]*np.power(sza, 2) + coeffs[3]*np.power(sza, 3))
hicut0 = (hicut0 * 0.01) + threshold['b2adj']
hicut0 = hicut0 * threshold['b2bias_adj']
midpt0 = hicut0 + (threshold['b2mid'] * threshold['b2bias_adj'])
locut0 = midpt0 + (threshold['b2lo'] * threshold['b2bias_adj'])
-
- if sunglint is True:
- # Keep in mind that band_n uses MODIS band numbers (2=0.86um and 7=2.1um)
- # For VIIRS would be 2=M7 (0.865um) and 7=M11 (2.25um)
- sunglint_thr = utils.get_sunglint_thresholds(refang, threshold['Sun_Glint'],
- band_n, sunglint_thr)
- locut1 = sunglint[0]
- midpt1 = sunglint[1]
- hicut1 = sunglint[2]
- else:
- locut1 = locut0
- midpt1 = midpt0
- hicut1 = hicut0
-
- cosvza = np.cos(vza*dtr)
- locut2 = locut1 * (1./np.power(cosvza, vzcpow[0]))
- midpt2 = midpt1 * (1./np.power(cosvza, vzcpow[0]))
- hicut2 = hicut1 * (1./np.power(cosvza, vzcpow[0]))
-
- corr_thr = [locut2, hicut2, 1.0, midpt2]
- confidence = utils.conf_test(rad, corr_thr)
-
- return confidence
+ 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]))
+
+ 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):
diff --git a/utils.py b/utils.py
index 5be3d78aad42764e0c96c6288027983243a5f1c5..044ec89fcccc1121005a4d13ada215d37b58c469 100644
--- a/utils.py
+++ b/utils.py
@@ -13,16 +13,16 @@ _test_thr = [5, 10, 15, 1, 1]
# over all pixels one by one.
def sunglint_scene(refang, sunglint_thr):
sunglint_flag = np.zeros(refang.shape)
- sunglint_flag[refang <= sunglint_thr['bounds'][3]] = 1
- sunglint_flag[refang <= sunglint_thr['bounds'][2]] = 2
+ 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):
+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
@@ -30,12 +30,15 @@ def get_sunglint_thresholds(refang, thresholds, band_n, sunglint):
# # dosgref[1] = mdcnf
# # sunglint[3] = doref2[3]
- if refang <= thresholds['bounds'][1]:
- sunglint = thresholds[f'{band}_0deg']
+ # if refang <= thresholds['bounds'][1]:
+ if sunglint_flag == 3:
+ sunglint_thr = np.full((4, len(refang)), thresholds[f'{band}_0deg'].T)
+ # sunglint_thr[:, :] = thresholds[f'{band}_0deg']
else:
- if (refang > thresholds['bounds'][1] and refang <= thresholds['bounds'][2]):
+ # 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]
@@ -43,22 +46,25 @@ def get_sunglint_thresholds(refang, thresholds, band_n, sunglint):
power = thresholds[f'{band}_10deg'][3]
conf_range = thresholds[f'{band}_10deg'][2]
- elif (refang > thresholds['bounds'][2] and refang <= thresholds['bounds'][3]):
+ # 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 = sunglint[1]
+ 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[1] = midpt
- sunglint[2] = midpt - conf_range
- sunglint[0] = midpt + conf_range
- sunglint[3] = power
+ sunglint_thr[1, :] = midpt
+ sunglint_thr[2, :] = midpt - conf_range
+ sunglint_thr[0, :] = midpt + conf_range
+ sunglint_thr[3, :] = power
- return sunglint
+ return sunglint_thr
def conf_test(rad=_test_rad, thr=_test_thr):
@@ -87,9 +93,7 @@ def conf_test(rad=_test_rad, thr=_test_thr):
beta = thr[1]
locut = thr[2]
power = thr[3]
- radshape = rad.shape
- rad = rad.reshape((rad.shape[0]*rad.shape[1]))
- c = np.zeros(rad.shape)
+ confidence = np.zeros(rad.shape)
if hicut > locut:
gamma = thr[0]
@@ -104,30 +108,28 @@ def conf_test(rad=_test_rad, thr=_test_thr):
range_ = 2. * (beta - alpha)
s1 = (rad[rad <= beta] - alpha) / range_
if flipped is False:
- c[rad <= beta] = coeff * np.power(s1, power)
+ confidence[rad <= beta] = coeff * np.power(s1, power)
if flipped is True:
- c[rad <= beta] = 1. - (coeff * np.power(s1, power))
+ 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:
- c[rad > beta] = 1. - (coeff * np.power(s1, power))
+ confidence[rad > beta] = 1. - (coeff * np.power(s1, power))
if flipped is True:
- c[rad > beta] = coeff * np.power(s1, power)
+ confidence[rad > beta] = coeff * np.power(s1, power)
# Rad outside alpha-gamma interval
if flipped is False:
- c[rad > gamma] = 1
- c[rad < alpha] = 0
+ confidence[rad > gamma] = 1
+ confidence[rad < alpha] = 0
if flipped is True:
- c[rad > gamma] = 0
- c[rad < alpha] = 1
+ confidence[rad > gamma] = 0
+ confidence[rad < alpha] = 1
- c[c > 1] = 1
- c[c < 0] = 0
-
- confidence = c.reshape(radshape)
+ confidence[confidence > 1] = 1
+ confidence[confidence < 0] = 0
return confidence