diff --git a/modules/icing/pirep_goes.py b/modules/icing/pirep_goes.py index 7b3c9b5547dac929b6e496a269cf765891acbe04..b8852348a05443af939d52272fb909f1650f0288 100644 --- a/modules/icing/pirep_goes.py +++ b/modules/icing/pirep_goes.py @@ -2323,7 +2323,7 @@ def run_icing_predict(clvrx_dir='/Users/tomrink/data/clavrx/RadC/', output_dir=h h5f.close() -def run_icing_predict_new(clvrx_dir='/Users/tomrink/data/clavrx/RadC/', output_dir=homedir, +def run_icing_predict_fcn(clvrx_dir='/Users/tomrink/data/clavrx/RadC/', output_dir=homedir, day_model_path=model_path_day, night_model_path=model_path_night, prob_thresh=0.5, satellite='GOES16', domain='CONUS', day_night='AUTO', l1b_andor_l2='both', use_flight_altitude=False, res_fac=1, use_nan=False): @@ -2586,3 +2586,137 @@ def run_icing_predict_image(clvrx_dir='/Users/tomrink/data/clavrx/RadC/', output print('Done: ', clvrx_str_time) h5f.close() + + +def run_icing_predict_image_fcn(clvrx_dir='/Users/tomrink/data/clavrx/RadC/', output_dir=homedir, + day_model_path=model_path_day, night_model_path=model_path_night, + prob_thresh=0.5, satellite='GOES16', domain='CONUS', day_night='AUTO', + l1b_andor_l2='BOTH', use_flight_altitude=True, res_fac=1, + extent=[-105, -70, 15, 50], + pirep_file='/Users/tomrink/data/pirep/pireps_202109200000_202109232359.csv', + obs_lons=None, obs_lats=None, obs_times=None, obs_alt=None, flight_level=None): + + if use_flight_altitude is True: + flight_levels = [0, 1, 2, 3, 4] + else: + flight_levels = [0] + + if pirep_file is not None: + ice_dict, no_ice_dict, neg_ice_dict = setup(pirep_file) + + alt_lo, alt_hi = 0.0, 15000.0 + if flight_level is not None: + alt_lo, alt_hi = flt_level_ranges[flight_level] + + day_train_params = get_training_parameters(day_night='DAY', l1b_andor_l2=l1b_andor_l2) + nght_train_params = get_training_parameters(day_night='NIGHT', l1b_andor_l2=l1b_andor_l2) + + if day_night == 'AUTO': + train_params = list(set(day_train_params + nght_train_params)) + elif day_night == 'DAY': + train_params = day_train_params + elif day_night == 'NIGHT': + train_params = nght_train_params + + if satellite == 'H08': + clvrx_ds = CLAVRx_H08(clvrx_dir) + else: + clvrx_ds = CLAVRx(clvrx_dir) + clvrx_files = clvrx_ds.flist + + for fidx, fname in enumerate(clvrx_files): + h5f = h5py.File(fname, 'r') + dto = clvrx_ds.get_datetime(fname) + ts = dto.timestamp() + clvrx_str_time = dto.strftime('%Y-%m-%d_%H:%M') + + dto, _ = get_time_tuple_utc(ts) + dto_0 = dto - datetime.timedelta(minutes=30) + dto_1 = dto + datetime.timedelta(minutes=30) + ts_0 = dto_0.timestamp() + ts_1 = dto_1.timestamp() + + if pirep_file is not None: + _, keep_lons, keep_lats, _ = time_filter_3(ice_dict, ts_0, ts_1, alt_lo, alt_hi) + elif obs_times is not None: + keep = np.logical_and(obs_times >= ts_0, obs_times < ts_1) + keep = np.where(keep, np.logical_and(obs_alt >= alt_lo, obs_alt < alt_hi), False) + keep_lons = obs_lons[keep] + keep_lats = obs_lats[keep] + else: + keep_lons = None + keep_lats = None + + data_dct, solzen, satzen, ll, cc = prepare_evaluate(h5f, name_list=train_params, satellite=satellite, domain=domain, offset=8) + num_elems = len(cc) + num_lines = len(ll) + + if fidx == 0: + nav = get_navigation(satellite, domain) + lons_2d, lats_2d, x_rad, y_rad = get_lon_lat_2d_mesh(nav, ll, cc) + + day_idxs = solzen < 80.0 + num_day_tiles = np.sum(day_idxs) + + nght_idxs = solzen > 100.0 + num_nght_tiles = np.sum(nght_idxs) + + # initialize output arrays + probs_2d_dct = {flvl: None for flvl in flight_levels} + preds_2d_dct = {flvl: None for flvl in flight_levels} + for flvl in flight_levels: + fd_preds = np.zeros(num_lines * num_elems, dtype=np.int8) + fd_preds[:] = -1 + fd_probs = np.zeros(num_lines * num_elems, dtype=np.float32) + fd_probs[:] = -1.0 + preds_2d_dct[flvl] = fd_preds + probs_2d_dct[flvl] = fd_probs + + if (day_night == 'AUTO' or day_night == 'DAY') and num_day_tiles > 0: + + preds_day_dct, probs_day_dct = run_evaluate_static_fcn(data_dct, day_model_path, + day_night='DAY', l1b_or_l2=l1b_andor_l2, + prob_thresh=prob_thresh, + use_flight_altitude=use_flight_altitude, + flight_levels=flight_levels) + for flvl in flight_levels: + preds = preds_day_dct[flvl] + probs = probs_day_dct[flvl] + fd_preds = preds_2d_dct[flvl] + fd_probs = probs_2d_dct[flvl] + fd_preds[day_idxs] = preds[day_idxs] + fd_probs[day_idxs] = probs[day_idxs] + + if (day_night == 'AUTO' or day_night == 'NIGHT') and num_nght_tiles > 0: + preds_nght_dct, probs_nght_dct = run_evaluate_static_fcn(data_dct, night_model_path, + day_night='NIGHT', l1b_or_l2=l1b_andor_l2, + prob_thresh=prob_thresh, + use_flight_altitude=use_flight_altitude, + flight_levels=flight_levels) + for flvl in flight_levels: + preds = preds_nght_dct[flvl] + probs = probs_nght_dct[flvl] + fd_preds = preds_2d_dct[flvl] + fd_probs = probs_2d_dct[flvl] + fd_preds[nght_idxs] = preds[nght_idxs] + fd_probs[nght_idxs] = probs[nght_idxs] + + for flvl in flight_levels: + fd_preds = preds_2d_dct[flvl] + fd_probs = probs_2d_dct[flvl] + preds_2d_dct[flvl] = fd_preds.reshape((num_lines, num_elems)) + probs_2d_dct[flvl] = fd_probs.reshape((num_lines, num_elems)) + + prob_s = [] + for flvl in flight_levels: + probs = probs_2d_dct[flvl] + prob_s.append(probs) + prob_s = np.stack(prob_s, axis=-1) + max_prob = np.max(prob_s, axis=2) + max_prob = np.where(max_prob < 0.5, np.nan, max_prob) + + make_icing_image(h5f, max_prob, None, None, clvrx_str_time, satellite, domain, + ice_lons_vld=keep_lons, ice_lats_vld=keep_lats, extent=extent) + + print('Done: ', clvrx_str_time) + h5f.close() \ No newline at end of file