snapshot...

modules/icing/pirep_goes.py

@@ -302,15 +302,52 @@ def analyze(ice_dct, no_ice_dct):
         print(dt_str[2:])
 
 
-def apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_phase, cld_opd, bt_11um, cld_mask):
+def run_qc(filename, filename_l1b):
+    f = h5py.File(filename, 'r')
+    icing_alt = f['icing_altitude'][:]
+    cld_top_hgt = f['cld_height_acha'][:, 10:30, 10:30]
+    cld_phase = f['cloud_phase'][:, 10:30, 10:30]
+    cld_opd = f['cld_opd_acha'][:, 10:30, 10:30]
+    cld_mask = f['cloud_mask'][:, 10:30, 10:30]
+
+    f_l1b = h5py.File(filename_l1b, 'r')
+    bt_11um = f_l1b['temp_11_0um_nom'][:, 10:30, 10:30]
+
+    mask = apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_phase, cld_opd, cld_mask, bt_11um)
+
+    f.close()
+    f_l1b.close()
+
+    cnt = 0
+    bts = []
+    phs = []
+    for i in range(len(mask)):
+        if (np.sum(mask[i]) / 400) > 0.15:
+            cnt += 1
+            bts.append((bt_11um[i,].flatten())[mask[i]])
+            phs.append((cld_phase[i,].flatten())[mask[i]])
+        #else:
+        #    bts.append((bt_11um[i,].flatten())[mask[i]])
+    print(cnt)
+
+    bts = np.concatenate(bts)
+    phs = np.concatenate(phs)
+    print(bts.shape)
+    print(np.histogram(bts, bins=20))
+    print(np.histogram(phs, bins=6))
+
+    return mask
+
+
+def apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_phase, cld_opd, cld_mask, bt_11um):
     opd_threshold = 2
     closeness = 100.0  # meters
     num_obs = len(icing_alt)
-    cld_mask = cld_mask[:, :, :].reshape((num_obs, -1))
-    cld_top_hgt = cld_top_hgt[:, :, :].reshape((num_obs, -1))
-    cld_phase = cld_phase[:, :, :].reshape((num_obs, -1))
-    cld_opd = cld_opd[:, :, :].reshape((num_obs, -1))
-    bt_11um = bt_11um[:, :, :].reshape((num_obs, -1))
+    cld_mask = cld_mask.reshape((num_obs, -1))
+    cld_top_hgt = cld_top_hgt.reshape((num_obs, -1))
+    cld_phase = cld_phase.reshape((num_obs, -1))
+    cld_opd = cld_opd.reshape((num_obs, -1))
+    bt_11um = bt_11um.reshape((num_obs, -1))
 
     mask = []
     for i in range(num_obs):
@@ -323,13 +360,13 @@ def apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_phase, cld_opd, bt_11um, c
         keep = np.where(keep, cld_top_hgt[i,] > icing_alt[i], False)
 
         keep = np.where(keep,
-            np.invert((cld_phase[i,] == 3) &
+            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)
 
         keep = np.where(keep, (cld_opd[i,] >= opd_threshold) & (cld_phase[i,] == 4) & (cld_top_hgt[i,] > icing_alt[i]), False)
 
-        keep = np.where(keep, np.invert((cld_phase[i,] == 3) & (cld_opd[i,] < 0.1) & (cld_top_hgt[i,] > icing_alt[i])), False)
+        keep = np.where(keep, np.invert((cld_phase[i,] == 4) & (cld_opd[i,] < 0.1) & (cld_top_hgt[i,] > icing_alt[i])), False)
 
         keep = np.where(keep, np.invert(bt_11um[i,] > 270.0), False)