snapshot...

modules/icing/pirep_goes.py

@@ -2,7 +2,7 @@ from icing.pireps import pirep_icing
 import numpy as np
 import pickle
 import os
-from util.util import get_time_tuple_utc, GenericException, add_time_range_to_filename
+from util.util import get_time_tuple_utc, GenericException, add_time_range_to_filename, is_night
 from aeolus.datasource import CLAVRx, GOESL1B
 from util.geos_nav import GEOSNavigation
 import h5py
@@ -35,7 +35,8 @@ ds_list = ['cld_height_acha', 'cld_geo_thick', 'cld_press_acha', 'sensor_zenith_
 ds_types = ['f4' for i in range(21)] + ['i4' for i in range(3)]
 
 
-a_clvr_file = '/home/rink/data/clavrx/clavrx_OR_ABI-L1b-RadC-M3C01_G16_s20190020002186.level2.nc'
+#a_clvr_file = '/home/rink/data/clavrx/clavrx_OR_ABI-L1b-RadC-M3C01_G16_s20190020002186.level2.nc'
+a_clvr_file = '/Users/tomrink/data/clavrx/clavrx_OR_ABI-L1b-RadC-M3C01_G16_s20190020002186.level2.nc'
 
 
 def setup():
@@ -593,7 +594,7 @@ def process_1(ice_dct, no_ice_dct, neg_ice_dct):
     return new_ice_dct, new_no_ice_dct, new_neg_ice_dct
 
 
-def run_qc(filename, filename_l1b, outfile, outfile_l1b):
+def run_qc(filename, filename_l1b, outfile, outfile_l1b, day_night='ANY'):
     f = h5py.File(filename, 'r')
     icing_alt = f['icing_altitude'][:]
     cld_top_hgt = f['cld_height_acha'][:, 10:30, 10:30]
@@ -601,33 +602,40 @@ def run_qc(filename, filename_l1b, outfile, outfile_l1b):
     cld_opd = f['cld_opd_acha'][:, 10:30, 10:30]
     cld_opd_dc = f['cld_opd_dcomp'][:, 10:30, 10:30]
     cld_mask = f['cloud_mask'][:, 10:30, 10:30]
+    sol_zen = 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]
 
     print('num pireps: ', len(icing_alt))
 
-    mask = apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_phase, cld_opd, cld_mask, bt_11um)
+    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))
 
-    bts = []
-    phs = []
-    opd = []
+    bts_pass = []
+    phs_pass = []
+    opd_pass = []
     opd_dc = []
-    keep_idxs =[]
+    keep_idxs = []
+    bts_reject = []
+    phs_reject = []
+    opd_reject = []
+
     for i in range(len(mask)):
         if (np.sum(mask[i]) / 400) > 0.20:
-            bts.append((bt_11um[i,].flatten())[mask[i]])
-            phs.append((cld_phase[i,].flatten())[mask[i]])
-            opd.append((cld_opd[i,].flatten())[mask[i]])
+            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]])
             #opd_dc.append(cld_opd_dc[i,].flatten())[mask[i]]
             keep_idxs.append(i)
-        #else:
-        #    bts.append((bt_11um[i,].flatten())[mask[i]])
-    print('num valid pireps: ', len(bts))
+        else:
+            bts_reject.append((bt_11um[i,].flatten())[mask[i]])
+
+    print('num valid pireps: ', len(bts_pass))
 
-    bts = np.concatenate(bts)
-    phs = np.concatenate(phs)
-    opd = np.concatenate(opd)
+    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))
@@ -659,10 +667,8 @@ def run_qc(filename, filename_l1b, outfile, outfile_l1b):
     f.close()
     f_l1b.close()
 
-    return mask
-
 
-def apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_phase, cld_opd, cld_mask, bt_11um):
+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
     closeness = 100.0  # meters
     num_obs = len(icing_alt)
@@ -672,17 +678,18 @@ def apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_phase, cld_opd, cld_mask,
     cld_opd = cld_opd.reshape((num_obs, -1))
     bt_11um = bt_11um.reshape((num_obs, -1))
 
-    skip = False
     mask = []
     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
+
         keep_0 = np.logical_or(cld_mask[i,] == 2, cld_mask[i,] == 3)  # cloudy
         keep_1 = np.invert(np.isnan(cld_top_hgt[i,]))
         keep_2 = np.invert(np.isnan(bt_11um[i,]))
         keep_3 = np.invert(np.isnan(cld_opd[i,]))
         keep = keep_0 & keep_1 & keep_2 & keep_3
-        if skip:
-            mask.append(keep)
-            continue
 
         keep = np.where(keep, cld_top_hgt[i,] > icing_alt[i], False)