diff --git a/modules/icing/pirep_goes.py b/modules/icing/pirep_goes.py
index 7443799958869c605ceaad35ca66d44006a57ecb..d7d6ab5f9ebb60ff2bced7b41a1569c35c5bde66 100644
--- a/modules/icing/pirep_goes.py
+++ b/modules/icing/pirep_goes.py
@@ -1069,6 +1069,141 @@ def run_qc(filename, filename_l1b, day_night='ANY', pass_thresh_frac=0.20, icing
return len(icing_alt), len(mask), len(keep_idxs)
+# def apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_opd, cld_mask, bt_11um, solzen, satzen, day_night='ANY'):
+#
+# if day_night == 'DAY':
+# opd_thick_threshold = 20
+# opd_thin_threshold = 1
+# elif day_night == 'NIGHT' or day_night == 'ANY':
+# opd_thick_threshold = 2
+# opd_thin_threshold = 0.1
+#
+# closeness_top = 100.0 # meters
+# closeness_bot = 100.0
+# max_depth = 3000.0
+#
+# num_obs = len(icing_alt)
+# cld_mask = cld_mask.reshape((num_obs, -1))
+# cld_top_hgt = cld_top_hgt.reshape((num_obs, -1))
+# cld_geo_dz = cld_geo_dz.reshape((num_obs, -1))
+# bt_11um = bt_11um.reshape((num_obs, -1))
+#
+# mask = []
+# idxs = []
+# num_tested = []
+# for i in range(num_obs):
+# if not check_oblique(satzen[i,]):
+# continue
+# if day_night == 'NIGHT' and not is_night(solzen[i,]):
+# continue
+# elif day_night == 'DAY' and not is_day(solzen[i,]):
+# continue
+# elif day_night == 'ANY':
+# pass
+# # if not (is_day(solzen[i,]) or is_night(solzen[i,])):
+# # continue
+#
+# keep_0 = np.logical_or(cld_mask[i,] == 2, cld_mask[i,] == 3) # cloudy
+# keep_1 = np.invert(np.isnan(cld_top_hgt[i,]))
+# keep_2 = np.invert(np.isnan(bt_11um[i,]))
+# keep_3 = np.invert(np.isnan(cld_geo_dz[i,]))
+# keep = keep_0 & keep_1 & keep_2 & keep_3
+# num_keep = np.sum(keep)
+# if num_keep == 0:
+# continue
+#
+# # Test 1
+# keep = np.where(keep, (cld_top_hgt[i,] + closeness_top) > icing_alt[i], False)
+# # keep = np.where(keep, (cld_top_hgt[i,] - cld_geo_dz[i,] - closeness_bot) < icing_alt[i], False)
+# keep = np.where(keep, (cld_top_hgt[i,] - max_depth) < icing_alt[i], False)
+#
+# # # Test2
+# # keep = np.where(keep,
+# # np.invert((cld_phase[i,] == 4) &
+# # np.logical_and(cld_top_hgt[i,]+closeness > icing_alt[i], cld_top_hgt[i,]-closeness < icing_alt[i])),
+# # False)
+# #
+# # # Test4
+# # keep = np.where(keep, np.invert((cld_phase[i,] == 4) & (cld_opd[i,] < opd_thin_threshold)), False)
+# #
+# # # Test5 and Test6
+# # keep = np.where(keep, np.logical_and(bt_11um[i,] > 228.0, bt_11um[i,] < 270.0), False)
+# #
+# # # Test3
+# # keep = np.where(keep, (cld_opd[i,] >= opd_thick_threshold) & (cld_phase[i,] == 4), False)
+#
+# mask.append(keep)
+# idxs.append(i)
+# num_tested.append(num_keep)
+#
+# return mask, idxs, num_tested
+
+
+# def apply_qc_no_icing_pireps(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_opd, cld_mask, bt_11um, solzen, satzen, day_night='ANY'):
+#
+# if day_night == 'DAY':
+# opd_thick_threshold = 20
+# opd_thin_threshold = 1
+# elif day_night == 'NIGHT' or day_night == 'ANY':
+# opd_thick_threshold = 2
+# opd_thin_threshold = 0.1
+#
+# closeness_top = 100.0 # meters
+# closeness_bot = 200.0
+# max_depth = 3000.0
+#
+# num_obs = len(icing_alt)
+# cld_mask = cld_mask.reshape((num_obs, -1))
+# cld_top_hgt = cld_top_hgt.reshape((num_obs, -1))
+# cld_geo_dz = cld_geo_dz.reshape((num_obs, -1))
+# bt_11um = bt_11um.reshape((num_obs, -1))
+#
+# mask = []
+# idxs = []
+# num_tested = []
+# for i in range(num_obs):
+# if not check_oblique(satzen[i,]):
+# continue
+# if day_night == 'NIGHT' and not is_night(solzen[i,]):
+# continue
+# elif day_night == 'DAY' and not is_day(solzen[i,]):
+# continue
+# elif day_night == 'ANY':
+# pass
+# # if not (is_day(solzen[i,]) or is_night(solzen[i,])):
+# # continue
+#
+# keep_0 = np.logical_or(cld_mask[i,] == 2, cld_mask[i,] == 3) # cloudy
+# keep_1 = np.invert(np.isnan(cld_top_hgt[i,]))
+# keep_2 = np.invert(np.isnan(bt_11um[i,]))
+# keep_3 = np.invert(np.isnan(cld_geo_dz[i,]))
+# keep = keep_0 & keep_1 & keep_2 & keep_3
+# num_keep = np.sum(keep)
+# if num_keep == 0:
+# continue
+#
+# keep = np.where(keep, (cld_top_hgt[i,] + closeness_top) > icing_alt[i], False)
+# # keep = np.where(keep, (cld_top_hgt[i,] - cld_geo_dz[i,] - closeness_bot) < icing_alt[i], False)
+# keep = np.where(keep, (cld_top_hgt[i,] - max_depth) < icing_alt[i], False)
+#
+# # keep = np.where(keep, np.logical_and(bt_11um[i,] > 228.0, bt_11um[i,] < 273.0), False)
+#
+# # Test6
+# # keep = np.where(keep, (bt_11um[i,] < 228.0), False)
+#
+# # Test1
+# # keep = np.where(keep, cld_top_hgt[i,] > icing_alt[i], False)
+#
+# # Test3
+# # keep = np.where(keep, np.invert((cld_opd[i,] >= opd_thick_threshold) & (cld_phase[i,] == 4) & (cld_top_hgt[i,] > icing_alt[i])), False)
+#
+# mask.append(keep)
+# idxs.append(i)
+# num_tested.append(num_keep)
+#
+# return mask, idxs, num_tested
+
+
def apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_opd, cld_mask, bt_11um, solzen, satzen, day_night='ANY'):
if day_night == 'DAY':
@@ -1078,8 +1213,7 @@ def apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_opd
opd_thick_threshold = 2
opd_thin_threshold = 0.1
- closeness_top = 100.0 # meters
- closeness_bot = 100.0
+ closeness_top = 200.0 # meters
max_depth = 3000.0
num_obs = len(icing_alt)
@@ -1103,6 +1237,10 @@ def apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_opd
# if not (is_day(solzen[i,]) or is_night(solzen[i,])):
# continue
+ # if not (500.0 < icing_alt[i] < 3000.0):
+ if not (icing_alt[i] < 3000.0):
+ continue
+
keep_0 = np.logical_or(cld_mask[i,] == 2, cld_mask[i,] == 3) # cloudy
keep_1 = np.invert(np.isnan(cld_top_hgt[i,]))
keep_2 = np.invert(np.isnan(bt_11um[i,]))
@@ -1114,23 +1252,9 @@ def apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_opd
# Test 1
keep = np.where(keep, (cld_top_hgt[i,] + closeness_top) > icing_alt[i], False)
- # keep = np.where(keep, (cld_top_hgt[i,] - cld_geo_dz[i,] - closeness_bot) < icing_alt[i], False)
keep = np.where(keep, (cld_top_hgt[i,] - max_depth) < icing_alt[i], False)
- # # Test2
- # keep = np.where(keep,
- # np.invert((cld_phase[i,] == 4) &
- # np.logical_and(cld_top_hgt[i,]+closeness > icing_alt[i], cld_top_hgt[i,]-closeness < icing_alt[i])),
- # False)
- #
- # # Test4
- # keep = np.where(keep, np.invert((cld_phase[i,] == 4) & (cld_opd[i,] < opd_thin_threshold)), False)
- #
- # # Test5 and Test6
- # keep = np.where(keep, np.logical_and(bt_11um[i,] > 228.0, bt_11um[i,] < 270.0), False)
- #
- # # Test3
- # keep = np.where(keep, (cld_opd[i,] >= opd_thick_threshold) & (cld_phase[i,] == 4), False)
+ keep = np.where(keep, np.logical_and(bt_11um[i,] > 228.0, bt_11um[i,] < 270.0), False)
mask.append(keep)
idxs.append(i)
@@ -1148,8 +1272,7 @@ def apply_qc_no_icing_pireps(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_
opd_thick_threshold = 2
opd_thin_threshold = 0.1
- closeness_top = 100.0 # meters
- closeness_bot = 200.0
+ closeness_top = 200.0 # meters
max_depth = 3000.0
num_obs = len(icing_alt)
@@ -1173,6 +1296,10 @@ def apply_qc_no_icing_pireps(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_
# if not (is_day(solzen[i,]) or is_night(solzen[i,])):
# continue
+ # if not (500.0 < icing_alt[i] < 3000.0):
+ if not (icing_alt[i] < 3000.0):
+ continue
+
keep_0 = np.logical_or(cld_mask[i,] == 2, cld_mask[i,] == 3) # cloudy
keep_1 = np.invert(np.isnan(cld_top_hgt[i,]))
keep_2 = np.invert(np.isnan(bt_11um[i,]))
@@ -1183,19 +1310,9 @@ def apply_qc_no_icing_pireps(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_
continue
keep = np.where(keep, (cld_top_hgt[i,] + closeness_top) > icing_alt[i], False)
- # keep = np.where(keep, (cld_top_hgt[i,] - cld_geo_dz[i,] - closeness_bot) < icing_alt[i], False)
keep = np.where(keep, (cld_top_hgt[i,] - max_depth) < icing_alt[i], False)
- # keep = np.where(keep, np.logical_and(bt_11um[i,] > 228.0, bt_11um[i,] < 273.0), False)
-
- # Test6
- # keep = np.where(keep, (bt_11um[i,] < 228.0), False)
-
- # Test1
- # keep = np.where(keep, cld_top_hgt[i,] > icing_alt[i], False)
-
- # Test3
- # keep = np.where(keep, np.invert((cld_opd[i,] >= opd_thick_threshold) & (cld_phase[i,] == 4) & (cld_top_hgt[i,] > icing_alt[i])), False)
+ keep = np.where(keep, np.logical_and(bt_11um[i,] > 228.0, bt_11um[i,] < 271.5), False)
mask.append(keep)
idxs.append(i)
@@ -2404,11 +2521,13 @@ def analyze(preds_file, labels, prob_avg, test_file):
print('report altitude (m): ', np.histogram(nda, bins=12))
+ flt_alt = nda.copy()
+
iint = np.where(iint == -1, 0, iint)
iint = np.where(iint != 0, 1, iint)
- nda[np.logical_and(nda >= 0, nda < 3400)] = 0
- nda[np.logical_and(nda >= 3400, nda < 15000)] = 1
+ nda[np.logical_and(nda >= 0, nda < 4000)] = 0
+ nda[np.logical_and(nda >= 4000, nda < 15000)] = 1
# nda[np.logical_and(nda >= 8000, nda < 15000)] = 2
#print(np.sum(nda == 0), np.sum(nda == 1), np.sum(nda == 2))
@@ -2437,6 +2556,7 @@ def analyze(preds_file, labels, prob_avg, test_file):
print('----------------------------------------------------------')
print('----------------------------------------------------------')
+ return flt_alt[iint == 0], flt_alt[iint == 1]
if prob_avg is None:
return
@@ -2495,6 +2615,8 @@ def analyze(preds_file, labels, prob_avg, test_file):
# print(np.nanmean(cld_dz[(nda == 2) & true_no_ice]), np.nanmean(cld_hgt[(nda == 2) & true_no_ice]), np.nanmean(cld_tmp[(nda == 2) & true_no_ice]))
# print('------------')
+ return flt_alt[iint == 0], flt_alt[iint == 1]
+
def get_training_parameters(day_night='DAY', l1b_andor_l2='both'):
if day_night == 'DAY':