started implementing some more tests. results are tracked in excel file

main.py

@@ -84,10 +84,171 @@ def main(*, data_path=_datapath, mod02=_fname_mod02, mod03=_fname_mod03,
 #                 'Polar_Ocean_Day': np.ones(viirs_data.M01.shape),
 #                 'Polar_Ocean_Night': np.ones(viirs_data.M01.shape)
 #                 }
-    cmin2 = np.ones(viirs_data.M01.shape)
-    cmin3 = np.ones(viirs_data.M01.shape)
-    cmin4 = np.ones(viirs_data.M01.shape)
-
+#    cmin2 = np.ones(viirs_data.M01.shape)
+#    cmin3 = np.ones(viirs_data.M01.shape)
+#    cmin4 = np.ones(viirs_data.M01.shape)
+
+    # ------------------- #
+    # ### Testing Setup ###
+    # ------------------- #
+    perform = {'11um BT Test': False,
+               'CO2 High Clouds Test': False,
+               'Water Vapor High Clouds Test': False,
+               'Surface Temperature Test': False,
+               'SST Test': False,
+               '8.6-11um BT Difference Test': False,
+               '11-12um BTD Thin Cirrus Test': True,
+               '11-4um BT Difference Test': False,
+               '7.3-11um BT Difference Mid-level Clouds': False,
+               'Water Vapor Cloud Test': False,
+               '11um BT Variability Test': False,
+               '11-4um BTD Oceanic Stratus': False,
+               'NIR Reflectance Test': False,
+               'Vis/NIR Ratio Test': False,
+               '1.6um or 2.1um NIR Reflectance Test': False,
+               'Visible Reflectance Test': False,
+               'GEMI Test': False,
+               '1.38um High Cloud Test': False,
+               '4-12um BTD Thin Cirrus Test': False
+               }
+
+    # --------------------- #
+    # ### Group 1 Tests ### #
+    # --------------------- #
+
+    if perform['11um BT Test'] is True:
+        # 11um BT Test
+        for scene_name in ['Ocean_Day', 'Ocean_Night', 'Polar_Ocean_Day', 'Polar_Ocean_Night']:
+            SceneType = CloudTests(viirs_data, scene_name, thresholds)
+            cmin_G1 = SceneType.single_threshold_test('11um_Test', 'M15', cmin_G1)
+
+    if perform['CO2 High Clouds Test'] is True:
+        # CO2 High Clouds Test
+        for scene_name in ['Land_Day', 'Land_Night', 'Land_Day_Coast', 'Land_Day_Desert',
+                           'Land_Day_Desert_Coast', 'Ocean_Day', 'Ocean_Night', 'Day_Snow', 'Night_Snow']:
+            SceneType = CloudTests(viirs_data, scene_name, thresholds)
+            cmin_G1 = SceneType.single_threshold_test('CO2_High_Clouds_tests', 'bad_data', cmin_G1)
+
+    if perform['Water Vapor High Clouds Test'] is True:
+        # Water Vapor High Clouds Test
+        for scene_name in ['Land_Day', 'Land_Night', 'Land_Day_Coast', 'Land_Day_Desert',
+                           'Land_Day_Desert_Coast', 'Ocean_Day', 'Ocean_Night', 'Polar_Day_Land',
+                           'Polar_Night_Land', 'Polar_Day_Coast', 'Polar_Day_Desert',
+                           'Polar_Day_Desert_Coast', 'Polar_Day_Snow', 'Polar_Night_Snow', 'Polar_Ocean_Day',
+                           'Polar_Ocean_Night', 'Day_Snow', 'Night_Snow', 'Antarctic_Day']:
+            SceneType = CloudTests(viirs_data, scene_name, thresholds)
+            cmin_G1 = SceneType.single_threshold_test('Water_Vapor_High_Clouds_tests', 'bad_data', cmin_G1)
+
+    if perform['Surface Temperature Test'] is True:
+        # Surface Temperature Test
+        # NOTE: this test requires the thresholds to be preprocessed
+        for scene_name in ['Land_Night', 'Polar_Night_Land']:
+            SceneType = CloudTests(viirs_data, scene_name, thresholds)
+            cmin_G1 = SceneType.single_threshold_test('Surface_Temperature_Test', 'M15-M16', cmin_G1)
+
+    if perform['SST Test'] is True:
+        # SST Test
+        for scene_name in ['Ocean_Day', 'Ocean_Night', 'Polar_Ocean_Day', 'Polar_Ocean_Night']:
+            SceneType = CloudTests(viirs_data, scene_name, thresholds)
+            cmin_G1 = SceneType.single_threshold_test('SST_Test', 'M15-M16', cmin_G1)
+
+    # --------------------- #
+    # ### Group 2 tests ### #
+    # --------------------- #
+
+    if perform['8.6-11um BT Difference Test'] is True:
+        # 8.6-11um BT Difference Test
+        for scene_name in ['Ocean_Day', 'Ocean_Night', 'Polar_Ocean_Day', 'Polar_Ocean_Night']:
+            SceneType = CloudTests(viirs_data, scene_name, thresholds)
+            cmin_G2 = SceneType.single_threshold_test('8.6-11um_Test', 'M15-M14', cmin_G2)
+
+    if perform['11-12um BTD Thin Cirrus Test'] is True:
+        # 11-12um BT BTD Transmissive Cirrus Test
+        # NOTE: some of the tests have some differences in how the thresholds are derived
+        # The commented list scene_name is the complete list, the one currently in use is to test that
+        # the template code works at least with a subset of scenes that have the same way of deriving the
+        # thresholds
+        # for scene_name in ['Land_Day', 'Land_Day_Coast', 'Land_Day_Desert', 'Land_Day_Desert_Coast',
+        #                    'Ocean_Day', 'Ocean_Night', 'Polar_Day_Land', 'Polar_Day_Coast',
+        #                    'Polar_Day_Desert', 'Polar_Day_Desert_Coast', 'Polar_Day_Snow', 'Polar_Night_Snow',
+        #                    'Polar_Ocean_Day', 'Polar_Ocean_Night', 'Day_Snow', 'Night_Snow']:
+        for scene_name in ['Land_Day', 'Land_Day_Coast', 'Land_Day_Desert', 'Land_Day_Desert_Coast',
+                           'Ocean_Day', 'Ocean_Night', 'Polar_Ocean_Day', 'Polar_Ocean_Night']:
+            SceneType = CloudTests(viirs_data, scene_name, thresholds)
+            cmin_G2 = SceneType.single_threshold_test('11-12um_Cirrus_Test', 'M15-M16', cmin_G2)
+
+    if perform['11-4um BT Difference Test'] is True:
+        # 11-4um BT Difference Test
+        for scene_name in ['Ocean_Night', 'Polar_Ocean_Night']:
+            pass
+
+    if perform['7.3-11um BT Difference Mid-level Clouds'] is True:
+        for scene_name in ['Land_Night', 'Polar_Night_Land', 'Polar_Night_Snow', 'Night_Snow']:
+            pass
+
+    if perform['Water Vapor Cloud Test'] is True:
+        for scene_name in ['Polar_Night_Ocean']:
+            pass
+
+    if perform['11um BT Variability Test'] is True:
+        for scene_name in ['Polar_Night_Ocean']:
+            pass
+
+    if perform['11-4um BTD Oceanic Stratus'] is True:
+        # 11-4um BT Difference for oceanic stratus (low emissivity water cloud) Test
+        for scene_name in ['Land_Day', 'Land_Day_Coast', 'Land_Day_Desert', 'Land_Day_Desert_Coast',
+                           'Ocean_Day', 'Ocean_Night', 'Polar_Day_Land', 'Polar_Day_Coast',
+                           'Polar_Day_Desert', 'Polar_Day_Desert_Coast', 'Polar_Day_Snow', 'Polar_Night_Snow',
+                           'Polar_Ocean_Day', 'Polar_Ocean_Night', 'Day_Snow', 'Night_Snow']:
+            SceneType = CloudTests(viirs_data, scene_name, thresholds)
+            cmin_G2 = SceneType.single_threshold_test('11-4um_BTD_Oceanic_Stratus_Test', 'M15-M13', cmin_G2)
+
+    # --------------------- #
+    # ### Group 3 tests ### #
+    # --------------------- #
+
+    if perform['NIR Reflectance Test'] is True:
+        for scene_name in ['Ocean_Day', 'Polar_Ocean_Day']:
+            pass
+
+    if perform['Vis/NIR Ratio Test'] is True:
+        for scene_name in ['Ocean_Day', 'Polar_Ocean_Day']:
+            pass
+
+    if perform['1.6um or 2.1um NIR Reflectance Test'] is True:
+        for scene_name in ['Ocean_Day', 'Polar_Ocean_Day']:
+            pass
+
+    if perform['Visible Reflectance Test'] is True:
+        for scene_name in ['Land_Day', 'Land_Day_Coast', 'Land_Day_Desert', 'Land_Day_Desert_Coast',
+                           'Polar_Day_Land', 'Polar_Day_Coast', 'Polar_Day_Desert', 'Polar_Day_Desert_Coast']:
+            pass
+
+    if perform['GEMI Test'] is True:
+        for scene_name in ['Land_Day_Desert', 'Polar_Day_Desert']:
+            pass
+
+    # --------------------- #
+    # ### Group 4 tests ### #
+    # --------------------- #
+
+    if perform['1.38um High Cloud Test'] is True:
+        for scene_name in ['Land_Day', 'Land_Day_Coast', 'Land_Day_Desert', 'Land_Day_Desert_Coast',
+                           'Ocean_Day', 'Polar_Day_Land', 'Polar_Day_Coast', 'Polar_Day_Desert',
+                           'Polar_Day_Desert_Coast', 'Polar_OCean_Day', 'Day_Snow']:
+            pass
+
+    # --------------------- #
+    # ### Group 5 tests ### #
+    # --------------------- #
+
+    if perform['4-12um BTD Thin Cirrus Test'] is True:
+        for scene_name in ['Land_Night', 'Polar_Night_Land', 'Polar_Night_Snow', 'Night_Snow']:
+            pass
+
+    return cmin_G2
+
+    '''
     Land_Day = CloudTests(viirs_data, 'Land_Day', thresholds)
     Land_Night = CloudTests(viirs_data, 'Land_Night', thresholds)
     Land_Day_Coast = CloudTests(viirs_data, 'Land_Day_Coast', thresholds)
@@ -235,6 +396,7 @@ def main(*, data_path=_datapath, mod02=_fname_mod02, mod03=_fname_mod03,
              lat=viirs_data.latitude.values, lon=viirs_data.longitude.values)
 
     return confidence
+    '''
 
 
 def test_main():

preprocess_thresholds.py

+import numpy as np
+import xarray as xr
+
+import ancillary_data as anc
+
+_dtr = np.pi/180
+
+
+# this case is written for the 11-12um Cirrus Test for scenes that follow pattern 1 (see note below)
+def prepare_thresholds(data, thresholds):
+
+    coeff_values = np.empty((data.M01.shape[0], data.M01.shape[1], 2))
+    coeff_values[:, :, 0] = np.full(data.M01.shape, thresholds['11-12um_Cirrus_Test']['coeffs'][0])
+    coeff_values[:, :, 1] = np.full(data.M01.shape, thresholds['11-12um_Cirrus_Test']['coeffs'][1])
+    cmult_values = np.full(data.M01.shape, thresholds['11-12um_Cirrus_Test']['cmult'])
+    adj_values = np.full(data.M01.shape, thresholds['11-12um_Cirrus_Test']['adj'])
+
+    thr_dict = {'coeffs': (['number_of_lines', 'number_of_pixels', 'z'], coeff_values),
+                'cmult': (['number_of_lines', 'number_of_pixels'], cmult_values),
+                'adj': (['number_of_lines', 'number_of_pixels'], adj_values)
+                }
+
+    return xr.Dataset(data_vars=thr_dict)
+
+
+def preproc(data, thresholds):
+    cosvza = np.cos(data.sensor_zenith * _dtr)
+    schi = (1/cosvza).where(cosvza > 0, 99.0)
+
+    schi = schi.values.reshape(np.prod(schi.shape))
+    m15 = data.M15.values.reshape(np.prod(data.M15.shape))
+    thr = anc.py_cithr(1, schi, m15)
+    thr = thr.reshape(data.M15.shape)
+    schi = schi.reshape(data.M15.shape)
+
+    # thr_xr = xr.Dataset(np.full(data.sensor_zenith.shape, thresholds['coeffs']),
+    #                     dims=('number_of_lines', 'number_of_pixels'))
+    thr_xr = prepare_thresholds(data, thresholds)
+
+    midpt = thr_xr.coeffs[:, :, 0].where((thr < 0.1) | (np.abs(schi-99) < 0.0001), thr)
+    locut = midpt + (thr_xr.cmult * midpt)
+    hicut = midpt - thr_xr.adj
+
+    thr_out = xr.DataArray(data=np.dstack((locut, midpt, hicut, np.ones(locut.shape), np.ones(locut.shape))),
+                           dims=('number_of_lines', 'number_of_pixels', 'z'))
+    return thr_out
+    # return locut, hicut, midpt
+
+
+# NOTE: About the 11-12um Cirrus Test
+# hicut is computed in different ways depending on the scene
+# 1. midpt - adj
+# - Land_Day
+# - Land_Day_Coast
+# - Land_Day_Desert
+# - Land_Day_Desert_Coast
+# - Ocean_Day
+# - Ocean_Night
+# - Polar_Day_Ocean
+# - Polar_Night_Ocean
+#
+# 2. midpt - (btd_thr * adj)
+# - Polar_Day_Land
+# - Polar_Day_Coast
+# - Polar_Day_Desert
+# - Polar_Day_Desert_Coast
+# - Polar_Day_Snow
+#
+# 3. Others
+# - Land_Night
+# - Polar_Night_Land
+# - Polar_Night_Snow
+# - Day_Snow
+# - Night_Snow

read_data.py

@@ -136,6 +136,12 @@ def get_data(file_names, sunglint_angle):
     viirs_data = read_data('viirs', f'{mod02}', f'{mod03}')
     viirs_data = read_ancillary_data(file_names, viirs_data)
 
+    # Include channel differences and ratios that are used in the tests
+    viirs_data['M15-M14'] = (('number_of_lines', 'number_of_pixels'),
+                             viirs_data.M15.values - viirs_data.M14.values)
+    viirs_data['M15-M16'] = (('number_of_lines', 'number_of_pixels'),
+                             viirs_data.M15.values - viirs_data.M16.values)
+
     scene_flags = scn.find_scene(viirs_data, sunglint_angle)
     scene = scn.scene_id(scene_flags)
 

scene.py

@@ -10,7 +10,7 @@ import ancillary_data as anc
 _scene_list = ['Ocean_Day', 'Ocean_Night', 'Land_Day', 'Land_Night', 'Snow_Day', 'Snow_Night', 'Coast_Day',
                'Desert_Day', 'Antarctic_Day', 'Polar_Day_Snow', 'Polar_Day_Desert', 'Polar_Day_Ocean',
                'Polar_Day_Desert_Coast', 'Polar_Day_Coast', 'Polar_Day_Land', 'Polar_Night_Snow',
-               'Polar_Night_Land', 'Polar_Night_Ocean', 'Land_Day_Desert_Coast', 'Land_Day_Coast']
+               'Polar_Night_Land', 'Polar_Night_Ocean', 'Land_Day_Desert_Coast', 'Land_Day_Coast', 'Desert']
 _flags = ['day', 'night', 'land', 'coast', 'sh_lake', 'sh_ocean', 'water', 'polar', 'sunglint',
           'greenland', 'high_elevation', 'antarctica', 'desert', 'visusd', 'vrused', 'map_snow', 'map_ice',
           'ndsi_snow', 'snow', 'ice', 'new_zealand', 'uniform']
