diff --git a/modules/icing/pirep_goes.py b/modules/icing/pirep_goes.py index 19ead2ffac068b8bf055be55b6f9d8e13b804091..1c2b716fe3b24492682926275c0acebe7c560785 100644 --- a/modules/icing/pirep_goes.py +++ b/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)