diff --git a/preprocess_thresholds.py b/preprocess_thresholds.py index 38b61a04dfc1f42e0cf957e370f58647ab5094f9..fd101d173457b1d46e4028207eb6bf676a33887c 100644 --- a/preprocess_thresholds.py +++ b/preprocess_thresholds.py @@ -429,6 +429,107 @@ 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'] + + gemi_thr = np.ones((data.M01.shape[0], data.M01.shape[1], 5)) + + idx = np.nonzero(data.ndvi < 0.1) + gemi_thr[idx[0], idx[1], :3] = thresh['gemi0'][:3] + idx = np.nonzero((data.ndvi >= 0.1) & (data.ndvi < 0.2)) + gemi_thr[idx[0], idx[1], :3] = thresh['gemi1'][:3] + idx = np.nonzero((data.ndvi >= 0.2) & (data.ndvi < 0.3)) + gemi_thr[idx[0], idx[1], :3] = thresh['gemi2'][:3] + + thr_out = xr.DataArray(data=np.dstack((gemi_thr[:, :, 0], gemi_thr[:, :, 1], gemi_thr[:, :, 2], + np.ones(gemi_thr[:, :, 0].shape), + np.ones(gemi_thr[:, :, 0].shape))), + dims=('number_of_lines', 'number_of_pixels', 'z')) + return thr_out + + +def bt11_4um_preproc(data, thresholds, scene_name): + thresh = thresholds[scene_name]['11-4um_BT_Difference_Test'] + c = thresh['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_out = np.vstack([hicut0, midpt0, locut0, locut1, midpt1, hicut1, + np.ones(hicut0.shape), np.ones(hicut0.shape)]) + return thr_out + + +def test_1_38um_preproc(data, thresholds, scene_name): + sza = data.solar_zenith.values + vza = data.sensor_zenith.values + thresh = thresholds[scene_name]['1.38um_High_Cloud_Test'] + c = thresh['coeffs'] + vzcpow = thresholds['VZA_correction']['vzcpow'][1] + + hicut0 = c[0] + c[1]*sza + c[2]*np.power(sza, 2) + c[3]*np.power(sza, 3) + \ + c[4]*np.power(sza, 4) + c[5]*np.power(sza, 5) + hicut0 = hicut0*0.01 + (np.maximum((sza/90.)*thresh['szafac']*thresh['adj'], thresh['adj'])) + midpt0 = hicut0 + 0.001 + locut0 = midpt0 + 0.001 + + cosvza = np.cos(vza*_dtr) + + locut0 = locut0 * (1/np.power(cosvza, vzcpow)) + midpt0 = midpt0 * (1/np.power(cosvza, vzcpow)) + hicut0 = hicut0 * (1/np.power(cosvza, vzcpow)) + + out_thr = xr.DataArray(data=np.dstack((locut0, midpt0, hicut0, np.ones(data.ndvi.shape), + np.ones(data.ndvi.shape))), + dims=('number_of_lines', 'number_of_pixels', 'z')) + return out_thr + # NOTE: 11-12um Cirrus Test # hicut is computed in different ways depending on the scene # 1. midpt - adj diff --git a/tests.py b/tests.py index 8e39ff93ff3c3da3038659847bc5319ff630f3c9..82df1fd97d2f0966ae220033c528c8acedc143eb 100644 --- a/tests.py +++ b/tests.py @@ -254,6 +254,16 @@ class CloudTests: self.scene_name = scene_name self.thresholds = thresholds + def test_11um(self, band): + confidence = np.ones(self.data.M01.shape) + if self.thresholds[self.scene_name]['11um_Test']['perform'] is True: + + confidence = conf.conf_test(data, band) + + cmin = np.fmin(cmin, confidence) + + return cmin + def single_threshold_test(self, test_name, band, cmin): if band == 'bad_data': @@ -292,16 +302,24 @@ class CloudTests: thr_xr['threshold'], self.data['M128'] = pt.vis_refl_thresholds(self.data, self.thresholds, self.scene_name) + elif test_name == '1.6_2.1um_NIR_Reflectance_Test': + corr_thr = pt.nir_refl(self.data, self.thresholds, self.scene_name) + thr_xr['threshold'] = (('number_of_lines', 'number_of_pixels', 'z'), corr_thr) + thr = np.ones((5,)) elif test_name == '4-12um_BTD_Thin_Cirrus_Test': thr_xr['threshold'] = pt.get_pn_thresholds(self.data, self.thresholds, self.scene_name, '4-12um_BTD_Thin_Cirrus_Test') thr = np.ones((5,)) - elif test_name == 'Vis/NIR_Ratio_Test': - thr_no_sunglint = np.array([self.thresholds[self.scene_name][test_name][i] for i in range(8)]) - thr_sunglint = np.array([self.thresholds['Sun_Glint']['snglnt'][i] for i in range(8)]) + elif test_name == 'GEMI_Test': + thr_xr['threshold'] = pt.GEMI_test(self.data, self.thresholds, self.scene_name) + thr = np.ones((5,)) + elif (test_name == '1.38um_High_Cloud_Test' and self.scene_name in ['Ocean_Day', 'Polar_Ocean_Day']): + thr_xr['threshold'] = pt.test_1_38um_preproc(self.data, self.thresholds, self.scene_name) + thr = np.ones((5,)) 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 test_name == 'SST_Test': @@ -313,19 +331,27 @@ class CloudTests: data['11um_var'] = data.M15 data['11um_var'].values[var != 9] = np.nan - if test_name == 'NIR_Reflectance_Test': - pass + if thr[4] == 1: + print('test running...') + confidence = conf_xr.conf_test(data, band) - # if thr[4] == 1: - if (test_name == '11-4um_Oceanic_Stratus_Test' and - self.scene_name in ['Land_Day_Desert', 'Land_Day_Desert_Coast', 'Polar_Day_Desert', - 'Polar_Day_Desert_Coast']): + cmin = np.fmin(cmin, confidence) + + return cmin + + def double_threshold_test(self, test_name, band, cmin): + data = self.data + if test_name == '11-4um_BT_Difference_Test': + thr = pt.bt11_4um_preproc(self.data, self.thresholds, self.scene_name) print('test running...') - confidence = conf.conf_test_dble(data['M15-M16'].values, thr) + confidence = conf.conf_test_dble(data['M15-M13'].values, thr) confidence = confidence.reshape(data.M01.shape) - elif test_name == 'Vis/NIR_Ratio_Test': - vrat = data.M07.values/data.M05.values + + if test_name == 'Vis/NIR_Ratio_Test': print('test running...') + thr_no_sunglint = np.array([self.thresholds[self.scene_name][test_name][i] for i in range(8)]) + thr_sunglint = np.array([self.thresholds['Sun_Glint']['snglnt'][i] for i in range(8)]) + vrat = data.M07.values/data.M05.values _dtr = np.pi/180.0 sza = data.sensor_zenith.values raz = data.relative_azimuth.values @@ -340,9 +366,15 @@ class CloudTests: confidence[idx] = conf.conf_test_dble(vrat[idx], thr_sunglint) confidence = confidence.reshape(data.M01.shape) - else: + + if (test_name == '11-4um_Oceanic_Stratus_Test' and + self.scene_name in ['Land_Day_Desert', 'Land_Day_Desert_Coast', 'Polar_Day_Desert', + 'Polar_Day_Desert_Coast']): + thr = np.array([self.thresholds[self.scene_name][test_name][i] for i in range(8)]) + print('test running...') - confidence = conf_xr.conf_test(data, band) + confidence = conf.conf_test_dble(data['M15-M16'].values, thr) + confidence = confidence.reshape(data.M01.shape) cmin = np.fmin(cmin, confidence) diff --git a/thresholds.mvcm.snpp.v0.0.1.yaml b/thresholds.mvcm.snpp.v0.0.1.yaml index fd7af279f4fce4a20b5ef0bb6e44b538cb0dc6d5..4ec8be03dac71f123c3edd71b5bb70804c016b3c 100644 --- a/thresholds.mvcm.snpp.v0.0.1.yaml +++ b/thresholds.mvcm.snpp.v0.0.1.yaml @@ -91,9 +91,10 @@ Land_Day_Desert: 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] + GEMI_Test: + gemi0: [0.085, 0.095, 0.115, 1.00, 1.0] + gemi1: [0.145, 0.170, 0.220, 1.00] + gemi2: [0.310, 0.335, 0.360, 1.00] lds_ref3_tpw : -10.00 ldsbt1 : 320.0 # lds11_12lcmult : 0.3 @@ -124,14 +125,18 @@ Ocean_Day: cmult: 0.3 adj: 1.25 11-4um_Oceanic_Stratus_Test: [-11.0, -9.0, -7.0, 1.0, 1.0] - 11um_Test: [267.0, 270.0, 273.0, 1.0, 1.0] + 11um_Test: + thr: [267.0, 270.0, 273.0, 1.0, 1.0] + 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] - do_b7coeffs : [1.01339303876915, -0.00277128813739, 0.00027804834484, -0.00001549681141, 0.00000016623006, 0.00000000000000] - do_b7adj : 0.0125 - do_b7mid : 0.0070 - do_b7lo : 0.0100 - do_b7pfm : 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] @@ -151,10 +156,11 @@ Ocean_Day: SST_Test: thr: [3.000, 2.500, 2.000, 1.0, 1.0] coeffs: [1.886, 0.938, 0.128, 1.094] - do_b26coeffs : [-0.2419, 0.0455, -0.0024, 0.000059029, -0.00000069964, 0.000000003204] - do_b26adj : 0.0010 - do_b26_szafac : 4.0 - do_b26pfm : 1.0 + 1.38um_High_Cloud_Test: + coeffs: [-0.2419, 0.0455, -0.0024, 0.000059029, -0.00000069964, 0.000000003204] + adj: 0.0010 + szafac: 4.0 + perform: True do_b6bias_adj : 0.98 do_b7bias_adj : 0.97 @@ -163,7 +169,6 @@ Ocean_Night: coeffs: [3.0, 1.0] cmult: 0.3 adj: 1.0 - no_intadj : 0.0 11-4um_Oceanic_Stratus_Test: [1.25, 1.00, 0.25, 1.0, 1.0] 11um_Test: [267.0, 270.0, 273.0, 2.0, 1.0] CO2_High_Clouds_Test: [222.0, 224.0, 226.0, 1.0, 1.0] @@ -174,11 +179,13 @@ Ocean_Night: SST_Test: thr: [3.000, 2.500, 2.000, 1.0, 1.0] coeffs: [1.8860, 0.9380, 0.1280, 1.094] - no11_4coef : [-0.7093, 0.1128, -0.1567] - no11_4hiad : [-1.5, 1.5] - no11_4mdad : [-2.5, 2.5] - no11_4load : [-4.0, 4.0] - no11_4_pfm : 1.0 + 11-4um_BT_Difference_Test: + coeffs: [-0.7093, 0.1128, -0.1567] + corr: 0.0 # no_intadj + locut_coeff: [-4.0, 4.0] # no11_4load + midpt_coeff: [-2.5, 2.5] # no11_4mdad + hicut_coeff: [-1.5, 1.5] # no11_4hiad + perform: True # no11_4_pfm Polar_Day_Land: 11-12um_Cirrus_Test: @@ -250,9 +257,10 @@ Polar_Day_Desert: 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] - pdsgemi2 : [0.310, 0.335, 0.360, 1.00] + GEMI_Test: + gemi0: [0.085, 0.110, 0.135, 1.00, 1.0] + gemi1: [0.170, 0.220, 0.270, 1.00] + gemi2: [0.310, 0.335, 0.360, 1.00] pds_ref3_tpw : -10.00 pdsbt1 : 320.0 @@ -343,6 +351,13 @@ Polar_Ocean_Day: bias: 0.96 midpt_coeff: 0.0200 locut_coeff: 0.0100 + 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 + perform_test: True 11-4um_Oceanic_Stratus_Test: [-11.0, -9.0, -7.0, 1.0, 1.0] 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] @@ -351,11 +366,14 @@ Polar_Ocean_Day: SST_Test: thr: [3.000, 2.500, 2.000, 1.0, 1.0] coeffs: [1.886, 0.938, 0.128, 1.094] - pdo_b7pfm : 1.0 pdo_b26pfm : 1.0 pdo_b2bias_adj : 0.96 - pdo_b6bias_adj : 0.98 pdo_b7bias_adj : 0.97 + 1.38um_High_Cloud_Test: + coeffs: [-0.2419, 0.0455, -0.0024, 0.000059029, -0.00000069964, 0.000000003204] + adj: 0.0010 + szafac: 4.0 + perform: True Polar_Ocean_Night: 11-12um_Cirrus_Test: @@ -363,9 +381,7 @@ Polar_Ocean_Night: cmult: 0.3 adj: 1.0 pnobt1 : 280.0 - pno_intadj : 0.0 11-4um_Oceanic_Stratus_Test: [1.25, 1.00, 0.25, 1.0, 1.0] - pno_11_4_pfm : 1.0 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] @@ -374,10 +390,14 @@ Polar_Ocean_Night: SST_Test: thr: [3.000, 2.500, 2.000, 1.0, 1.0] coeffs: [1.8860, 0.9380, 0.1280, 1.094] - pno11_4hiad : [-1.5, 1.5] - pno11_4mdad : [-2.5, 2.5] - pno11_4load : [-4.0, 4.0] - pno11_4coef : [-0.7093, 0.1128, -0.1567] + 11-4um_BT_Difference_Test: + coeffs: [-0.7093, 0.1128, -0.1567] + corr: 0.0 # pno_intadj + locut_coeff: [-4.0, 4.0] # pno11_4load + midpt_coeff: [-2.5, 2.5] # pno11_4mdad + hicut_coeff: [-1.5, 1.5] # pno11_4hiad + perform: True # pno11_4_pfm + Day_Snow: 11-12um_Cirrus_Test: