diff --git a/modules/icing/pirep_goes.py b/modules/icing/pirep_goes.py
index 123f8a88ce6234554ee4958d93feedd14f5fc6bd..4ac49859b09e1ab09243200bb5989c0d0f7a3714 100644
--- a/modules/icing/pirep_goes.py
+++ b/modules/icing/pirep_goes.py
@@ -199,7 +199,7 @@ def get_grid_values(h5f, grid_name, j_c, i_c, half_width, num_j=None, num_i=None
return grd_vals
-def create_file(filename, data_dct, ds_list, ds_types, lon_c, lat_c, time_s, fl_alt_s, icing_intensity, unq_ids):
+def create_file(filename, data_dct, ds_list, ds_types, lon_c, lat_c, time_s, fl_alt_s, icing_intensity, unq_ids, mask=None):
h5f_expl = h5py.File(a_clvr_file, 'r')
h5f = h5py.File(filename, 'w')
@@ -234,6 +234,14 @@ def create_file(filename, data_dct, ds_list, ds_types, lon_c, lat_c, time_s, fl_
unq_ids_ds = h5f.create_dataset('unique_id', data=unq_ids, dtype='i4')
unq_ids_ds.attrs.create('long_name', data='ID mapping to PIREP icing dictionary: see pireps.py')
+ if mask is not None:
+ mask = mask.astype(np.byte)
+ mask_ds = h5f.create_dataset('FOV_mask', data=mask, dtype=np.byte)
+ mask_ds.attrs.create('long_name', data='The FOVs which pass the cloudy icing report test')
+ mask_ds.dims[0].label = 'time'
+ mask_ds.dims[1].label = 'y'
+ mask_ds.dims[2].label = 'x'
+
# copy relevant attributes
for ds_name in ds_list:
h5f_ds = h5f[ds_name]
@@ -935,8 +943,10 @@ def analyze2(filename, filename_l1b):
# --------------------------------------------
-x_a = 10
-x_b = 30
+# x_a = 10
+# x_b = 30
+x_a = 12
+x_b = 28
y_a = x_a
y_b = x_b
nx = ny = (x_b - x_a)
@@ -1025,24 +1035,34 @@ def run_qc(filename, filename_l1b, day_night='ANY', pass_thresh_frac=0.20, cloud
mask, idxs, num_tested = apply_qc_icing_pireps_v3(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 = []
+ mask_keep_idxs = []
for i in range(len(mask)):
- 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])
+ # keep_idxs.append(idxs[i])
+ frac = np.sum(mask[i]) / num_tested[i]
+ if icing:
+ if frac > pass_thresh_frac:
+ keep_idxs.append(idxs[i])
+ mask_keep_idxs.append(i)
+ elif frac > pass_thresh_frac:
+ keep_idxs.append(idxs[i])
+ mask_keep_idxs.append(i)
+
+ if len(keep_idxs) == 0:
+ f.close()
+ f_l1b.close()
+ return len(icing_alt), len(mask), len(keep_idxs)
print('day_night, icing, all, valid, pass: ', day_night, icing, len(icing_alt), len(mask), len(keep_idxs))
print('-----------------------')
keep_idxs = np.array(keep_idxs)
+ mask_keep_idxs = np.array(mask_keep_idxs)
+ mask = np.concatenate(mask).reshape((-1, ny, nx))
+ mask = mask[mask_keep_idxs,]
data_dct = {}
for didx, ds_name in enumerate(ds_list):
- data_dct[ds_name] = f[ds_name][keep_idxs,]
+ data_dct[ds_name] = f[ds_name][keep_idxs, y_a:y_b, x_a:x_b]
lon_c = f['longitude'][keep_idxs]
lat_c = f['latitude'][keep_idxs]
@@ -1055,17 +1075,17 @@ def run_qc(filename, filename_l1b, day_night='ANY', pass_thresh_frac=0.20, cloud
fbase, fext = os.path.splitext(fname)
outfile = path + '/' + fbase + '_' + 'QC' + '_' + day_night + fext
- create_file(outfile, data_dct, ds_list, ds_types, lon_c, lat_c, time_s, fl_alt_s, ice_int_s, unq_ids)
+ create_file(outfile, data_dct, ds_list, ds_types, lon_c, lat_c, time_s, fl_alt_s, ice_int_s, unq_ids, mask)
data_dct = {}
for didx, ds_name in enumerate(l1b_ds_list):
- data_dct[ds_name] = f_l1b[ds_name][keep_idxs]
+ data_dct[ds_name] = f_l1b[ds_name][keep_idxs, y_a:y_b, x_a:x_b]
path, fname = os.path.split(filename_l1b)
fbase, fext = os.path.splitext(fname)
outfile_l1b = path + '/' + fbase + '_' + 'QC' + '_' + day_night + fext
- create_file(outfile_l1b, data_dct, l1b_ds_list, l1b_ds_types, lon_c, lat_c, time_s, fl_alt_s, ice_int_s, unq_ids)
+ create_file(outfile_l1b, data_dct, l1b_ds_list, l1b_ds_types, lon_c, lat_c, time_s, fl_alt_s, ice_int_s, unq_ids, mask)
f.close()
f_l1b.close()
@@ -1396,8 +1416,8 @@ def apply_qc_icing_pireps_exp2(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cl
def apply_qc_icing_pireps_v3(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 = 1500.0
+ closeness_top = 50.0 # meters
+ max_cld_depth = 800.0
max_altitude = 4500.0
max_cld_altitude = 4000.0
@@ -1440,12 +1460,12 @@ def apply_qc_icing_pireps_v3(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_
# if icing_alt[i,] > max_altitude:
# continue
- cth_mean = np.nanmean(cld_top_hgt[i,])
- if not (closeness_top+cth_mean > icing_alt[i,] > cth_mean - max_cld_depth):
- continue
+ # cth_mean = np.nanmean(cld_top_hgt[i,])
+ # if not (closeness_top+cth_mean > icing_alt[i,] > cth_mean - max_cld_depth):
+ # 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_hgt[i,] + closeness_top) > icing_alt[i,],
+ (cld_top_hgt[i,] - max_cld_depth) < 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)
@@ -1783,8 +1803,10 @@ def fov_extract(icing_files, no_icing_files, trnfile='/home/rink/fovs_l1b_train.
def tile_extract(icing_files, no_icing_files, trnfile='/home/rink/tiles_train.h5', vldfile='/home/rink/tiles_valid.h5', tstfile='/home/rink/tiles_test.h5', L1B_or_L2='L1B',
cld_mask_name='cloud_mask', augment=False, do_split=True, has_test=True):
# 16x16
- n_a, n_b = 12, 28
- m_a, m_b = 12, 28
+ # n_a, n_b = 12, 28
+ # m_a, m_b = 12, 28
+ n_a, n_b = 0, 16
+ m_a, m_b = 0, 16
# 10x10
# n_a, n_b = 15, 25
# m_a, m_b = 15, 25