@@ -430,4 +430,7 @@ def scene_id(scene_flag):
                      (scene_flag['water'] == 1))
     scene['Polar_Night_Ocean'][idx] = 1
 
+    idx = np.nonzero(scene_flag['desert'] == 1)
+    scene['Desert'][idx] = 1
+
     return scene

tests.py

@@ -6,6 +6,7 @@ from numpy.lib.stride_tricks import sliding_window_view
 import utils
 import conf
 import conf_xr
+import preprocess_thresholds as pt
 
 
 # ############## GROUP 1 TESTS ############## #
@@ -212,7 +213,8 @@ class CloudMaskTests(object):
         pass
 
 
-class CloudTests:
+# old class, doesn't use xarray much
+class CloudTests_old:
 
     def __init__(self, scene_ids, scene_name, thresholds):
         self.scene = scene_ids
@@ -240,7 +242,7 @@ class CloudTests:
 
 
 # new class to try to use xarray more extensively
-class CloudTests_new:
+class CloudTests:
 
     def __init__(self, data, scene_name, thresholds):
         self.data = data
@@ -249,11 +251,19 @@ class CloudTests_new:
 
     def single_threshold_test(self, test_name, band, cmin):
 
+        if band == 'bad_data':
+            return cmin
+        print(f'Running test "{test_name}" for "{self.scene_name}"')
         # preproc_thresholds()
         thr = np.array(self.thresholds[self.scene_name][test_name])
         thr_xr = xr.Dataset()
