reworked conf_test() to account for thresholds being 2d arrays in some cases

conf.py

+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
+
+
+

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
 

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

tests.py

 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):

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