diff --git a/modules/icing/pirep_goes.py b/modules/icing/pirep_goes.py
index 31ec0e11f559c8a730c1d749d3a7e63aaf4ec4e5..37c78d7ca99b2162a8d0ababa1e6b73797c6115c 100644
--- a/modules/icing/pirep_goes.py
+++ b/modules/icing/pirep_goes.py
@@ -1019,6 +1019,7 @@ def run_qc(filename, filename_l1b, day_night='ANY', pass_thresh_frac=0.20, icing
cld_top_hgt = f['cld_height_acha'][:, y_a:y_b, x_a:x_b]
cld_geo_dz = f['cld_geo_thick'][:, y_a:y_b, x_a:x_b]
cld_phase = f['cloud_phase'][:, y_a:y_b, x_a:x_b]
+ cld_temp = f['cld_temp_acha'][:, y_a:y_b, x_a:x_b]
if day_night == 'DAY':
cld_opd = f['cld_opd_dcomp'][:, y_a:y_b, x_a:x_b]
@@ -1033,20 +1034,21 @@ def run_qc(filename, filename_l1b, day_night='ANY', pass_thresh_frac=0.20, icing
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, day_night=day_night)
+ 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)
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, day_night=day_night)
+ 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)
keep_idxs = []
for i in range(len(mask)):
- # frac = np.sum(mask[i]) / nx_x_ny
- frac = np.sum(mask[i]) / num_tested[i]
- if icing:
- if frac > pass_thresh_frac:
- keep_idxs.append(idxs[i])
- elif frac > pass_thresh_frac:
- keep_idxs.append(idxs[i])
+ keep_idxs.append(idxs[i])
+ # # frac = np.sum(mask[i]) / nx_x_ny
+ # frac = np.sum(mask[i]) / num_tested[i]
+ # if icing:
+ # if frac > pass_thresh_frac:
+ # keep_idxs.append(idxs[i])
+ # elif frac > pass_thresh_frac:
+ # keep_idxs.append(idxs[i])
print('day_night, icing, all, valid, pass: ', day_night, icing, len(icing_alt), len(mask), len(keep_idxs))
print('-----------------------')
@@ -1220,7 +1222,7 @@ def run_qc(filename, filename_l1b, day_night='ANY', pass_thresh_frac=0.20, icing
# 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', cloud_fraction=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_temp, day_night='ANY', cloud_fraction=0.5):
if day_night == 'DAY':
opd_thick_threshold = 20
@@ -1230,9 +1232,9 @@ def apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_opd
opd_thin_threshold = 0.1
closeness_top = 100.0 # meters
- max_cld_depth = 3000.0
+ max_cld_depth = 1500.0
max_altitude = 5000.0
- max_cld_altitude = 6000.0
+ max_cld_altitude = 5000.0
num_obs = len(icing_alt)
cld_mask = cld_mask.reshape((num_obs, -1))
@@ -1255,8 +1257,8 @@ 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 (icing_alt[i] < max_altitude):
- # continue
+ if not (icing_alt[i] < max_altitude):
+ 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,]))
@@ -1269,13 +1271,23 @@ def apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_opd
if num_keep / nx_x_ny < cloud_fraction: # At least this fraction cloudy
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,] + closeness_top) > icing_alt[i], False)
# keep = np.where(keep, cld_top_hgt[i,] < max_cld_altitude, False)
# keep = np.where(keep, (cld_top_hgt[i,] - max_cld_depth) < icing_alt[i], False)
# keep = np.where(keep, np.logical_and(bt_11um[i,] > 228.0, bt_11um[i,] < 270.0), False)
- keep = np.where(keep, bt_11um[i,] < 275.0, False)
+ # keep = np.where(keep, bt_11um[i,] < 275.0, False)
+
+ med_cld_hgt = np.median(cld_top_hgt[i, keep])
+ if icing_alt[i] > med_cld_hgt:
+ continue
+
+ med_cld_tmp = np.median(cld_temp[i, keep])
+ if med_cld_tmp > 274.5:
+ continue
+
+ if icing_alt < med_cld_hgt - max_cld_depth:
+ continue
mask.append(keep)
idxs.append(i)
@@ -1284,7 +1296,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_no_icing_pireps(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_opd, cld_mask, bt_11um, solzen, satzen, day_night='ANY', cloudy_frac=0.5):
+def apply_qc_no_icing_pireps(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_opd, cld_mask, bt_11um, solzen, satzen, cld_temp, day_night='ANY', cloudy_frac=0.5):
if day_night == 'DAY':
opd_thick_threshold = 20
@@ -1294,9 +1306,9 @@ def apply_qc_no_icing_pireps(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_
opd_thin_threshold = 0.1
closeness_top = 100.0 # meters
- max_cld_depth = 3000.0
+ max_cld_depth = 1500.0
max_altitude = 5000.0
- max_cld_altitude = 6000.0
+ max_cld_altitude = 5000.0
num_obs = len(icing_alt)
cld_mask = cld_mask.reshape((num_obs, -1))
@@ -1319,8 +1331,8 @@ 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 (icing_alt[i] < max_altitude):
- # continue
+ if not (icing_alt[i] < max_altitude):
+ 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,]))
@@ -1333,12 +1345,23 @@ def apply_qc_no_icing_pireps(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_
if num_keep / nx_x_ny < cloudy_frac: # At least this fraction cloudy
continue
- keep = np.where(keep, (cld_top_hgt[i,] + closeness_top) > icing_alt[i], False)
+ # keep = np.where(keep, (cld_top_hgt[i,] + closeness_top) > icing_alt[i], False)
# keep = np.where(keep, cld_top_hgt[i,] < max_cld_altitude, False)
# keep = np.where(keep, (cld_top_hgt[i,] - max_cld_depth) < icing_alt[i], False)
# keep = np.where(keep, np.logical_and(bt_11um[i,] > 228.0, bt_11um[i,] < 270.0), False)
- keep = np.where(keep, bt_11um[i,] < 275.0, False)
+ # keep = np.where(keep, bt_11um[i,] < 275.0, False)
+
+ med_cld_hgt = np.median(cld_top_hgt[i, keep])
+ if icing_alt[i] > med_cld_hgt:
+ continue
+
+ med_cld_tmp = np.median(cld_temp[i, keep])
+ if med_cld_tmp > 274.5:
+ continue
+
+ if icing_alt < med_cld_hgt - max_cld_depth:
+ continue
mask.append(keep)
idxs.append(i)