diff --git a/modules/icing/pirep_goes.py b/modules/icing/pirep_goes.py
index a642cf0081ac030c50ac5a48dd8ea2a47874b283..49e908e41c8d00f41cbdeac8a753f6515460955b 100644
--- a/modules/icing/pirep_goes.py
+++ b/modules/icing/pirep_goes.py
@@ -1024,11 +1024,11 @@ def run_qc(filename, filename_l1b, day_night='ANY', pass_thresh_frac=0.20, cloud
bt_11um = f_l1b['temp_11_0um_nom'][:, y_a:y_b, x_a:x_b]
if icing:
- # mask, idxs, num_tested = apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_opd, cld_mask, bt_11um, sol_zen, sat_zen, cld_temp, day_night=day_night, cloud_frac=cloudy_frac)
- mask, idxs, num_tested = apply_qc_icing_pireps_exp1(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_opd, cld_mask, bt_11um, sol_zen, sat_zen, cld_temp, day_night=day_night, cloud_frac=cloudy_frac)
+ # mask, idxs, num_tested = apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_opd, cld_mask, bt_11um, sol_zen, sat_zen, cld_temp, day_night=day_night, cloudy_frac=cloudy_frac)
+ mask, idxs, num_tested = apply_qc_icing_pireps_exp1(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_opd, cld_mask, bt_11um, sol_zen, sat_zen, cld_temp, day_night=day_night, cloudy_frac=cloudy_frac)
else:
# mask, idxs, num_tested = apply_qc_no_icing_pireps(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_opd, cld_mask, bt_11um, sol_zen, sat_zen, cld_temp, day_night=day_night, cloudy_frac=cloudy_frac)
- mask, idxs, num_tested = apply_qc_icing_pireps_exp1(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_opd, cld_mask, bt_11um, sol_zen, sat_zen, cld_temp, day_night=day_night, cloud_frac=cloudy_frac)
+ mask, idxs, num_tested = apply_qc_icing_pireps_exp1(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_opd, cld_mask, bt_11um, sol_zen, sat_zen, cld_temp, day_night=day_night, cloudy_frac=cloudy_frac)
keep_idxs = []
@@ -1213,7 +1213,7 @@ def run_qc(filename, filename_l1b, day_night='ANY', pass_thresh_frac=0.20, cloud
# 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, cld_top_temp, day_night='ANY', cloud_frac=0.5):
+def apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_opd, cld_mask, bt_11um, solzen, satzen, cld_top_temp, day_night='ANY', cloudy_frac=0.5):
if day_night == 'DAY':
opd_thick_threshold = 20
@@ -1292,7 +1292,7 @@ def apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_opd
return mask, idxs, num_tested
-def apply_qc_icing_pireps_exp1(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_opd, cld_mask, bt_11um, solzen, satzen, cld_top_temp, day_night='ANY', cloud_frac=0.5):
+def apply_qc_icing_pireps_exp1(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_opd, cld_mask, bt_11um, solzen, satzen, cld_top_temp, day_night='ANY', cloudy_frac=0.5):
if day_night == 'DAY':
opd_thick_threshold = 20
@@ -1349,14 +1349,7 @@ def apply_qc_icing_pireps_exp1(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cl
return mask, idxs, num_tested
-def apply_qc_icing_pireps_exp2(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_opd, cld_mask, bt_11um, solzen, satzen, cld_top_temp, day_night='ANY', cloud_frac=0.5):
-
- 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
+def apply_qc_icing_pireps_exp2(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_opd, cld_mask, bt_11um, solzen, satzen, cld_top_temp, day_night='ANY', cloudy_frac=0.5):
closeness_top = 200.0 # meters
max_cld_depth = 1000.0
@@ -1485,6 +1478,56 @@ def apply_qc_no_icing_pireps(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_
return mask, idxs, num_tested
+def apply_qc_no_icing_pireps_exp2(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_opd, cld_mask, bt_11um, solzen, satzen, cld_top_temp, day_night='ANY', cloudy_frac=0.5):
+
+ closeness_top = 200.0 # meters
+ max_cld_depth = 1000.0
+ max_altitude = 4000.0
+ max_cld_altitude = 4000.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))
+ cld_top_temp = cld_top_temp.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
+
+ # keep_0 = np.logical_or(cld_mask[i,] == 2, cld_mask[i,] == 3) # cloudy
+ keep_0 = cld_mask[i, ] == 3 # Confident 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
+ if num_keep / nx_x_ny < cloud_frac: # At least this fraction cloudy
+ continue
+
+ keep = np.where(keep, np.logical_and(cld_top_hgt[i,] > icing_alt[i,],
+ (cld_top_hgt[i,] - 3000.0) < icing_alt[i,]), False)
+ # keep = np.where(keep, np.logical_and(cld_top_temp[i,] > 228.0, cld_top_temp[i,] < 274.0), False)
+
+ mask.append(keep)
+ idxs.append(i)
+ num_tested.append(num_keep)
+
+ return mask, idxs, num_tested
+
+
def fov_extract(icing_files, no_icing_files, trnfile='/home/rink/fovs_l1b_train.h5', tstfile='/home/rink/fovs_l1b_test.h5', L1B_or_L2='L1B', split=0.2):
ice_times = []
icing_int_s = []