diff --git a/preprocess_thresholds.py b/preprocess_thresholds.py
index b0c96f7fa374254b322f4d24118fc5b299a692fc..9d51feffd470b63c389485fd7f57f5f90da8091f 100644
--- a/preprocess_thresholds.py
+++ b/preprocess_thresholds.py
@@ -127,18 +127,18 @@ def thresholds_NIR(data, thresholds, scene, test_name, scene_idx):
return corr_thr
-def preproc_surf_temp(data, thresholds):
- thr_sfc1 = thresholds['Surface_Temperature_Test_1']
- thr_sfc2 = thresholds['Surface_Temperature_Test_2']
- thr_df1 = thresholds['Surface_Temperature_Test_df1']
- thr_df2 = thresholds['Surface_Temperature_Test_df2']
+def thresholds_surface_temperature(data, thresholds, scene_idx):
+ # def preproc_surf_temp(data, thresholds):
+ thr_sfc1 = thresholds['desert_thr']
+ thr_sfc2 = thresholds['regular_thr']
+ thr_df1 = thresholds['channel_diff_11-12um_thr']
+ thr_df2 = thresholds['channel_diff_11-4um_thr']
max_vza = 70.13 # This values is set based on sensor. Check mask_processing_constants.h for MODIS value
- rs = np.prod(data.M15.shape)
- df1 = (data.M15 - data.M16).values.reshape(rs)
- df2 = (data.M15 - data.M13).values.reshape(rs)
- desert_flag = data.Desert.values.reshape(rs)
- thresh = np.ones((rs, )) * thr_sfc1
+ df1 = (data.M15 - data.M16).values[scene_idx].ravel()
+ df2 = (data.M15 - data.M13).values[scene_idx].ravel()
+ desert_flag = data.Desert.values[scene_idx].ravel()
+ thresh = np.ones(df1.shape) * thr_sfc1
idx = np.where((df1 >= thr_df1[0]) | ((df1 < thr_df1[0]) & ((df2 <= thr_df2[0]) | (df2 >= thr_df2[1]))))
thresh[idx] = thr_sfc2
@@ -149,15 +149,14 @@ def preproc_surf_temp(data, thresholds):
idx = np.where(df1 >= thr_df1[1])
midpt[idx] = thresh[idx] + 2.0*df1[idx]
- corr = np.power(data.sensor_zenith.values/max_vza, 4) * 3.0
- midpt = midpt.reshape(corr.shape) + corr
+ corr = np.power(data.sensor_zenith.values[scene_idx].ravel()/max_vza, 4) * 3.0
+ midpt = midpt + corr
locut = midpt + 2.0
hicut = midpt - 2.0
- 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'))
+ thr_out = np.dstack((locut, midpt, hicut, np.ones(locut.shape), np.ones(locut.shape)))
- return thr_out
+ return np.squeeze(thr_out.T)
# This function is currently not used
@@ -455,18 +454,17 @@ def gemi_thresholds(data, thresholds, scene_name, scene_idx):
return gemi_thr.T
-def bt11_4um_preproc(data, thresholds, scene_name):
- thresh = thresholds[scene_name]['11-4um_BT_Difference_Test']
- c = thresh['coeffs']
+def bt_diff_11_4um_thresholds(data, threshold):
+ c = threshold['coeffs']
tpw = data.geos_tpw.values
- thr = c[0] + thresh['corr'] + c[1]*tpw + c[2]*np.power(tpw, 2)
- hicut0 = (thr + thresh['hicut_coeff'][0]).reshape(1, np.prod(tpw.shape))
- hicut1 = (thr + thresh['hicut_coeff'][1]).reshape(1, np.prod(tpw.shape))
- midpt0 = (hicut0 + thresh['midpt_coeff'][0]).reshape(1, np.prod(tpw.shape))
- midpt1 = (hicut1 + thresh['midpt_coeff'][1]).reshape(1, np.prod(tpw.shape))
- locut0 = (hicut0 + thresh['locut_coeff'][0]).reshape(1, np.prod(tpw.shape))
- locut1 = (hicut1 + thresh['locut_coeff'][1]).reshape(1, np.prod(tpw.shape))
+ thr = c[0] + threshold['corr'] + c[1]*tpw + c[2]*np.power(tpw, 2)
+ hicut0 = (thr + threshold['hicut_coeff'][0]).reshape(1, np.prod(tpw.shape))
+ hicut1 = (thr + threshold['hicut_coeff'][1]).reshape(1, np.prod(tpw.shape))
+ midpt0 = (hicut0 + threshold['midpt_coeff'][0]).reshape(1, np.prod(tpw.shape))
+ midpt1 = (hicut1 + threshold['midpt_coeff'][1]).reshape(1, np.prod(tpw.shape))
+ locut0 = (hicut0 + threshold['locut_coeff'][0]).reshape(1, np.prod(tpw.shape))
+ locut1 = (hicut1 + threshold['locut_coeff'][1]).reshape(1, np.prod(tpw.shape))
thr_out = np.vstack([hicut0, midpt0, locut0, locut1, midpt1, hicut1,
np.ones(hicut0.shape), np.ones(hicut0.shape)])
diff --git a/read_data.py b/read_data.py
index 24ec026fecf7b2fc9096db09577f3351b054f2bb..96f61b17dbd4075f7064b96f4ca511a72cb85ac3 100644
--- a/read_data.py
+++ b/read_data.py
@@ -150,24 +150,37 @@ def get_data(file_names: Dict[str, str],
mod02 = file_names['MOD02']
mod03 = file_names['MOD03']
- img02 = file_names['IMG02']
- img03 = file_names['IMG03']
+
+ if hires is True:
+ img02 = file_names['IMG02']
+ img03 = file_names['IMG03']
if hires is False:
viirs_data = read_data('viirs', f'{mod02}', f'{mod03}')
viirs_data = read_ancillary_data(file_names, viirs_data)
- m01 = viirs_data.M05.values
- m02 = viirs_data.M07.values
- r1 = 2.0 * (np.power(m02, 2.0) - np.power(m01, 2.0)) + (1.5 * m02) + (0.5 * m01)
- r2 = m02 + m01 + 0.5
- r3 = r1 / r2
- gemi = r3 * (1.0 - 0.25*r3) - ((m01 - 0.125) / (1.0 - m01))
-
- idx = np.nonzero((viirs_data.M05.values < -99) | (viirs_data.M05.values > 2))
- viirs_data['M05'].values[idx] = _bad_data
- idx = np.nonzero((viirs_data.M07.values < -99) | (viirs_data.M07.values > 2))
- viirs_data['M07'].values[idx] = _bad_data
+ if (('M05' in viirs_data) and ('M07' in viirs_data)):
+ m01 = viirs_data.M05.values
+ m02 = viirs_data.M07.values
+ r1 = 2.0 * (np.power(m02, 2.0) - np.power(m01, 2.0)) + (1.5 * m02) + (0.5 * m01)
+ r2 = m02 + m01 + 0.5
+ r3 = r1 / r2
+ gemi = r3 * (1.0 - 0.25*r3) - ((m01 - 0.125) / (1.0 - m01))
+ else:
+ gemi = np.full((viirs_data.M15.shape), _bad_data)
+
+ if 'M05' in viirs_data:
+ idx = np.nonzero((viirs_data.M05.values < -99) | (viirs_data.M05.values > 2))
+ viirs_data['M05'].values[idx] = _bad_data
+ else:
+ viirs_data['M05'] = (('number_of_lines', 'number_of_pixels'),
+ np.full(viirs_data.M15.shape, _bad_data))
+ if 'M07' in viirs_data:
+ idx = np.nonzero((viirs_data.M07.values < -99) | (viirs_data.M07.values > 2))
+ viirs_data['M07'].values[idx] = _bad_data
+ else:
+ viirs_data['M07'] = (('number_of_lines', 'number_of_pixels'),
+ np.full(viirs_data.M15.shape, _bad_data))
idx = np.nonzero((viirs_data.M12.values < 0) | (viirs_data.M12.values > 1000))
viirs_data['M12'].values[idx] = _bad_data
@@ -180,7 +193,7 @@ def get_data(file_names: Dict[str, str],
idx = np.nonzero((viirs_data.M16.values < 0) | (viirs_data.M16.values > 1000))
viirs_data['M16'].values[idx] = _bad_data
- # Compute channel differences and ratios that are used in the tests
+ # Compute channel differences and ratios that are used in the tests
viirs_data['M14-M15'] = (('number_of_lines', 'number_of_pixels'),
viirs_data.M14.values - viirs_data.M15.values)
viirs_data['M15-M16'] = (('number_of_lines', 'number_of_pixels'),
@@ -194,7 +207,7 @@ def get_data(file_names: Dict[str, str],
viirs_data['GEMI'] = (('number_of_lines', 'number_of_pixels'), gemi)
# temp value to force the code to work
- viirs_data['M128'] = (('number_of_lines', 'number_of_pixels'), np.zeros(viirs_data.M01.shape))
+ viirs_data['M128'] = (('number_of_lines', 'number_of_pixels'), np.zeros(viirs_data.M15.shape))
else:
viirs_data = read_data('viirs', f'{img02}', f'{img03}')
diff --git a/spectral_tests.py b/spectral_tests.py
index 772448beaa4e3a5b6696cff501ea236883cc0724..d08c6001af6822ac4d8933db8f614a12c80edb8e 100644
--- a/spectral_tests.py
+++ b/spectral_tests.py
@@ -76,11 +76,28 @@ class CloudTests(object):
@run_if_test_exists_for_scene
def surface_temperature_test(self,
- band31: str,
- band32: str,
+ band: str,
+ viirs_data: xr.Dataset,
cmin: np.ndarray,
- test_name: str = 'Surface_Temperature_test') -> np.ndarray:
- pass
+ test_name: str = 'Surface_Temperature_Test') -> np.ndarray:
+
+ confidence = np.ones(self.data[band].shape)
+ qa_bit = np.zeros(self.data[band].shape)
+ test_bit = np.zeros(self.data[band].shape)
+ threshold = self.thresholds[self.scene_name][test_name]
+
+ if (threshold['perform'] is True and self.pixels_in_scene is True):
+ qa_bit[self.scene_idx] = 1
+ print(f'Testing "{self.scene_name}"\n')
+ rad = self.data[band].values[self.scene_idx]
+ sfcdif = viirs_data.geos_sfct.values[self.scene_idx] - rad
+ # need to write the test_bit here
+ thr = preproc.thresholds_surface_temperature(viirs_data, threshold, self.scene_idx)
+ confidence[self.scene_idx] = conf.conf_test_new(sfcdif, thr)
+
+ cmin = np.fmin(cmin, confidence)
+
+ return cmin, test_bit
@run_if_test_exists_for_scene
def sst_test(self,
@@ -115,7 +132,7 @@ class CloudTests(object):
cmin = np.fmin(cmin, confidence)
# return cmin, np.abs(1-test_bit)*qa_bit
- return sfcdif, test_bit
+ return cmin, test_bit
@run_if_test_exists_for_scene
def bt_diff_86_11um(self,
@@ -139,7 +156,7 @@ class CloudTests(object):
cmin = np.fmin(cmin, confidence)
- return cmin, np.abs(1-test_bit)*qa_bit
+ return cmin, test_bit # np.abs(1-test_bit)*qa_bit
@run_if_test_exists_for_scene
def test_11_12um_diff(self,
@@ -169,6 +186,36 @@ class CloudTests(object):
# return cmin, np.abs(1-test_bit)*qa_bit
return cmin, test_bit
+ @run_if_test_exists_for_scene
+ def bt_difference_11_4um_test(self,
+ band: str,
+ cmin: np.ndarray,
+ test_name: str = '11-4um_BT_Difference_Test') -> np.ndarray:
+
+ confidence = np.ones(self.data.M01.shape)
+ qa_bit = np.zeros(self.data[band].shape)
+ test_bit = np.zeros(self.data[band].shape)
+ threshold = self.thresholds[self.scene_name][test_name]
+
+ if (threshold['perform'] is True and self.pixels_in_scene is True):
+ qa_bit[self.scene_idx] = 1
+ thr = preproc.bt_diff_11_4um_thresholds(self.data, threshold)
+
+ # CONTINUE FROM HERE...
+
+
+ @run_if_test_exists_for_scene
+ def midlevel_cloud_test():
+ pass
+
+ @run_if_test_exists_for_scene
+ def water_vapor_cloud_test():
+ pass
+
+ @run_if_test_exists_for_scene
+ def variability_11um_test():
+ pass
+
@run_if_test_exists_for_scene
def oceanic_stratus_11_4um_test(self,
band: str,
@@ -205,7 +252,7 @@ class CloudTests(object):
cmin = np.fmin(cmin, confidence)
- return cmin, np.abs(1-test_bit)*qa_bit
+ return cmin, test_bit # np.abs(1-test_bit)*qa_bit
@run_if_test_exists_for_scene
def nir_reflectance_test(self,
diff --git a/thresholds.mvcm.snpp.v0.0.1.yaml b/thresholds.mvcm.snpp.v0.0.1.yaml
index 5f5817e059e597aeae0a67b8f54210d3febe261d..e8f59044f4b192e4f94bf6320701eb16aebe1ad2 100644
--- a/thresholds.mvcm.snpp.v0.0.1.yaml
+++ b/thresholds.mvcm.snpp.v0.0.1.yaml
@@ -45,14 +45,20 @@ Land_Night:
nl_11_4h : [6.5, 6.0, 5.5, 1.0]
nl_11_4m : [-0.5, 6.0, 0.5, 1.0]
nl_11_4_pfm : 1.0
- Surface_Temperature_Test_df1: [-0.2, 1.0]
- Surface_Temperature_Test_df2: [-0.5, 1.0]
- Surface_Temperature_Test_difference: [-0.2, 1.0, -0.5, 1.0] # <- this merges the previous two arrays
+ Surface_Temperature_Test:
+ desert_thr: 20.0
+ regular_thr: 12.0
+ channel_diff_11-12um_thr: [-0.2, 1.0]
+ channel_diff_11-4um_thr: [-0.5, 1.0]
+ perform: True
+ # Surface_Temperature_Test_df1: [-0.2, 1.0]
+ # Surface_Temperature_Test_df2: [-0.5, 1.0]
+ # Surface_Temperature_Test_difference: [-0.2, 1.0, -0.5, 1.0] # <- this merges the previous two arrays
bt_diff_bounds : [1.0, -1.0]
- Surface_Temperature_Test_1: 20.0 # | might be worth figuring out if we can
- Surface_Temperature_Test_2: 12.0 # | merge these three coefficients
- Surface_Temperature_Test_pfm: 1.0 # __|
- Surface_Temperature_Test: [20.0, 12.0, 1.0] # <- First attempt to merge the three values above
+ # Surface_Temperature_Test_1: 20.0 # | might be worth figuring out if we can
+ # Surface_Temperature_Test_2: 12.0 # | merge these three coefficients
+ # Surface_Temperature_Test_pfm: 1.0 # __|
+ # Surface_Temperature_Test: [20.0, 12.0, 1.0] # <- First attempt to merge the three values above
nlbt1 : 270.0
nl_lat : 30.0
nl_ndvi : 0.25
@@ -256,13 +262,19 @@ Polar_Night_Land:
bt1: 270.0
perform: True
pnlbt2 : 270.0
- Surface_Temperature_Test_df1: [0.0, 1.0]
- Surface_Temperature_Test_df2: [-0.5, 1.0]
- Surface_Temperature_Test_difference: [-0.2, 1.0, -0.5, 1.0] # <- this merges the previous two arrays
- Surface_Temperature_Test_1: 20.0 # | might be worth figuring out if we can
- Surface_Temperature_Test_2: 12.0 # | merge these three coefficients
- Surface_Temperature_Test_pfm: 1.0 # __|
- Surface_Temperature_Test: [20.0, 12.0, 1.0] # <- First attempt to merge the three values above
+ Surface_Temperature_Test:
+ desert_thr: 20.0
+ regular_thr: 12.0
+ channel_diff_11-12um_thr: [0, 1.0]
+ channel_diff_11-4um_thr: [-0.5, 1.0]
+ perform: True
+ # Surface_Temperature_Test_df1: [0.0, 1.0]
+ # Surface_Temperature_Test_df2: [-0.5, 1.0]
+ # Surface_Temperature_Test_difference: [-0.2, 1.0, -0.5, 1.0] # <- this merges the previous two arrays
+ # Surface_Temperature_Test_1: 20.0 # | might be worth figuring out if we can
+ # Surface_Temperature_Test_2: 12.0 # | merge these three coefficients
+ # Surface_Temperature_Test_pfm: 1.0 # __|
+ # Surface_Temperature_Test: [20.0, 12.0, 1.0] # <- First attempt to merge the three values above
pnl_11_4_pfm : 1.0
pnl_7_11_pfm : 1.0
pnl_4_12_pfm : 1.0