-        thr_xr['threshold'] = (('number_of_lines', 'number_of_pixels', 'z'),
-                               np.ones((self.data[band].shape[0], self.data[band].shape[1], 5))*thr)
+        if test_name == '11-12um_Cirrus_Test':
+            thr_xr['threshold'] = pt.preproc(self.data, self.thresholds[self.scene_name])
+            thr = np.ones((5,))  # This is only temporary to force the logic of the code
+                                 # I need to find a better solution at some point
+        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)
         data = xr.Dataset(self.data, coords=thr_xr)
 
         if thr[4] == 1:
@@ -265,9 +275,29 @@ class CloudTests_new:
         return cmin
 
 
-# single_threshold_test('11BT', 'M15', cmin)
-# single_threshold_test('12-11BT', 'M16-M15', cmin)
-# single_threshold_test('12BT', 'M16', cmin)
+def preproc_thresholds(thresholds, data):
+    thr = np.array(thresholds)
+    thr_xr = xr.Dataset()
+    thr_xr['tresholds'] = (('number_of_lines', 'number_of_pixels', 'z'),
+                           np.ones((data['M01'].shape[0], data['M01'].shape[1], 5))*thr)
+
+    nl_sfct1 = thresholds['Land_Night']['Surface_Temperature_Test'][0]
+#    nl_sfct2 = thresholds['Land_Night']['Surface_Temperature_Test'][1]
+#    nlsfct_pfm = thresholds['Land_Night']['Surface_Temperature_Test'][2]
+    nl_df1 = thresholds['Land_Night']['Surface_Temperature_Test_difference'][0:2]
+    nl_df2 = thresholds['Land_Night']['Surface_Temperature_Test_difference'][2:]
+
+#    df1 = data.M15 - data.M16
+#    df2 = data.M15 - data.M13
+    thr_xr = thr_xr.where(data.desert != 1, nl_sfct1)
+    thr_xr = thr_xr.where((data['M15-M16'] > nl_df1[0]) |
+                          ((data['M15-M16'] < nl_df1[0]) &
+                           ((data['M15-M13'] <= nl_df2[0]) | (data['M15-M13'] >= nl_df2[1]))),
+                          nl_sfct1[0])
+
+    data = xr.Dataset(data, coords=thr_xr)
+
+    return data
 
 
 def single_threshold_test(test, rad, threshold):