snapshot...

modules/icing/pirep_goes.py

@@ -607,10 +607,9 @@ def run_qc(filename, filename_l1b, outfile, outfile_l1b, day_night='ANY'):
     f_l1b = h5py.File(filename_l1b, 'r')
     bt_11um = f_l1b['temp_11_0um_nom'][:, 10:30, 10:30]
 
-    print('num pireps: ', len(icing_alt))
-
-    mask = apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_phase, cld_opd, cld_mask, bt_11um, sol_zen, day_night=day_night)
-    print('num pireps: ', len(mask))
+    print('num pireps all: ', len(icing_alt))
+    mask, idxs = apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_phase, cld_opd, cld_mask, bt_11um, sol_zen, day_night=day_night)
+    print('num pireps day_night: ', len(mask), day_night)
 
     bts_pass = []
     phs_pass = []
@@ -623,26 +622,15 @@ def run_qc(filename, filename_l1b, outfile, outfile_l1b, day_night='ANY'):
 
     for i in range(len(mask)):
         if (np.sum(mask[i]) / 400) > 0.20:
-            bts_pass.append((bt_11um[i,].flatten())[mask[i]])
-            phs_pass.append((cld_phase[i,].flatten())[mask[i]])
-            opd_pass.append((cld_opd[i,].flatten())[mask[i]])
+            bts_pass.append((bt_11um[idxs[i],].flatten())[mask[i]])
+            phs_pass.append((cld_phase[idxs[i],].flatten())[mask[i]])
+            opd_pass.append((cld_opd[idxs[i],].flatten())[mask[i]])
             #opd_dc.append(cld_opd_dc[i,].flatten())[mask[i]]
-            keep_idxs.append(i)
+            keep_idxs.append(idxs[i])
         else:
-            bts_reject.append((bt_11um[i,].flatten())[mask[i]])
-
-    print('num valid pireps: ', len(bts_pass))
-
-    bts = np.concatenate(bts_pass)
-    phs = np.concatenate(phs_pass)
-    opd = np.concatenate(opd_pass)
-    #opd_dc = np.concatenate(opd_dc)
-
-    print(np.histogram(bts, bins=10))
-    print(np.histogram(opd, bins=10))
-    #print(np.histogram(opd_dc, bins=10))
-    print(np.histogram(phs, bins=6))
+            bts_reject.append((bt_11um[idxs[i],].flatten())[mask[i]])
 
+    print('num valid pireps: ', len(keep_idxs))
     keep_idxs = np.array(keep_idxs)
 
     data_dct = {}
@@ -656,17 +644,35 @@ def run_qc(filename, filename_l1b, outfile, outfile_l1b, day_night='ANY'):
     ice_int_s = f['icing_intensity'][keep_idxs]
     unq_ids = f['unique_id'][keep_idxs]
 
+    path, fname = os.path.split(filename)
+    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)
 
     data_dct = {}
     for didx, ds_name in enumerate(l1b_ds_list):
         data_dct[ds_name] = f_l1b[ds_name][keep_idxs]
 
+    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)
 
     f.close()
     f_l1b.close()
 
+    bts = np.concatenate(bts_pass)
+    phs = np.concatenate(phs_pass)
+    opd = np.concatenate(opd_pass)
+    #opd_dc = np.concatenate(opd_dc)
+
+    print(np.histogram(bts, bins=10))
+    print(np.histogram(opd, bins=10))
+    #print(np.histogram(opd_dc, bins=10))
+    print(np.histogram(phs, bins=6))
+
 
 def apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_phase, cld_opd, cld_mask, bt_11um, solzen, day_night='ANY'):
     opd_threshold = 2
@@ -679,11 +685,14 @@ def apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_phase, cld_opd, cld_mask,
     bt_11um = bt_11um.reshape((num_obs, -1))
 
     mask = []
+    idxs = []
     for i in range(num_obs):
-        if day_night == 'NIGHT' and not is_night(solzen[i,]):
-            continue
-        elif day_night == 'DAY' and is_night(solzen[i,]):
-            continue
+        if day_night == 'NIGHT':
+            if not is_night(solzen[i,]):
+                continue
+        elif day_night == 'DAY':
+            if 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,]))
@@ -707,5 +716,6 @@ def apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_phase, cld_opd, cld_mask,
         keep = np.where(keep, np.invert(bt_11um[i,] < 228.0), False)
 
         mask.append(keep)
+        idxs.append(i)
 
-    return mask
+    return mask, idxs