finally implemented the Visible_Thresholds_Test. It definitely needs a closer...

finally implemented the Visible_Thresholds_Test. It definitely needs a closer look to make sure everything works as intended

main.py

@@ -195,7 +195,7 @@ def main(*, data_path=_datapath, mod02=_fname_mod02, mod03=_fname_mod03,
             pass
 
     if perform['11um BT Variability Test'] is True:
-        for scene_name in ['Polar_Night_Ocean']:
+        for scene_name in ['Ocean_Night', 'Polar_Night_Ocean']:
             pass
 
     if perform['11-4um BTD Oceanic Stratus'] is True:

preprocess_thresholds.py

@@ -4,6 +4,8 @@ import xarray as xr
 import ancillary_data as anc
 import utils
 
+from numpy.lib.stride_tricks import sliding_window_view
+
 _dtr = np.pi/180
 
 
@@ -171,10 +173,142 @@ def preproc_surf_temp(data, thresholds):
     return thr_out
 
 
-def get_b1_thresholds():
-    # fill_ndvi[0] is for fill_ndvi1
-    # fill_ndvi[1] is for fill_ndvi2
-    pass
+def var_11um(data, thresholds):
+    rad = data.M15.values
+    var = np.zeros((rad.shape[0], rad.shape[1], 9))
+    var_thr = thresholds['Daytime_Ocean_Spatial_Variability']['dovar11']
+
+    test = sliding_window_view(np.pad(rad, [1, 1], mode='constant'), (3, 3)) - np.expand_dims(rad, (2, 3))
+
+    var[np.abs(test).reshape(rad.shape[0], rad.shape[1], 9) < var_thr] = 1
+    var = var.sum(axis=2)
+
+    return var
+
+
+def get_b1_thresholds(data, thresholds):
+    ndvi = data.ndvi.values.reshape(data.ndvi.shape[0]*data.ndvi.shape[1])
+    sctang = data.scattering_angle.values.reshape(data.ndvi.shape[0]*data.ndvi.shape[1])
+
+    # this is hardcoded in the function
+    delta_ndvi_bin = 0.1
+
+    des_ndvi = thresholds['Misc']['des_ndvi']
+    thr_adj_fac_desert = thresholds['Misc']['adj_fac_desert']
+    thr_adj_fac_land = thresholds['Misc']['adj_fac_land']
+    ndvi_bnd1 = thresholds['Misc']['ndvi_bnd1']
+    ndvi_bnd2 = thresholds['Misc']['ndvi_bnd2']
+    fill_ndvi = thresholds['Misc']['fill_ndvi']
+
+    coeff1 = np.array(thresholds['Coeffs_Band1_land_thresh']).reshape(10, 3, 4)
+    coeff2 = np.zeros((10, 3, 4))
+    coeff2[:3, :, :] = np.array(thresholds['Coeffs_Band8_land_thresh']).reshape(3, 3, 4)
+    coeff = np.stack((coeff1, coeff2))
+
+    indvi = np.zeros(ndvi.shape)
+    indvi[ndvi >= ndvi_bnd2] = 9
+
+    x, y2 = np.zeros(ndvi.shape), np.zeros(ndvi.shape)
+
+    # this is equivalent to interp=1 in the C code
+    idx = np.nonzero((ndvi >= ndvi_bnd1) & (ndvi < ndvi_bnd2))
+
+    indvi[idx] = (ndvi[idx]/delta_ndvi_bin) - 0.5
+    indvi[ndvi < 0] = 0
+
+    x1 = delta_ndvi_bin*indvi + delta_ndvi_bin/2.0
+    x2 = x1 + delta_ndvi_bin
+    x[idx] = (ndvi[idx] - x1[idx])/(x2[idx] - x1[idx])
+    x = np.clip(x, 0, 1)
+
+    indvi = np.array(indvi, dtype=np.int)
+    thr = np.empty((ndvi.shape[0], 4))
+    thr_adj = np.empty((ndvi.shape[0], 4))
+
+    for i in range(3):
+        y1 = coeff[0, indvi, i, 0] + coeff[0, indvi, i, 1]*sctang + \
+             coeff[0, indvi, i, 2]*sctang**2 + coeff[0, indvi, i, 3]*sctang**3
+        des = np.nonzero(ndvi < des_ndvi)
+        y1[des] = coeff[1, indvi[des], i, 0] + coeff[1, indvi[des], i, 1]*sctang[des] + \
+            coeff[1, indvi[des], i, 2]*sctang[des]**2 + coeff[1, indvi[des], i, 3]*sctang[des]**3
+
+        y2[idx] = coeff[0, indvi[idx], i, 0] + \
+            coeff[0, indvi[idx], i, 1]*sctang[idx] + \
+            coeff[0, indvi[idx], i, 2]*sctang[idx]**2 + \
+            coeff[0, indvi[idx], i, 3]*sctang[idx]**3
+        idxdes = np.nonzero((ndvi >= ndvi_bnd1) & (ndvi < ndvi_bnd2) & (ndvi < des_ndvi))
+        y2[idxdes] = coeff[0, indvi[idxdes], i, 0] + \
+            coeff[0, indvi[idxdes], i, 1]*sctang[idxdes] + \
+            coeff[0, indvi[idxdes], i, 2]*sctang[idxdes]**2 + \
+            coeff[0, indvi[idxdes], i, 3]*sctang[idxdes]**3
+
+        thr[:, i] = (1.0 - x) + (x + y2)
+        thr_adj[:, i] = thr[:, i] * thr_adj_fac_desert
+        thr_adj[ndvi >= des_ndvi, i] = thr[ndvi >= des_ndvi, i] * thr_adj_fac_land
+
+    hicut = ((thr[:, 0] + thr_adj[:, 0])/100)  # .reshape(data.ndvi.shape)
+    midpt = ((thr[:, 1] + thr_adj[:, 1])/100)  # .reshape(data.ndvi.shape)
+    locut = ((thr[:, 2] + thr_adj[:, 2])/100)  # .reshape(data.ndvi.shape)
+
+    idx = np.nonzero((ndvi >= fill_ndvi[0]) | (ndvi <= fill_ndvi[1]))
+    hicut[idx] = -999
+    midpt[idx] = -999
+    locut[idx] = -999
+
+#    out_thr = xr.DataArray(data=np.dstack((locut, midpt, hicut, np.ones(data.ndvi.shape),
+#                                           np.full(data.ndvi.shape, 2))),
+#                           dims=('number_of_lines', 'number_of_pixels', 'z'))
+#
+#    return out_thr
+    return locut, midpt, hicut
+
+
+def vis_refl_thresholds(data, thresholds, scene):
+
+    locut, midpt, hicut = get_b1_thresholds(data, thresholds)
+    bias_adj = thresholds[scene]['Visible_Reflectance_Test']['adj']
+    ndvi = data.ndvi.values.reshape(data.ndvi.shape[0]*data.ndvi.shape[1])
+    m01 = data.M05.values.reshape(data.ndvi.shape[0]*data.ndvi.shape[1])
+    m02 = data.M07.values.reshape(data.ndvi.shape[0]*data.ndvi.shape[1])
+    m08 = data.M01.values.reshape(data.ndvi.shape[0]*data.ndvi.shape[1])
+
+    m128 = m01
+
+    b1_locut = locut * bias_adj
+    b1_midpt = midpt * bias_adj
+    b1_hicut = hicut * bias_adj
+
+    if ((scene == 'Land_Day_Desert') | (scene == 'Land_Day_Desert_Coast')):
+        ndvi_desert_thr = thresholds[scene]['Visible_Reflectance_Test']['ndvi_thr']
+        idx = np.nonzero(ndvi < ndvi_desert_thr)
+
+        b1_locut[idx] = locut[idx]
+        b1_midpt[idx] = midpt[idx]
+        b1_hicut[idx] = hicut[idx]
+        m128[idx] = m08[idx]
+
+    b1_power = np.full(b1_locut.shape, 2)
+
+    idx = np.nonzero(locut == -999)
+    b1_locut[idx] = thresholds[scene]['Visible_Reflectance_Test']['thr'][0]
+    b1_midpt[idx] = thresholds[scene]['Visible_Reflectance_Test']['thr'][1]
+    b1_hicut[idx] = thresholds[scene]['Visible_Reflectance_Test']['thr'][2]
+    b1_power[idx] = thresholds[scene]['Visible_Reflectance_Test']['thr'][3]
+    m128[idx] = m02[idx]
+
+    cosvza = np.cos(data.sensor_zenith.values * _dtr).reshape(ndvi.shape)
+
+    vzcpow = thresholds['VZA_correction']['vzcpow'][0]
+    b1_locut = (b1_locut * (1.0 / np.power(cosvza, vzcpow))).reshape(data.ndvi.shape)
+    b1_midpt = (b1_midpt * (1.0 / np.power(cosvza, vzcpow))).reshape(data.ndvi.shape)
+    b1_hicut = (b1_hicut * (1.0 / np.power(cosvza, vzcpow))).reshape(data.ndvi.shape)
+
+    out_thr = xr.DataArray(data=np.dstack((b1_locut, b1_midpt, b1_hicut, np.ones(data.ndvi.shape),
+                                           b1_power.reshape(data.ndvi.shape))),
+                           dims=('number_of_lines', 'number_of_pixels', 'z'))
+    out_rad = xr.DataArray(data=m128.reshape(data.M01.shape), dims=('number_of_lines', 'number_of_pixels'))
+
+    return out_thr, out_rad
 
 
 # NOTE: 11-12um Cirrus Test

tests.py

@@ -280,6 +280,10 @@ class CloudTests:
         elif test_name == 'NIR_Reflectance_Test':
             corr_thr = pt.preproc_nir(self.data, self.thresholds, self.scene_name)
             thr_xr['threshold'] = (('number_of_lines', 'number_of_pixels', 'z'), corr_thr)
+        elif test_name == 'Visible_Reflectance_Test':
+            thr_xr['threshold'], self.data['M128'] = pt.vis_refl_thresholds(self.data,
+                                                                            self.thresholds,
+                                                                            self.scene_name)
         else:
             thr_xr['threshold'] = (('number_of_lines', 'number_of_pixels', 'z'),
                                    np.ones((self.data[band].shape[0], self.data[band].shape[1], 5))*thr)
@@ -289,6 +293,11 @@ class CloudTests:
             data['sfcdif'] = (('number_of_lines', 'number_of_pixels'),
                               pt.preproc_sst(data, self.thresholds[self.scene_name][test_name]).values)
             band = 'sfcdif'
+        if test_name == '11um_Variability_Test':
+            var = pt.var_11um(self.data, self.thresholds)
+            data['11um_var'] = data.M15
+            data['11um_var'].values[var != 9] = np.nan
+
         if test_name == 'NIR_Reflectance_Test':
             pass
 

thresholds.mvcm.snpp.v0.0.1.yaml

@@ -12,12 +12,13 @@ Land_Day:
   11-4um_Oceanic_Stratus_Test: [-16.0, -14.0, -12.0, 1.0, 1.0]
   CO2_High_Clouds_Test: [222.0, 224.0, 226.0, 1.0, 1.0]
   Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0]
