diff --git a/main.py b/main.py
index 32e01b3b63ba372dcd18ccb283775688b448a14c..10b5f1f34571df28759321e47ff3dab12ce0a523 100644
--- a/main.py
+++ b/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:
diff --git a/preprocess_thresholds.py b/preprocess_thresholds.py
index fc7fc7fef76a9a928c30d6d20846d555c0f2f771..daf69b5b113944a046e22d2c9c0ef7a1c0ec8b4f 100644
--- a/preprocess_thresholds.py
+++ b/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
diff --git a/tests.py b/tests.py
index 8e79c3612125d50198ea085691704e5a4d0e65be..2debb8b829696f2342b3126a8ac278f157f426ea 100644
--- a/tests.py
+++ b/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
diff --git a/thresholds.mvcm.snpp.v0.0.1.yaml b/thresholds.mvcm.snpp.v0.0.1.yaml
index 1f495c88f6ab2f20361d0148d9c4d89654c6a216..4debff3cecad2fe10b8770bf7e859f9283249e00 100644
--- a/thresholds.mvcm.snpp.v0.0.1.yaml
+++ b/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]