diff --git a/.gitignore b/.gitignore
index 47273c4d3484caf6cf5b143c116469bc481d7e96..fbe75cf7eca557150c62dcfece25ad234b6365cf 100644
--- a/.gitignore
+++ b/.gitignore
@@ -4,6 +4,8 @@
*.npy
*.npz
+*.vscode
+
viz.py
build
diff --git a/conf.py b/conf.py
index 161789bc4e009e41532e79d95990a765a75b8260..ecda20600c1b6648892f0e0783eb4f804dc86478 100644
--- a/conf.py
+++ b/conf.py
@@ -23,8 +23,7 @@ def test_dble(flipped=False):
def conf_test_new(rad: np.ndarray,
thr: np.ndarray) -> np.ndarray:
- '''
- Assuming a linear function between min and max confidence level, the plot below shows
+ """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.
@@ -41,19 +40,18 @@ def conf_test_new(rad: np.ndarray,
e 0________/
| alpha
--------- radiance ---------->
- '''
+ """
radshape = rad.shape
rad = rad.ravel()
thr = np.array(thr)
-
if thr.ndim == 1:
thr = np.full((rad.shape[0], 4), thr[:4]).T
- coeff = np.power(2, (thr[3] - 1))
- hicut = thr[2, :]
- beta = thr[1, :]
+ coeff = np.power(2, (thr[3, :] - 1))
locut = thr[0, :]
+ beta = thr[1, :]
+ hicut = thr[2, :]
power = thr[3, :]
confidence = np.zeros(rad.shape)
@@ -185,7 +183,7 @@ def conf_test_dble(rad, coeffs):
coeffs = np.array(coeffs)
radshape = rad.shape
- rad = rad.reshape(np.prod(radshape))
+ rad = rad.ravel()
confidence = np.zeros(rad.shape)
alpha1, gamma1 = np.empty(rad.shape), np.empty(rad.shape)
@@ -267,10 +265,9 @@ def conf_test_dble(rad, coeffs):
confidence[(alpha1 - gamma1 <= 0) & ((rad > gamma1) | (rad < gamma2))] = 0
confidence[(alpha1 - gamma1 <= 0) & (rad <= alpha1) & (rad >= alpha2)] = 1
- confidence[confidence > 1] = 1
- confidence[confidence < 0] = 0
+ confidence = np.clip(confidence, 0, 1)
- return confidence
+ return confidence.reshape(radshape)
if __name__ == "__main__":
diff --git a/main.py b/main.py
index eba3e025a5cdce92818443e2d972c2d2330c70c9..2e0268febd8f5863d7b6c203d39f72868252bf68 100644
--- a/main.py
+++ b/main.py
@@ -1,17 +1,15 @@
-import ruamel_yaml as yml
+try:
+ from ruamel import yaml as yml
+except ImportError:
+ import ruamel_yaml as yml
+
import numpy as np
-# import xarray as xr
from glob import glob
import read_data as rd
-# import scene as scn
from tests import CloudTests
-# import tests
-# import ocean_day_tests as odt
-# import restoral
-
# #################################################################### #
# TEST CASE
# data:
diff --git a/preprocess_thresholds.py b/preprocess_thresholds.py
index 3f4ed92354e2bf530ad520a306964340ab584d81..413dbc68e231079918c969d41d506d16c709c53c 100644
--- a/preprocess_thresholds.py
+++ b/preprocess_thresholds.py
@@ -5,122 +5,149 @@ import ancillary_data as anc
import utils
from numpy.lib.stride_tricks import sliding_window_view
+from typing import Dict
_dtr = np.pi/180
_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):
+def prepare_11_12um_thresholds(thresholds: np.ndarray,
+ dim1: int) -> Dict:
- 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'])
- if 'bt1' in list(thresholds['11-12um_Cirrus_Test']):
- bt1 = np.full(data.M01.shape, thresholds['11-12um_Cirrus_Test']['bt1'])
+ coeff_values = np.empty((dim1, 2))
+ coeff_values[:, 0] = np.full(dim1, thresholds['coeffs'][0])
+ coeff_values[:, 1] = np.full(dim1, thresholds['coeffs'][1])
+ cmult_values = np.full(dim1, thresholds['cmult'])
+ adj_values = np.full(dim1, thresholds['adj'])
+ if 'bt1' in list(thresholds):
+ bt1 = np.full(dim1, thresholds['bt1'])
else:
- bt1 = np.full(data.M01.shape, -999)
- if 'lat' in list(thresholds['11-12um_Cirrus_Test']):
- lat = np.full(data.M01.shape, thresholds['11-12um_Cirrus_Test']['lat'])
+ bt1 = np.full(dim1, -999)
+ if 'lat' in list(thresholds):
+ lat = np.full(dim1, thresholds['lat'])
else:
- lat = np.full(data.M01.shape, -999)
- 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),
- 'bt1': (['number_of_lines', 'number_of_pixels'], bt1),
- 'lat': (['number_of_lines', 'number_of_pixels'], lat),
+ lat = np.full(dim1, -999)
+ thr_dict = {'coeffs': coeff_values,
+ 'cmult': cmult_values,
+ 'adj': adj_values,
+ 'bt1': bt1,
+ 'lat': lat,
}
- return xr.Dataset(data_vars=thr_dict)
-# 3. Others
-# - Land_Night
-# - Polar_Night_Land
-# - Polar_Night_Snow
-# - Day_Snow
-# - Night_Snow
+ return thr_dict
-def preproc(data, thresholds, scene):
- cosvza = np.cos(data.sensor_zenith * _DTR)
- schi = (1/cosvza).where(cosvza > 0, 99.0)
+# function was called preproc
+def thresholds_11_12um(data: xr.Dataset,
+ thresholds: np.ndarray,
+ scene: str,
+ scene_idx: np.ndarray) -> np.ndarray:
+ cosvza = np.cos(data.sensor_zenith[scene_idx].values.ravel() * _DTR)
+ schi = np.full(cosvza.shape, 99.0)
+ schi[cosvza > 0] = 1/cosvza
+ schi = np.array(schi, dtype=np.float32) # this is because the C function expects a float
- schi = schi.values.reshape(np.prod(schi.shape))
- m15 = data.M15.values.reshape(np.prod(data.M15.shape))
+ m15 = data.M15.values[scene_idx].ravel()
+ latitude = data.latitude[scene_idx].values.ravel
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)
+ thr_dict = prepare_11_12um_thresholds(thresholds, m15.shape[0])
- midpt = thr_xr.coeffs[:, :, 0].where((thr < 0.1) | (np.abs(schi-99) < 0.0001), thr)
- locut = midpt + (thr_xr.cmult * midpt)
+ midpt = np.full(m15.shape[0], thr)
+ idx = np.nonzero((thr < 0.1) | (np.abs(schi-99) < 0.0001))
+ midpt[idx] = thr_dict['coeffs'][idx, 0]
+ locut = midpt + (thr_dict['cmult'] * midpt)
if scene in ['Land_Day', 'Land_Day_Coast', 'Land_Day_Desert', 'Land_Day_Desert_Coast',
- 'Ocean_Day', 'Ocean_Night', 'Polar_Ocean_Day', 'Polar_Ocean_Night']:
- hicut = midpt - thr_xr.adj
+ 'Ocean_Day', 'Ocean_Night', 'Polar_Day_Ocean', 'Polar_Night_Ocean']:
+ hicut = midpt - thr_dict['adj']
elif scene in ['Polar_Day_Land', 'Polar_Day_Coast', 'Polar_Day_Desert',
'Polar_Day_Desert_Coast', 'Polar_Day_Snow']:
- hicut = midpt - (thr_xr.adj * midpt)
+ hicut = midpt - (thr_dict['adj'] * midpt)
elif scene in ['Land_Night', 'Polar_Night_Land', 'Polar_Night_Snow', 'Day_Snow', 'Night_Snow']:
_coeffs = {'Land_Night': 0.3, 'Polar_Night_Land': 0.3, 'Polar_Night_Snow': 0.3,
'Day_Snow': 0.0, 'Night_Snow': 0.3}
midpt = midpt - (_coeffs[scene] * locut)
if scene in ['Polar_Night_Land', 'Polar_Night_Snow', 'Night_Snow']:
- hicut = (midpt - (0.2 * locut)).where(data.M15 < thr_xr.bt1, midpt - 1.25)
+ hicut = np.full(m15.shape, midpt - 1.25)
+ idx = np.nonzero(m15 < thr_dict['bt1'])
+ hicut[idx] = midpt[idx] - (0.2 * locut[idx])
elif scene in ['Land_Night']:
- hicut = -0.1 - np.power(90.0 - np.abs(data.latitude)/60, 4) * 1.15
- hicut = hicut.where((data.M15 < thr_xr.bt1) & (data.latitude > thr_xr.lat), 1.25)
+ hicut = np.full(m15.shape, 1.25)
+ idx = np.nonzero((m15 < thr_dict['bt1']) & (latitude > thr_dict['lat']))
+ hicut[idx] = -0.1 - np.power(90.0 - np.abs(latitude[idx])/60, 4) * 1.15
elif scene in ['Day_Snow']:
- hicut = locut - (thr_xr.cmult * locut)
+ hicut = locut - (thr_dict['cmult'] * locut)
else:
print('Scene not recognized\n')
else:
print('Scene not recognized\n')
- 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
+ thr_out = np.dstack((locut, midpt, hicut, np.ones(locut.shape), np.ones(locut.shape)))
+ return np.squeeze(thr_out.T)
-def preproc_sst(data, thresholds):
- m31c = data.M15 - 273.16
- m32c = data.M16 - 273.16
- m31c_m32c = m31c - m32c
- sstc = data.geos_sfct - 273.16
- cosvza = np.cos(data.sensor_zenith*_DTR)
+def thresholds_NIR(data, thresholds, scene, test_name, scene_idx):
- a = thresholds['coeffs']
+ sza = data.solar_zenith.values[scene_idx].ravel()
+ band_n = 2
+ # NOTE: the visud condition in the C code is equivalent to having sza <= 85
+ # For the time being the visud filtering is not implemented
- modsst = 273.16 + a[0] + a[1]*m31c + a[2]*m31c_m32c*sstc + a[3]*m31c_m32c*((1/cosvza) - 1)
- sfcdif = data.geos_sfct - modsst
+ c = np.array(thresholds[scene][test_name]['coeffs'])
+ vzcpow = thresholds['VZA_correction']['vzcpow'][0]
+ refang = data.sunglint_angle.values[scene_idx].ravel()
+ sunglint_thresholds = thresholds['Sun_Glint']
+ sunglint_flag = utils.sunglint_scene(refang, sunglint_thresholds)
+ nir_thresh = thresholds[scene][test_name]
- return sfcdif
+ if test_name == 'NIR_Reflectance_Test':
+ hicut = c[0] + c[1]*sza + c[2]*np.power(sza, 2) + c[3]*np.power(sza, 3)
+ elif test_name == '1.6_2.1um_NIR_Reflectance_Test':
+ hicut = c[0] + c[1]*sza + c[2]*np.power(sza, 2) + c[3]*np.power(sza, 3) + c[4]*np.power(sza, 4)
+ else:
+ pass
+ hicut = (hicut * 0.01) + nir_thresh['adj']
+ hicut = (hicut * nir_thresh['bias'])
+ midpt = hicut + (nir_thresh['midpt_coeff'] * nir_thresh['bias'])
+ locut = midpt + (nir_thresh['locut_coeff'] * nir_thresh['bias'])
+ thr = np.array([locut, midpt, hicut, nir_thresh['thr'][3]*np.ones(refang.shape)])
+
+ cosvza = np.cos(data.sensor_zenith.values[scene_idx]*_DTR).ravel()
+
+ corr_thr = np.zeros((4, refang.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):
+ if len(refang[sunglint_flag == flag]) > 0:
+ dosgref = utils.get_sunglint_thresholds(refang, sunglint_thresholds, band_n, flag, thr)
+ corr_thr[:3, sunglint_flag == flag] = dosgref[:3, sunglint_flag == flag] * \
+ (1./np.power(cosvza[sunglint_flag == flag], vzcpow))
+ corr_thr[3, sunglint_flag == flag] = dosgref[3, sunglint_flag == flag]
+ return corr_thr
-def preproc_nir(data, thresholds, scene):
- sza = data.solar_zenith
+def nir_refl(data, thresholds, scene_name):
+ sza = data.solar_zenith.values
band_n = 2
- # NOTE: the visud condition in the C code is equivalent to having sza <= 85
- # For the time being the visud filtering is not implemented
- a = np.array(thresholds[scene]['NIR_Reflectance_Test']['coeffs'])
+ thr_nir = thresholds[scene_name]['1.6_2.1um_NIR_Reflectance_Test']
+ coeffs = thr_nir['coeffs']
vzcpow = thresholds['VZA_correction']['vzcpow'][0]
+
refang = data.sunglint_angle.values.reshape(np.prod(data.sunglint_angle.shape))
+
sunglint_thresholds = thresholds['Sun_Glint']
sunglint_flag = utils.sunglint_scene(refang, sunglint_thresholds).reshape(refang.shape)
- nir_thresh = thresholds[scene]['NIR_Reflectance_Test']
- hicut0 = a[0] + a[1]*sza + a[2]*np.power(sza, 2) + a[3]*np.power(sza, 3)
- hicut0 = (hicut0 * 0.01) + nir_thresh['adj']
- hicut0 = (hicut0 + nir_thresh['bias']).values.reshape(refang.shape)
- midpt0 = hicut0 + (nir_thresh['midpt_coeff'] * nir_thresh['bias'])
- locut0 = midpt0 + (nir_thresh['locut_coeff'] * nir_thresh['bias'])
- thr = np.array([locut0, midpt0, hicut0, nir_thresh['thr'][3]*np.ones(refang.shape)])
+ hicut0 = np.array(coeffs[0] + coeffs[1]*sza + coeffs[2]*np.power(sza, 2) +
+ coeffs[3]*np.power(sza, 3) + coeffs[4]*np.power(sza, 4))
+ hicut0 = (hicut0*0.01) + thr_nir['adj']
+ hicut0 = (hicut0*thr_nir['bias']).reshape(refang.shape)
+ midpt0 = hicut0 + (thr_nir['midpt_coeff']*thr_nir['bias'])
+ locut0 = midpt0 + (thr_nir['locut_coeff']*thr_nir['bias'])
+ thr = np.array([locut0, midpt0, hicut0, np.ones(refang.shape)])
cosvza = np.cos(data.sensor_zenith*_dtr).values.reshape(refang.shape)
corr_thr = np.zeros((4, refang.shape[0]))
@@ -174,6 +201,22 @@ def preproc_surf_temp(data, thresholds):
return thr_out
+# This function is currently not used
+def preproc_sst(data, thresholds):
+ m31c = data.M15 - 273.16
+ m32c = data.M16 - 273.16
+ m31c_m32c = m31c - m32c
+ sstc = data.geos_sfct - 273.16
+ cosvza = np.cos(data.sensor_zenith*_DTR)
+
+ a = thresholds['coeffs']
+
+ modsst = 273.16 + a[0] + a[1]*m31c + a[2]*m31c_m32c*sstc + a[3]*m31c_m32c*((1/cosvza) - 1)
+ sfcdif = data.geos_sfct - modsst
+
+ return sfcdif
+
+
def var_11um(data, thresholds):
rad = data.M15.values
var = np.zeros((rad.shape[0], rad.shape[1], 9))
@@ -430,46 +473,6 @@ def vis_refl_thresholds(data, thresholds, scene):
return out_thr, out_rad
-def nir_refl(data, thresholds, scene_name):
- sza = data.solar_zenith.values
- band_n = 2
-
- thr_nir = thresholds[scene_name]['1.6_2.1um_NIR_Reflectance_Test']
- coeffs = thr_nir['coeffs']
- vzcpow = thresholds['VZA_correction']['vzcpow'][0]
-
- refang = data.sunglint_angle.values.reshape(np.prod(data.sunglint_angle.shape))
-
- sunglint_thresholds = thresholds['Sun_Glint']
- sunglint_flag = utils.sunglint_scene(refang, sunglint_thresholds).reshape(refang.shape)
-
- hicut0 = np.array(coeffs[0] + coeffs[1]*sza + coeffs[2]*np.power(sza, 2) +
- coeffs[3]*np.power(sza, 3) + coeffs[4]*np.power(sza, 4))
- hicut0 = (hicut0*0.01) + thr_nir['adj']
- hicut0 = (hicut0*thr_nir['bias']).reshape(refang.shape)
- midpt0 = hicut0 + (thr_nir['midpt_coeff']*thr_nir['bias'])
- locut0 = midpt0 + (thr_nir['locut_coeff']*thr_nir['bias'])
- thr = np.array([locut0, midpt0, hicut0, np.ones(refang.shape)])
-
- cosvza = np.cos(data.sensor_zenith*_dtr).values.reshape(refang.shape)
- corr_thr = np.zeros((4, refang.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):
- if len(refang[sunglint_flag == flag]) > 0:
- dosgref = utils.get_sunglint_thresholds(refang, sunglint_thresholds, band_n, flag, thr)
- corr_thr[:3, sunglint_flag == flag] = dosgref[:3, sunglint_flag == flag] * \
- (1./np.power(cosvza[sunglint_flag == flag], vzcpow))
- corr_thr[3, sunglint_flag == flag] = dosgref[3, sunglint_flag == flag]
-
- corr_thr = np.transpose(corr_thr.reshape((4, sza.shape[0], sza.shape[1])), (1, 2, 0))
-
- return corr_thr
-
-
def GEMI_test(data, thresholds, scene_name):
thresh = thresholds[scene_name]['GEMI_Test']
diff --git a/read_data.py b/read_data.py
index 3507048d199e0e5d9a6838f4f604c3c8d96c98de..18e3ea4401835b02a5d8d4be99d8e105ba7fa30d 100644
--- a/read_data.py
+++ b/read_data.py
@@ -157,6 +157,8 @@ def get_data(file_names: Dict[str, str],
viirs_data.M15.values - viirs_data.M13.values)
viirs_data['M13-M16'] = (('number_of_lines', 'number_of_pixels'),
viirs_data.M13.values - viirs_data.M16.values)
+ viirs_data['M07-M05ratio'] = (('number_of_lines', 'number_of_pixels'),
+ viirs_data.M07.values / viirs_data.M05.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 474a2407c5541cf4cdc5d98d998da36fc0bc2291..49c747ef06b03c24fc7057254e12c2f69a5f05bb 100644
--- a/scene.py
+++ b/scene.py
@@ -1,5 +1,8 @@
import numpy as np
-import ruamel_yaml as yml
+try:
+ from ruamel import yaml as yml
+except ImportError:
+ import ruamel_yaml as yml
from glob import glob
@@ -7,7 +10,7 @@ import read_data as rd
import ancillary_data as anc
# lsf: land sea flag
-_scene_list = ['Ocean_Day', 'Ocean_Night', 'Land_Day', 'Land_Night', 'Snow_Day', 'Snow_Night', 'Coast_Day',
+_scene_list = ['Ocean_Day', 'Ocean_Night', 'Land_Day', 'Land_Night', 'Day_Snow', 'Night_Snow', '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', 'Desert']
@@ -343,13 +346,13 @@ def scene_id(scene_flag):
idx = np.nonzero((scene_flag['day'] == 1) &
((scene_flag['ice'] == 1) | (scene_flag['snow'] == 1)) &
(scene_flag['polar'] == 1) & (scene_flag['antarctica'] == 1))
- scene['Snow_Day'][idx] = 1
+ scene['Day_Snow'][idx] = 1
# Snow Night
idx = np.nonzero((scene_flag['night'] == 1) &
((scene_flag['ice'] == 1) | (scene_flag['snow'] == 1)) &
(scene_flag['polar'] == 1) & (scene_flag['antarctica'] == 1))
- scene['Snow_Night'][idx] = 1
+ scene['Night_Snow'][idx] = 1
# Land Day Coast
idx = np.nonzero((scene_flag['land'] == 1) & (scene_flag['day'] == 1) &
diff --git a/tests.py b/tests.py
index 2adb710ff6a92315dce9c3e7944a21ae8508eb37..f4c817a872307578b5cd3dcff28c9bb896a7441c 100644
--- a/tests.py
+++ b/tests.py
@@ -4,6 +4,8 @@ import xarray as xr
from numpy.lib.stride_tricks import sliding_window_view
from typing import Dict
+import functools
+
import utils
import conf
import conf_xr
@@ -19,7 +21,7 @@ importlib.reload(pt)
# ############## GROUP 1 TESTS ############## #
-def test_11um(rad, threshold):
+def test_11um_old(rad, threshold):
radshape = rad.shape
rad = rad.reshape(np.prod(radshape))
@@ -236,11 +238,23 @@ class CloudTests(object):
self.thresholds = thresholds
self.scene_idx = tuple(np.nonzero(data[scene_name] == 1))
+ def run_if_test_exists_for_scene(func):
+ @functools.wraps(func)
+ def wrapper(self, *args, test_name, **kwargs):
+ if test_name not in self.thresholds[self.scene_name]:
+ print('Not running test for this scene')
+ # returns cmin. This could be changed into a keyworded argument for readability
+ return args[-1]
+ return func(self, *args, test_name, **kwargs)
+ return wrapper
+
+ @run_if_test_exists_for_scene
def test_11um(self,
band: str,
- cmin: np.ndarray) -> np.ndarray:
- confidence = np.ones(self.data.M01.shape)
- threshold = self.thresholds[self.scene_name]['11um_Test']
+ cmin: np.ndarray,
+ test_name: str = '11um_Test') -> np.ndarray:
+ confidence = np.ones(self.data[band].shape)
+ threshold = self.thresholds[self.scene_name][test_name]
if threshold['perform'] is True:
print(f'Testing "{self.scene_name}"\n')
@@ -251,16 +265,19 @@ class CloudTests(object):
return cmin
+ @run_if_test_exists_for_scene
def sst_test(self,
- band: str,
- cmin: np.ndarray) -> np.ndarray:
+ band31: str,
+ band32: str,
+ cmin: np.ndarray,
+ test_name: str = 'SST_Test') -> np.ndarray:
confidence = np.ones(self.data.M01.shape)
- threshold = self.thresholds[self.scene_name]['SST_Test']
+ threshold = self.thresholds[self.scene_name][test_name]
if threshold['perform'] is True:
- m31 = self.data.M15.values - 273.16
- bt_diff = self.data['M15-M16'].values
+ m31 = self.data[band31].values - 273.16
+ bt_diff = self.data[band31].values - self.data[band32].values
sst = self.data.geos_sfct.values - 273.16
cosvza = np.cos(self.data.sensor_zenith.values*_DTR)
@@ -276,12 +293,14 @@ class CloudTests(object):
return cmin
+ @run_if_test_exists_for_scene
def bt_diff_86_11um(self,
band: str,
- cmin: np.ndarray) -> np.ndarray:
+ cmin: np.ndarray,
+ test_name: str = '8.6-11um_Test') -> np.ndarray:
confidence = np.ones(self.data.M01.shape)
- threshold = self.thresholds[self.scene_name]['8.6-11um_Test']
+ threshold = self.thresholds[self.scene_name][test_name]
if threshold['perform'] is True:
print(f'Testing "{self.scene_name}"\n')
@@ -292,6 +311,21 @@ class CloudTests(object):
return cmin
+ def run_tests_temp(self):
+ cmin_G1 = np.ones(self.data.M01.shape)
+ cmin_G2 = np.ones(self.data.M01.shape)
+
+ # Group 1
+ cmin_G1 = self.test_11um('M15', cmin_G1, test_name='11um_Test')
+ cmin_G1 = self.sst_test('M15', 'M16', cmin_G1, test_name='SST_Test')
+
+ # Group 2
+ cmin_G2 = self.bt_diff_86_11um('M14-M15', cmin_G2, '8.6-11um_Test')
+
+ cmin = cmin_G1 * cmin_G2
+
+ return cmin
+
def single_threshold_test(self, test_name, band, cmin):
if band == 'bad_data':
@@ -409,6 +443,30 @@ class CloudTests(object):
return cmin
+class ComputeTests(CloudTests):
+
+ def __init__(self,
+ data: xr.Dataset,
+ scene_name: str,
+ thresholds: Dict) -> None:
+ super().__init__(data, scene_name, thresholds)
+
+ def run_tests(self):
+ cmin_G1 = np.ones(self.data.M01.shape)
+ cmin_G2 = np.ones(self.data.M01.shape)
+
+ # Group 1
+ cmin_G1 = self.test_11um('M15', cmin_G1, test_name='11um_Test')
+ cmin_G1 = self.sst_test('M15', 'M16', cmin_G1, test_name='SST_Test')
+
+ # Group 2
+ cmin_G2 = self.bt_diff_86_11um('M14-M15', cmin_G2, test_name='8.6-11um_Test')
+
+ cmin = cmin_G1 * cmin_G2
+
+ return cmin
+
+
def preproc_thresholds(thresholds, data):
thr = np.array(thresholds)
thr_xr = xr.Dataset()
diff --git a/thresholds.mvcm.snpp.v0.0.1.yaml b/thresholds.mvcm.snpp.v0.0.1.yaml
index 0cffd7aefe37013ef480d61d2dc15e1bfd665a22..bd850c9b64de384895265f48a100fda56e9396bd 100644
--- a/thresholds.mvcm.snpp.v0.0.1.yaml
+++ b/thresholds.mvcm.snpp.v0.0.1.yaml
@@ -8,6 +8,7 @@ Land_Day:
coeffs: [3.0, 1.0]
cmult: 0.3
adj: 1.25
+ perform: True
dl_ref3_tpw : -10.00
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]
@@ -27,6 +28,7 @@ Land_Night:
adj: 0
bt1: 270.0 # WILL FIND A BETTER NAME AT SOME POINT, MAYBE
lat: 30.0
+ perform: True
4-12um_BTD_Thin_Cirrus_Test:
thr: [15.0, 10.0, 5.00, 1.0, 1.0]
CO2_High_Clouds_Test: [222.0, 224.0, 226.0, 1.0, 1.0]
@@ -58,6 +60,7 @@ Land_Day_Coast:
coeffs: [3.0, 1.0]
cmult: 0.3
adj: 1.25
+ perform: True
dl_ref3_tpw_t2 : -10.00
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]
@@ -75,7 +78,8 @@ Land_Day_Desert:
coeffs: [3.5, 1.0]
cmult: 0.3
adj: 1.25
- 11-4um_Oceanic_Stratus_Test: [-28.0, -26.0, -24.0, -2.0, 0.0, 2.0, 1.0, 1.0] # this preplaces lds11_4hi and
+ perform: True
+ 11-4um_Oceanic_Stratus_Test: [-28.0, -26.0, -24.0, -2.0, 0.0, 2.0, 1.0, 1.0] # this replaces lds11_4hi and
# lds11_4lo. The values are
# sorted left to right so that
# they follow the scheme:
@@ -106,6 +110,7 @@ Land_Day_Desert_Coast:
coeffs: [3.5, 1.0]
cmult: 0.3
adj: 1.25
+ perform: True
11-4um_Oceanic_Stratus_Test: [-23.0, -21.0, -19.0, -2.0, 0.0, 2.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]
@@ -127,38 +132,43 @@ Ocean_Day:
thr: [3.000, 2.500, 2.000, 1.0, 1.0]
coeffs: [1.886, 0.938, 0.128, 1.094]
perform: True
+ 8.6-11um_Test:
+ thr: [-0.50, -1.00, -1.50, 1.0, 1.0]
+ perform: True
11-12um_Cirrus_Test:
coeffs: [3.0, 1.0]
cmult: 0.3
adj: 1.25
- 8.6-11um_Test:
- thr: [-0.50, -1.00, -1.50, 1.0, 1.0]
perform: True
11-4um_Oceanic_Stratus_Test: [-11.0, -9.0, -7.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]
- 1.6_2.1um_NIR_Reflectance_Test:
- coeffs: [1.01339303876915, -0.00277128813739, 0.00027804834484, -0.00001549681141, 0.00000016623006, 0.00000000000000]
- locut_coeff: 0.0100
- midpt_coeff: 0.0070
- adj: 0.0125
- bias: 0.98 # this corresponds to do_b6bias_adj
- perform_test: True
NIR_Reflectance_Test:
thr: [0.062, 0.043, 0.029, 1.0, 1.0]
coeffs: [1.7291, 0.0715, -0.0026, 0.000025889]
adj: 0.0050
bias: 0.96
+ perform: True
midpt_coeff: 0.0200
locut_coeff: 0.0100
+ Vis/NIR_Ratio_Test:
+ thr: [0.837, 0.889, 0.934, 1.001, 1.154, 1.205, 1.0, 1.0] # This replace dovrathi and dovratlo
+ perform: True # The values are sorted left to right
+ # so that they follow the scheme:
+ # Hi-Mid-Lo--Lo-Mid-Hi
+ # Last two values are power and switch
+ # to turn the test on/off
+ 1.6_2.1um_NIR_Reflectance_Test:
+ thr: [0.0, 0.0, 0.0, 1.0, 1.0] # This is a placeholder to force consistency with NIR_Reflectance_Test.
+ # In both cases the only value used is thr[3] AFAIK
+ coeffs: [1.01339303876915, -0.00277128813739, 0.00027804834484, -0.00001549681141, 0.00000016623006, 0.00000000000000]
+ bias: 0.98 # this corresponds to do_b6bias_adj
+ adj: 0.0125
+ locut_coeff: 0.0100
+ midpt_coeff: 0.0070
+ perform: True
+ 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]
# vnir_ratio_hi : [1.001, 1.154, 1.205, 1.0, 1.0]
# vnir_ratio_lo : [0.934, 0.889, 0.837, 1.0]
- Vis/NIR_Ratio_Test: [0.837, 0.889, 0.934, 1.001, 1.154, 1.205, 1.0, 1.0] # This replace dovrathi and dovratlo
- # The values are sorted left to right
- # so that they follow the scheme:
- # Hi-Mid-Lo--Lo-Mid-Hi
- # Last two values are power and switch
- # to turn the test on/off
1.38um_High_Cloud_Test:
coeffs: [-0.2419, 0.0455, -0.0024, 0.000059029, -0.00000069964, 0.000000003204]
adj: 0.0010
@@ -182,6 +192,7 @@ Ocean_Night:
coeffs: [3.0, 1.0]
cmult: 0.3
adj: 1.0
+ perform: True
11-4um_Oceanic_Stratus_Test: [1.25, 1.00, 0.25, 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]
@@ -200,6 +211,7 @@ Polar_Day_Land:
coeffs: [3.0, 1.0]
cmult: 0.3
adj: 0.3
+ perform: True
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]
@@ -215,6 +227,7 @@ Polar_Night_Land:
cmult: 0
adj: 0 # I NEED TO WORK ON THIS
bt1: 270.0
+ perform: True
pnlbt2 : 270.0
Surface_Temperature_Test_df1: [0.0, 1.0]
Surface_Temperature_Test_df2: [-0.5, 1.0]
@@ -246,6 +259,7 @@ Polar_Day_Coast:
coeffs: [3.0, 1.0]
cmult: 0 # I NEED TO WORK ON THIS
adj: 0 # I NEED TO WORK ON THIS
+ perform: True
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]
@@ -259,6 +273,7 @@ Polar_Day_Desert:
coeffs: [3.5, 1.0]
cmult: 0 # I NEED TO WORK ON THIS
adj: 0 # I NEED TO WORK ON THIS
+ perform: True
11-4um_Oceanic_Stratus_Test: [-22.0, -20.0, -18.0, -2.0, 0.0, 2.0, 1.0, 1.0]
Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0]
Visible_Reflectance_Test:
@@ -277,6 +292,7 @@ Polar_Day_Desert_Coast:
coeffs: [3.5, 1.0]
cmult: 0 # I NEED TO WORK ON THIS
adj: 0 # I NEED TO WORK ON THIS
+ perform: True
11-4um_Oceanic_Stratus_Test: [-23.0, -21.0, -19.0, -2.0, 0.0, 2.0, 1.0, 1.0]
Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0]
Visible_Reflectance_Test:
@@ -291,6 +307,7 @@ Polar_Day_Snow:
coeffs: [3.0, 1.0]
cmult: 0 # I NEED TO WORK ON THIS
adj: 0 # I NEED TO WORK ON THIS
+ perform: True
dps_ref3_tpw : 0.75
dpsbt1 : 230.0
Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0]
@@ -310,6 +327,7 @@ Polar_Night_Snow:
cmult: 0
adj: 0 # I NEED TO WORK ON THIS
bt1: 270.0
+ perform: True
4-12um_BTD_Thin_Cirrus_Test:
low: [4.50, 4.00, 3.50, 1.0] # pn_4_12l
mid1: [4.00, 2.00, 0.500, 1.0] # pn_4_12m1
@@ -362,6 +380,7 @@ Polar_Day_Ocean:
coeffs: [3.0, 1.0]
cmult: 0.3
adj: 1.25
+ perform: True
NIR_Reflectance_Test:
thr: [0.062, 0.043, 0.029, 1.0, 1.0]
coeffs: [1.7291, 0.0715, -0.0026, 0.000025889]
@@ -369,16 +388,20 @@ Polar_Day_Ocean:
bias: 0.96
midpt_coeff: 0.0200
locut_coeff: 0.0100
+ perform: True
+ Vis/NIR_Ratio_Test:
+ thr: [0.837, 0.889, 0.934, 1.001, 1.154, 1.205, 1.0, 1.0]
+ perform: True
1.6_2.1um_NIR_Reflectance_Test:
+ thr: [0.0, 0.0, 0.0, 1.0, 1.0]
coeffs: [1.01339303876915, -0.00277128813739, 0.00027804834484, -0.00001549681141, 0.00000016623006, 0.00000000000000]
+ bias: 0.98
+ adj: 0.0125
locut_coeff: 0.0100
midpt_coeff: 0.0070
- adj: 0.0125
- bias: 0.98
- perform_test: True
+ perform: True
11-4um_Oceanic_Stratus_Test: [-11.0, -9.0, -7.0, 1.0, 1.0]
Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0]
- Vis/NIR_Ratio_Test: [0.837, 0.889, 0.934, 1.001, 1.154, 1.205, 1.0, 1.0]
pdo_b26pfm : 1.0
pdo_b2bias_adj : 0.96
pdo_b7bias_adj : 0.97
@@ -403,6 +426,7 @@ Polar_Night_Ocean:
coeffs: [3.0, 1.0]
cmult: 0.3
adj: 1.0
+ perform: True
pnobt1 : 280.0
11-4um_Oceanic_Stratus_Test: [1.25, 1.00, 0.25, 1.0, 1.0]
Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0]
@@ -422,6 +446,7 @@ Day_Snow:
coeffs: [3.0, 1.0]
cmult: 0.3
adj: 0 # I NEED TO WORK ON THIS
+ perform: True
ds_ref3_tpw : 0.75
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]
@@ -437,6 +462,7 @@ Night_Snow:
cmult: 0
adj: 0 # I NEED TO WORK ON THIS
bt1: 270.0
+ perform: True
nsbt2 : 270.0
nsbt3 : 230.0
ns11_4lo : [0.70, 0.60, 0.50, 1.0, 1.0]