From c007a8b29f26650801f13175ea81df70568edf80 Mon Sep 17 00:00:00 2001
From: tomrink <rink@ssec.wisc.edu>
Date: Fri, 10 May 2024 10:39:00 -0500
Subject: [PATCH] snapshot...
---
modules/icing/util.py | 201 ++++++++++++++++++++++++++++++++++++++++++
1 file changed, 201 insertions(+)
diff --git a/modules/icing/util.py b/modules/icing/util.py
index 86b6feea..06fcc7cf 100644
--- a/modules/icing/util.py
+++ b/modules/icing/util.py
@@ -11,6 +11,7 @@ from aeolus.datasource import CLAVRx, CLAVRx_VIIRS, CLAVRx_H08, CLAVRx_H09
import h5py
import datetime
from netCDF4 import Dataset
+import joblib
import tensorflow as tf
import os
# from scipy.signal import medfilt2d
@@ -332,6 +333,59 @@ def prepare_evaluate(h5f, name_list, satellite='GOES16', domain='FD', res_fac=1,
return grd_dct_n, solzen, satzen, ll, cc
+def prepare_evaluate1x1(h5f, name_list, satellite='GOES16', domain='FD', res_fac=1, offset=0):
+ w_x = 1
+ w_y = 1
+ i_0 = 0
+ j_0 = 0
+ s_x = w_x // res_fac
+ s_y = w_y // res_fac
+
+ geos, xlen, xmin, xmax, ylen, ymin, ymax = get_cartopy_crs(satellite, domain)
+ if satellite == 'H08':
+ xlen = taiwan_lenx
+ ylen = taiwan_leny
+ i_0 = taiwan_i0
+ j_0 = taiwan_j0
+ elif satellite == 'H09':
+ xlen = taiwan_lenx
+ ylen = taiwan_leny
+ i_0 = taiwan_i0
+ j_0 = taiwan_j0
+
+ n_x = (xlen // s_x)
+ n_y = (ylen // s_y)
+
+ ll = [(offset+j_0) + j*s_y for j in range(n_y)]
+ cc = [(offset+i_0) + i*s_x for i in range(n_x)]
+
+ grd_dct_n = {name: [] for name in name_list}
+
+ cnt_a = 0
+ for ds_name in name_list:
+ fill_value, fill_value_name = get_fill_attrs(ds_name)
+ gvals = get_grid_values(h5f, ds_name, j_0, i_0, None, num_j=ylen, num_i=xlen, fill_value_name=fill_value_name, fill_value=fill_value)
+ print(gvals.shape)
+ gvals = np.expand_dims(gvals, axis=0)
+ print(gvals.shape)
+ print('--------------------')
+ if gvals is not None:
+ grd_dct_n[ds_name] = gvals
+ cnt_a += 1
+
+ if cnt_a > 0 and cnt_a != len(name_list):
+ raise GenericException('weirdness')
+
+ solzen = get_grid_values(h5f, 'solar_zenith_angle', j_0, i_0, None, num_j=ylen, num_i=xlen)
+ satzen = get_grid_values(h5f, 'sensor_zenith_angle', j_0, i_0, None, num_j=ylen, num_i=xlen)
+ cldmsk = get_grid_values(h5f, 'cloud_mask', j_0, i_0, None, num_j=ylen, num_i=xlen)
+ solzen = solzen[0:n_y*s_y:s_y, 0:n_x*s_x:s_x]
+ satzen = satzen[0:n_y*s_y:s_y, 0:n_x*s_x:s_x]
+ cldmsk = cldmsk[0:n_y*s_y:s_y, 0:n_x*s_x:s_x]
+
+ return grd_dct_n, solzen, satzen, cldmsk, ll, cc
+
+
flt_level_ranges_str = {k: None for k in range(6)}
flt_level_ranges_str[0] = '0_2000'
flt_level_ranges_str[1] = '2000_4000'
@@ -1422,3 +1476,150 @@ def run_icing_predict_image_fcn(clvrx_dir='/Users/tomrink/data/clavrx/RadC/', ou
print('Done: ', clvrx_str_time)
h5f.close()
+
+def run_icing_predict_image_1x1(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, use_max_cth_level=False,
+ 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, obs_intensity=None):
+
+ import deeplearning.icing_fcn as icing_fcn
+ # model_module = icing_fcn
+ #
+ # if day_model_path is not None:
+ # day_model = model_module.load_model(day_model_path, day_night='DAY', l1b_andor_l2=l1b_andor_l2,
+ # use_flight_altitude=use_flight_altitude)
+ # if night_model_path is not None:
+ # night_model = model_module.load_model(night_model_path, day_night='NIGHT', l1b_andor_l2=l1b_andor_l2,
+ # use_flight_altitude=use_flight_altitude)
+
+ # load parameter stats and model from disk
+ stdSclr = joblib.load('/home/rink/stdSclr_4.pkl')
+ model = joblib.load('/home/rink/icing_gbm.pkl')
+
+ if use_flight_altitude is True:
+ flight_levels = [0, 1, 2, 3, 4]
+ if use_max_cth_level:
+ flight_levels = [-1]
+ 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
+ train_params = ['cld_temp_acha', 'supercooled_cloud_fraction', 'cld_reff_dcomp', 'cld_opd_dcomp']
+
+ if satellite == 'H08':
+ clvrx_ds = CLAVRx_H08(clvrx_dir)
+ elif satellite == 'H09':
+ clvrx_ds = CLAVRx_H09(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, cldmsk, ll, cc = prepare_evaluate1x1(h5f, name_list=train_params, satellite=satellite, domain=domain, offset=8)
+ num_elems = len(cc)
+ num_lines = len(ll)
+
+ day_idxs = solzen < 80.0
+ day_idxs = day_idxs.flatten()
+ num_day_tiles = np.sum(day_idxs)
+
+ nght_idxs = solzen > 100.0
+ nght_idxs = nght_idxs.flatten()
+ 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 = icing_fcn.run_evaluate_static_2(day_model, data_dct, 1,
+ prob_thresh=prob_thresh,
+ flight_levels=flight_levels)
+ for flvl in flight_levels:
+ preds = preds_day_dct[flvl].flatten()
+ probs = probs_day_dct[flvl].flatten()
+ 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 = icing_fcn.run_evaluate_static_2(night_model, data_dct, 1,
+ prob_thresh=prob_thresh,
+ flight_levels=flight_levels)
+ for flvl in flight_levels:
+ preds = preds_nght_dct[flvl].flatten()
+ probs = probs_nght_dct[flvl].flatten()
+ 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()
+
+
--
GitLab