snapshot...

modules/icing/pirep_goes.py

 from icing.pireps import pirep_icing
 import numpy as np
 import pickle
+import matplotlib.pyplot as plt
 import os
 from util.util import get_time_tuple_utc, GenericException, add_time_range_to_filename, is_night
 from aeolus.datasource import CLAVRx, GOESL1B
@@ -594,6 +595,35 @@ def process_1(ice_dct, no_ice_dct, neg_ice_dct):
     return new_ice_dct, new_no_ice_dct, new_neg_ice_dct
 
 
+def analyze2(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_dc = f['cld_opd_dcomp'][:, 10:30, 10:30]
+    cld_opd = f['cld_opd_acha'][:, 10:30, 10:30]
+    solzen = f['solar_zenith_angle'][:, 10:30, 10:30]
+
+    f_l1b = h5py.File(filename_l1b, 'r')
+    bt_11um = f_l1b['temp_11_0um_nom'][:, 10:30, 10:30]
+
+    cld_opd = cld_opd.flatten()
+    cld_opd_dc = cld_opd_dc.flatten()
+    solzen = solzen.flatten()
+
+    keep1 = np.invert(np.isnan(cld_opd))
+    keep2 = np.invert(np.isnan(solzen))
+    keep = keep1 & keep2
+    cld_opd = cld_opd[np.invert(np.isnan(cld_opd))]
+    cld_opd_dc = cld_opd_dc[keep]
+    solzen = solzen[keep]
+
+    plt.hist(cld_opd, bins=20)
+    plt.show()
+    plt.hist(cld_opd_dc, bins=20)
+    plt.show()
+
+
 def run_qc(filename, filename_l1b, day_night='ANY'):
     f = h5py.File(filename, 'r')
     icing_alt = f['icing_altitude'][:]
@@ -610,7 +640,9 @@ def run_qc(filename, filename_l1b, day_night='ANY'):
     bt_11um = f_l1b['temp_11_0um_nom'][:, 10:30, 10:30]
 
     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 = []
@@ -687,7 +719,7 @@ def run_qc(filename, filename_l1b, day_night='ANY'):
 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
     if day_night == 'DAY':
-        opd_threshold = 5
+        opd_threshold = 20
     closeness = 100.0  # meters
     num_obs = len(icing_alt)
     cld_mask = cld_mask.reshape((num_obs, -1))