diff --git a/main.py b/main.py
index e7851d4ff3dc633fc677a84e546db4b765e6b591..307186ca166d8907d9daf3d98fe2f384c347ae28 100644
--- a/main.py
+++ b/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():
diff --git a/preprocess_thresholds.py b/preprocess_thresholds.py
new file mode 100644
index 0000000000000000000000000000000000000000..af2251dd0a4dbefd5c27c803f36295788a609860
--- /dev/null
+++ b/preprocess_thresholds.py
@@ -0,0 +1,74 @@
+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
diff --git a/read_data.py b/read_data.py
index e95ec8693e0ed78810513c0cf1b991c6f341d5cd..087c32b4dda9ea021e3fd4cc6c9f13ce7bf7d804 100644
--- a/read_data.py
+++ b/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)
diff --git a/scene.py b/scene.py
index ddfac8812c820ed6f7443e53b6c13595027aefc9..474a2407c5541cf4cdc5d98d998da36fc0bc2291 100644
--- a/scene.py
+++ b/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
diff --git a/tests.py b/tests.py
index 39860e56c22cbf01147a4614a394f0a13065c975..fecbe3d4500d3fbd4591f03ab7ed4b2557aeccbb 100644
--- a/tests.py
+++ b/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):