From c350b0bc7b4039d76fb26f931e169909499ca8f9 Mon Sep 17 00:00:00 2001
From: tomrink <rink@ssec.wisc.edu>
Date: Thu, 15 Dec 2022 17:21:33 -0600
Subject: [PATCH] snapshot..
---
modules/icing/pirep_goes.py | 98 +++++++++++++++++++++----------------
1 file changed, 56 insertions(+), 42 deletions(-)
diff --git a/modules/icing/pirep_goes.py b/modules/icing/pirep_goes.py
index 85589207..6e90718d 100644
--- a/modules/icing/pirep_goes.py
+++ b/modules/icing/pirep_goes.py
@@ -906,31 +906,48 @@ def process_boeing(ice_dct, no_ice_dct):
def analyze2(filename, filename_l1b):
f = h5py.File(filename, 'r')
- icing_alt = f['icing_altitude'][:]
- cld_top_hgt = f['cld_height_acha'][:, 12:28, 12:28]
- cld_phase = f['cloud_phase'][:, 12:28, 12:28]
- cld_opd_dc = f['cld_opd_dcomp'][:, 12:28, 12:28]
- cld_opd = f['cld_opd_acha'][:, 12:28, 12:28]
- solzen = f['solar_zenith_angle'][:, 12:28, 12:28]
-
- f_l1b = h5py.File(filename_l1b, 'r')
- bt_11um = f_l1b['temp_11_0um_nom'][:, 12:28, 12:28]
-
- bt_11um = bt_11um.flatten()
- cld_opd = cld_opd.flatten()
- cld_opd_dc = cld_opd_dc.flatten()
- solzen = solzen.flatten()
-
- keep1 = np.invert(np.isnan(cld_opd))
- keep2 = np.invert(np.isnan(solzen))
- keep = keep1 & keep2
- cld_opd = cld_opd[np.invert(np.isnan(cld_opd))]
- cld_opd_dc = cld_opd_dc[keep]
- solzen = solzen[keep]
-
- bt_11um = bt_11um[np.invert(np.isnan(bt_11um))]
+ iint = f['icing_intensity'][:]
+ icing_alt = f['flight_altitude'][:]
+ cld_top_hgt = f['cld_height_acha'][:, :, :]
+ cld_top_tmp = f['cld_temp_acha'][:, :, :]
+ sc_cld_frac = f['supercooled_cloud_fraction'][:, :, :]
+ # cld_phase = f['cloud_phase'][:, :, :]
+ # cld_opd_dc = f['cld_opd_dcomp'][:, :, :]
+ # cld_opd = f['cld_opd_acha'][:, :, :]
+ # solzen = f['solar_zenith_angle'][:, :, :]
+
+ cth = np.reshape(cld_top_hgt, [iint.size, -1])
+ ctt = np.reshape(cld_top_tmp, [iint.size, -1])
+ scf = np.reshape(sc_cld_frac, [iint.size, -1])
+ cth = np.nanmedian(cth, axis=1)
+ ctt = np.nanmedian(ctt, axis=1)
+ scf = np.nanmedian(scf, axis=1)
+ #cth = np.nanmax(cth, axis=1)
+ #ctt = np.nanmax(ctt, axis=1)
+ #cth = np.nanmin(cth, axis=1)
+ #ctt = np.nanmin(ctt, axis=1)
+ print(cth.shape, ctt.shape)
+
+ keep_0 = icing_alt < 6000.0
+ keep_1 = cth < 6000.0
+ keep_2 = ctt < 274.5
+ keep_3 = np.logical_and(icing_alt > cth - 1000.0, icing_alt < cth + 1000.0)
+ # keep = keep_0 & keep_1 & keep_2 & keep_3
+ keep = keep_0 & keep_1 & keep_2
+ icing_alt = icing_alt[keep]
+ iint = iint[keep]
+ cth = cth[keep]
+ ctt = ctt[keep]
+ scf = scf[keep]
+ print(iint.shape, cth.shape)
+ iint = np.where(iint == -1, 0, iint)
+ iint = np.where(iint != 0, 1, iint)
+ no_ice = iint == 0
+ ice = iint == 1
+ print('num no ice: ', np.sum(no_ice))
+ print('num ice: ', np.sum(ice))
- return cld_opd, bt_11um
+ return icing_alt[no_ice], icing_alt[ice], cth[no_ice], cth[ice], ctt[no_ice], ctt[ice], scf[no_ice], scf[ice]
# --------------------------------------------
@@ -1213,8 +1230,9 @@ def apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_opd
opd_thin_threshold = 0.1
closeness_top = 100.0 # meters
- max_depth = 3000.0
+ max_cld_depth = 3000.0
max_altitude = 4000.0
+ max_cld_altitude = 6000.0
num_obs = len(icing_alt)
cld_mask = cld_mask.reshape((num_obs, -1))
@@ -1251,11 +1269,11 @@ def apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_opd
# Test 1
keep = np.where(keep, (cld_top_hgt[i,] + closeness_top) > icing_alt[i], False)
- # keep = np.where(keep, cld_top_hgt[i,] < max_altitude, False)
- # keep = np.where(keep, (cld_top_hgt[i,] - max_depth) < icing_alt[i], False)
+ # keep = np.where(keep, cld_top_hgt[i,] < max_cld_altitude, False)
+ # keep = np.where(keep, (cld_top_hgt[i,] - max_cld_depth) < icing_alt[i], False)
# keep = np.where(keep, np.logical_and(bt_11um[i,] > 228.0, bt_11um[i,] < 270.0), False)
- keep = np.where(keep, bt_11um[i,] < 276.0, False)
+ keep = np.where(keep, bt_11um[i,] < 275.0, False)
mask.append(keep)
idxs.append(i)
@@ -1274,8 +1292,9 @@ def apply_qc_no_icing_pireps(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_
opd_thin_threshold = 0.1
closeness_top = 100.0 # meters
- max_depth = 3000.0
+ max_cld_depth = 3000.0
max_altitude = 4000.0
+ max_cld_altitude = 6000.0
num_obs = len(icing_alt)
cld_mask = cld_mask.reshape((num_obs, -1))
@@ -1315,7 +1334,7 @@ def apply_qc_no_icing_pireps(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_
# keep = np.where(keep, (cld_top_hgt[i,] - max_depth) < icing_alt[i], False)
# keep = np.where(keep, np.logical_and(bt_11um[i,] > 228.0, bt_11um[i,] < 270.0), False)
- keep = np.where(keep, bt_11um[i,] < 276.0, False)
+ keep = np.where(keep, bt_11um[i,] < 275.0, False)
mask.append(keep)
idxs.append(i)
@@ -1590,7 +1609,6 @@ def tile_extract(icing_files, no_icing_files, trnfile='/home/rink/tiles_train.h5
ice_lon_s = np.array(ice_lon_s)
ice_lat_s = np.array(ice_lat_s)
ice_flt_alt_s = np.array(ice_flt_alt_s)
- num_ice = icing_int_s.shape[0]
# No icing ------------------------------------------------------------
num_no_ice = 0
@@ -1636,17 +1654,9 @@ def tile_extract(icing_files, no_icing_files, trnfile='/home/rink/tiles_train.h5
icing_lats = np.concatenate([ice_lat_s, no_ice_lat_s])
icing_alt = np.concatenate([ice_flt_alt_s, no_ice_flt_alt_s])
- # icing_intensity = icing_int_s
- # icing_times = ice_time_s
- # icing_lons = ice_lon_s
- # icing_lats = ice_lat_s
- # icing_alt = ice_flt_alt_s
-
data_dct = {}
for ds_name in params:
data_dct[ds_name] = np.concatenate([icing_data_dct[ds_name], no_icing_data_dct[ds_name]])
- # for ds_name in params:
- # data_dct[ds_name] = icing_data_dct[ds_name]
# do sort -------------------------------------
ds_indexes = np.argsort(icing_times)
@@ -1660,6 +1670,8 @@ def tile_extract(icing_files, no_icing_files, trnfile='/home/rink/tiles_train.h5
if do_split:
trn_idxs, vld_idxs, tst_idxs = split_data(icing_times)
+ # Below for now test data, just train and valid
+ # trn_idxs = np.concatenate([trn_idxs, tst_idxs])
else:
trn_idxs = np.arange(icing_intensity.shape[0])
tst_idxs = None
@@ -1745,6 +1757,7 @@ def tile_extract(icing_files, no_icing_files, trnfile='/home/rink/tiles_train.h5
write_file(trnfile, params, param_types, trn_data_dct, trn_icing_intensity, trn_icing_times, trn_icing_lons, trn_icing_lats, trn_icing_alt)
if do_split:
+ # --- Comment out this block for not test data, just valid and train
tst_data_dct = {}
for ds_name in params:
tst_data_dct[ds_name] = data_dct[ds_name][tst_idxs,]
@@ -1765,6 +1778,7 @@ def tile_extract(icing_files, no_icing_files, trnfile='/home/rink/tiles_train.h5
tst_icing_alt = tst_icing_alt[ds_indexes]
write_file(tstfile, params, param_types, tst_data_dct, tst_icing_intensity, tst_icing_times, tst_icing_lons, tst_icing_lats, tst_icing_alt)
+ # -------------------------------------------------------------------
vld_data_dct = {}
for ds_name in params:
@@ -2537,8 +2551,8 @@ def analyze(preds_file, labels, prob_avg, test_file):
iint = np.where(iint == -1, 0, iint)
iint = np.where(iint != 0, 1, iint)
- nda[np.logical_and(nda >= 0, nda < 4000)] = 0
- nda[np.logical_and(nda >= 4000, nda < 15000)] = 1
+ nda[np.logical_and(nda >= 0, nda < 2500)] = 0
+ nda[np.logical_and(nda >= 2500, nda < 4000)] = 1
# nda[np.logical_and(nda >= 8000, nda < 15000)] = 2
#print(np.sum(nda == 0), np.sum(nda == 1), np.sum(nda == 2))
@@ -2567,7 +2581,7 @@ def analyze(preds_file, labels, prob_avg, test_file):
print('----------------------------------------------------------')
print('----------------------------------------------------------')
- return flt_alt[iint == 0], flt_alt[iint == 1]
+ # return flt_alt[iint == 0], flt_alt[iint == 1]
if prob_avg is None:
return
--
GitLab