-  vis_refl_test  : [0.207, 0.169, 0.132, 1.0, 1.0] # used to be dlref1
+  Visible_Reflectance_Test:
+    thr: [0.207, 0.169, 0.132, 1.0, 1.0] # used to be dlref1
+    adj: 0.94 # used to be dl_b1bias_adj
   1.38um_High_Cloud_Test: [0.0375,  0.0250, 0.0125,  1.0, 1.0]
   dlvrat         : [1.80, 1.85, 1.90, 1.0, 1.0]
 #  dl11_12lcmult  : 0.3
 #  dl11_12hcadj   : 1.25
-  b1bias_adj     : 0.94
 
 Land_Night:
   11-12um_Cirrus_Test:
@@ -59,7 +60,9 @@ Land_Day_Coast:
   11-4um_Oceanic_Stratus_Test: [-16.0, -14.0, -12.0, 1.0, 1.0]
   CO2_High_Clouds_Test: [222.0, 224.0, 226.0, 1.0, 1.0]
   Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0]
-  dlref1_t2      : [0.207, 0.169, 0.132, 1.0, 1.0]
+  Visible_Reflectance_Test:
+    thr: [0.207, 0.169, 0.132, 1.0, 1.0]
+    adj: 0.94
   1.38um_High_Cloud_Test: [0.0375,  0.0250, 0.0125,  1.0, 1.0]
 # dl11_12lcmult_t2 : 0.3
 # dl11_12hcadj_t2  : 1.25
