Skip to content
Snippets Groups Projects
Commit c007a8b2 authored by tomrink's avatar tomrink
Browse files

snapshot...

parent 61038eee
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 201 additions and 0 deletions
modules/icing/util.py
+ 201
− 0
View file @ c007a8b2
......@@ -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()
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment