diff --git a/main.py b/main.py
index 2836c5971f522439a795a51028b078582f7aeac6..29acf57277e0e8608d5b4e0593fe54316bb4f297 100644
--- a/main.py
+++ b/main.py
@@ -1,9 +1,32 @@
-import tests
import ruamel_yaml as yml
import numpy as np
+import read_data as rd
+import tests
+
def main():
+ datapath = '/ships19/hercules/pveglio/neige_data/snpp_test_input'
+ fname_l1b = 'VNP02MOD.A2014213.1548.001.2017301015346.uwssec.bowtie_restored_scaled.nc'
+ 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}')
+
+ 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[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'])
+
+ return confidence
+
+
+def test_main():
rad1 = [[255, 260, 265, 248, 223],
[278, 285, 270, 268, 256],
@@ -13,7 +36,6 @@ def main():
[280, 281, 272, 270, 267]]
thresh_file = '/home/pveglio/mvcm_leo/thresholds/new_thresholds.mvcm.snpp.v1.0.0.yaml'
-
with open(thresh_file) as f:
text = f.read()
@@ -23,12 +45,9 @@ def main():
confidence = np.zeros((2, rad1.shape[0], rad1.shape[1]))
thresholds = yml.safe_load(text)
- thr_11um = thresholds['Daytime_Ocean']
- thr_11_4diff = thresholds['Daytime_Ocean']
-
- confidence[0, :, :] = tests.test_11um(rad1, thr_11um['bt11'])
- confidence[1, :, :] = tests.test_11_4diff(rad1, rad2, thr_11_4diff['test11_4lo'])
+ confidence[0, :, :] = tests.test_11um(rad1, thresholds['Daytime_Ocean'])
+ confidence[1, :, :] = tests.test_11_4diff(rad1, rad2, thresholds['Daytime_Ocean'])
print(f'Confidence[0,:,:]: \n {confidence[0, :, :]}')
print(f'Confidence[1,:,:]: \n {confidence[1, :, :]}')
@@ -37,4 +56,4 @@ def main():
if __name__ == "__main__":
- main()
+ test_main()
diff --git a/read_data.py b/read_data.py
index 90037500d114439738d51760654cb6f9b1a62bad..036847406c3b61b314104e19b377806265d99c5f 100644
--- a/read_data.py
+++ b/read_data.py
@@ -1,8 +1,11 @@
-from netCDF4 import Dataset
+# from netCDF4 import Dataset
import xarray as xr
import numpy as np
+_DTR = np.pi/180.
+_RTD = 180./np.pi
+
def read_data(sensor: str, l1b_filename: str, geo_filename: str):
@@ -28,6 +31,26 @@ 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)
+
+ in_data['relative_azimuth'] = (('number_of_lines', 'number_of_pixels'), relazi)
+ in_data['sunglint_angle'] = (('number_of_lines', 'number_of_pixels'), sunglint)
+
return in_data
+def relative_azimuth_angle(sensor_azimuth: float, solar_azimuth: float) -> float:
+ rel_azimuth = np.abs(180. - np.abs(sensor_azimuth - solar_azimuth))
+ return rel_azimuth
+
+
+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))
+ sunglint_angle = np.arccos(cossna) * _RTD
+ return sunglint_angle
+
+
+def correct_reflectances():
+ pass
diff --git a/tests.py b/tests.py
index a5070ee1d7adbd3a8813235be3c6e606a1beea82..5564ea7e334850130b4972ef5172cf65b9f78e3b 100644
--- a/tests.py
+++ b/tests.py
@@ -1,29 +1,86 @@
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]
+import utils
def test_11um(rad, threshold):
- #if (~np.isnan(rad) or threshold[4] == 1.0):
- confidence = conf_test(rad, threshold)
+ thr = threshold['bt11']
+ 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)
return confidence
def test_11_4diff(rad1, rad2, threshold):
- raddiff = rad1 - rad2;
- confidence = conf_test(raddiff, threshold)
+ day = True
+ sunglint = False
+ thr = threshold['test11_4lo']
+ raddiff = rad1 - rad2
+
+ if day is True and sunglint is False:
+ confidence = utils.conf_test(raddiff, thr)
+
+ return confidence
+
+
+def nir_refl_test(rad, threshold):
+
+ sunglint = False
+ sza = 0
+ refang = 0
+ vza = 0
+ dtr = 0
+ band_n = 2
+ vzcpow = 0
+
+ sunglint_thr = np.zeros((4,))
+ # ref2 [5]
+ # b2coeffs [4]
+ # b2mid [1]
+ # b2bias_adj [1]
+ # b2lo [1]
+ # vzcpow [3] (in different place)
+
+ coeffs = threshold['b2coeffs']
+ hicut0 = 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
-def simple_threshold_test(rad, threshold):
- return conf_test(rad, threshold)
+def vis_nir_ratio_test(rad1, rad2, threshold):
+ pass
class CloudMaskTests(object):
@@ -51,177 +108,12 @@ class CloudMaskTests(object):
pass
-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]
- radshape = rad.shape
- rad = rad.reshape((rad.shape[0]*rad.shape[1]))
- c = 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:
- c[rad <= beta] = coeff * np.power(s1, power)
- if flipped is True:
- c[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))
- if flipped is True:
- c[rad > beta] = coeff * np.power(s1, power)
-
- # Rad outside alpha-gamma interval
- if flipped is False:
- c[rad > gamma] = 1
- c[rad < alpha] = 0
- if flipped is True:
- c[rad > gamma] = 0
- c[rad < alpha] = 1
-
- c[c > 1] = 1
- c[c < 0] = 0
-
- confidence = c.reshape(radshape)
-
- 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 test():
rad = np.random.randint(50, size=[4, 8])
- #coeffs = [5, 42, 20, 28, 15, 35, 1]
- #coeffs = [20, 28, 5, 42, 15, 35, 1]
+ # coeffs = [5, 42, 20, 28, 15, 35, 1]
+ # coeffs = [20, 28, 5, 42, 15, 35, 1]
coeffs = [35, 15, 20, 1, 1]
- #confidence = conf_test_dble(rad, coeffs)
+ # confidence = conf_test_dble(rad, coeffs)
confidence = test_11um(rad, coeffs)
print(rad)
print('\n')
diff --git a/utils.py b/utils.py
new file mode 100644
index 0000000000000000000000000000000000000000..5be3d78aad42764e0c96c6288027983243a5f1c5
--- /dev/null
+++ b/utils.py
@@ -0,0 +1,225 @@
+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'][3]] = 1
+ sunglint_flag[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):
+
+ band = f'band{band_n}'
+
+# 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]:
+ sunglint = thresholds[f'{band}_0deg']
+
+ else:
+
+ if (refang > thresholds['bounds'][1] and refang <= thresholds['bounds'][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]):
+ lo_ang = thresholds['bounds'][2]
+ hi_ang = thresholds['bounds'][3]
+ lo_ang_val = thresholds[f'{band}_20deg'][0]
+ hi_ang_val = sunglint[1]
+ power = thresholds[f'{band}_20deg'][3]
+ conf_range = thresholds[f'{band}_20deg'][2]
+
+ 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
+
+ return sunglint
+
+
+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]
+ radshape = rad.shape
+ rad = rad.reshape((rad.shape[0]*rad.shape[1]))
+ c = 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:
+ c[rad <= beta] = coeff * np.power(s1, power)
+ if flipped is True:
+ c[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))
+ if flipped is True:
+ c[rad > beta] = coeff * np.power(s1, power)
+
+ # Rad outside alpha-gamma interval
+ if flipped is False:
+ c[rad > gamma] = 1
+ c[rad < alpha] = 0
+ if flipped is True:
+ c[rad > gamma] = 0
+ c[rad < alpha] = 1
+
+ c[c > 1] = 1
+ c[c < 0] = 0
+
+ confidence = c.reshape(radshape)
+
+ 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