@@ -74,14 +77,16 @@ Land_Day_Desert:
   lds11_4lo      : [-28.0, -26.0, -24.0, 1.00]
   CO2_High_Clouds_Test: [222.0, 224.0, 226.0, 1.0, 1.0]
   Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0]
-  ldsref2        : [0.326, 0.288, 0.250, 1.00, 1.0]
+  Visible_Reflectance_Test:
+    thr: [0.326, 0.288, 0.250, 1.00, 1.0]
+    adj: 0.94
+    ndvi_thr: 0.25  # this used to be lds_ndvi
   1.38um_High_Cloud_Test: [0.0375, 0.0250, 0.0125, 1.00, 1.0]
   ldsgemi0       : [0.085, 0.095, 0.115, 1.00, 1.0]
   ldsgemi1       : [0.145, 0.170, 0.220, 1.00]
   ldsgemi2       : [0.310, 0.335, 0.360, 1.00]
   lds_ref3_tpw   : -10.00
   ldsbt1         : 320.0
-  lds_ndvi       : 0.25
 # lds11_12lcmult : 0.3
 # lds11_12hcadj  : 1.25
 
@@ -95,11 +100,13 @@ Land_Day_Desert_Coast:
   lds11_4lo_c    : [-23.0, -21.0, -19.0, 1.00, 1.0]
   CO2_High_Clouds_Test: [222.0, 224.0, 226.0, 1.0, 1.0]
   Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0]
-  ldsref2_c      : [0.326, 0.288, 0.250, 1.00, 1.0]
+  Visible_Reflectance_Test:
+    thr: [0.326, 0.288, 0.250, 1.00, 1.0]
+    adj: 0.94
+    ndvi_thr: 0.25
   1.38um_High_Cloud_Test: [0.0375, 0.0250, 0.0125,  1.00, 1.0]
   lds_ref3_tpw_c : -10.00
   ldsbt1_c       : 320.0
-  lds_ndvi_c     : 0.25
 # lds11_12lcmult_c : 0.3
 # lds11_12hcadj_c  : 1.25
 
@@ -154,7 +161,7 @@ Ocean_Night:
   CO2_High_Clouds_Test: [222.0, 224.0, 226.0, 1.0, 1.0]
   Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0]
   no86_73        : [14.0,  15.0,  16.0,  1.0, 1.0]
-  11um_var       : [3.0,   6.0,   7.0,   1.0, 1.0]
+  11um_Variability_Test: [3.0, 6.0, 7.0, 1.0, 1.0]
   8.6-11um_Test: [-0.5, -1.0,  -1.5,   1.0, 1.0]
   SST_Test:
     thr: [3.000, 2.500, 2.000, 1.0, 1.0]
@@ -172,8 +179,10 @@ Polar_Day_Land:
     adj: 0.3
   pdl_ref3_tpw   : -10.00
   11-4um_Oceanic_Stratus_Test: [-16.0, -14.0, -12.0, 1.0, 1.0]
-  pdlh20         : [215.0, 220.0, 225.0, 1.0, 1.0]
-  pdlref1        : [0.207, 0.169, 0.132, 1.0, 1.0]
+  pdlh20       : [215.0, 220.0, 225.0, 1.0, 1.0]
+  Visible_Reflectance_Test:
+    thr: [0.207, 0.169, 0.132, 1.0, 1.0]
+    adj: 0.94
   1.38um_High_Cloud_Test: [0.0375,  0.0250, 0.0125,  1.0, 1.0]
 
 Polar_Night_Land:
@@ -209,7 +218,9 @@ Polar_Day_Coast:
   pdl_ref3_tpw_t2 : -10.00
   11-4um_Oceanic_Stratus_Test: [-16.0, -14.0, -12.0, 1.0, 1.0]
   Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0]
-  pdlref1_t2      : [0.207, 0.169, 0.132, 1.0, 1.0]
+  Visible_Reflectance_Test:
+    thr: [0.207, 0.169, 0.132, 1.0, 1.0]
+    adj: 0.94
   1.38um_High_Cloud_Test: [0.0375,  0.0250, 0.0125,  1.0, 1.0]
 
 Polar_Day_Desert:
@@ -221,7 +232,9 @@ Polar_Day_Desert:
   pds11_4lo      : [-22.0, -20.0, -18.0, 1.00]
   pds11_4_pfm    : 1.0
   Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0]
-  pdsref2        : [0.326, 0.288, 0.250, 1.00, 1.0]
+  Visible_Reflectance_Test:
+    thr: [0.326, 0.288, 0.250, 1.00, 1.0]
+    adj: 0.94
   1.38um_High_Cloud_Test: [0.0375, 0.0250, 0.0125,  1.00, 1.0]
   pdsgemi0       : [0.085, 0.110, 0.135, 1.00, 1.0]
   pdsgemi1       : [0.170, 0.220, 0.270, 1.00]
@@ -238,7 +251,9 @@ Polar_Day_Desert_Coast:
   pds11_4lo_c    : [-23.0, -21.0, -19.0, 1.00]
   pds11_4_c_pfm  : 1.0
   Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0]
-  pdsref2_c      : [0.326, 0.288, 0.250, 1.00, 1.0]
+  Visible_Reflectance_Test:
+    thr: [0.326, 0.288, 0.250, 1.00, 1.0]
+    adj: 0.94
   1.38um_High_Cloud_Test: [0.0375, 0.0250, 0.0125,  1.00, 1.0]
   pds_ref3_tpw_c : -10.00
   pdsbt1_c       : 320.0
@@ -341,7 +356,7 @@ Polar_Ocean_Night:
   11um_Test: [267.0, 270.0, 273.0, 1.0, 1.0]
   Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0]
   pno86_73       : [14.0,  15.0,  16.0,  1.0, 1.0]
-  pno_11var      : [3.0,   6.0,   7.0,   1.0, 1.0]
+  11um_Variability_Test: [3.0,   6.0,   7.0,   1.0, 1.0]
   8.6-11um_Test: [-0.5,  -1.0,  -1.5,   1.0, 1.0]
   SST_Test:
     thr: [3.000, 2.500, 2.000, 1.0, 1.0]