From eeb6600b19fe61ae4154fe85c6a95d1c172c4a44 Mon Sep 17 00:00:00 2001
From: tomrink <rink@ssec.wisc.edu>
Date: Thu, 17 Dec 2020 21:19:30 -0600
Subject: [PATCH] snapshot...
---
modules/aeolus/aeolus_amv.py | 117 ++++++++++++++++++++++++++++++-----
1 file changed, 101 insertions(+), 16 deletions(-)
diff --git a/modules/aeolus/aeolus_amv.py b/modules/aeolus/aeolus_amv.py
index 4a80d503..e809be80 100644
--- a/modules/aeolus/aeolus_amv.py
+++ b/modules/aeolus/aeolus_amv.py
@@ -171,7 +171,8 @@ class OPS(AMVFiles):
'pressure': ('hPa', 'f4'), 'wind_speed': ('m s-1', 'f4'), 'wind_direction': ('degrees', 'f4')}
def get_navigation(self):
- return GEOSNavigation(sub_lon=-75.2)
+ # return GEOSNavigation(sub_lon=-75.2) ?
+ return GEOSNavigation(sub_lon=-75.0)
def get_datetime(self, pathname):
fname = os.path.split(pathname)[1]
@@ -471,19 +472,19 @@ def run_best_fit_all():
amv_prod_list = []
for k, file in enumerate(raob_files):
raob_dct, ts = get_raob_dict_cdf(raob_dir+file)
- # m_d = match_amvs_to_raobs(raob_dct, ts, amv_files=amv_files)
- # amvs_list, bf_list, raob_match_list, bf_gfs_list = run_best_fit(m_d, raob_dct, gfs_dir+gfs_files[k],
- # amvs_list=amvs_list, bf_list=bf_list, raob_match_list=raob_match_list, bf_gfs_list=bf_gfs_list)
- # prd_dct = get_product_at_locs(m_d, ts, prd_files, amv_prod_list=amv_prod_list)
- # analyze2(amvs_list, bf_list, raob_match_list, bf_gfs_list, amv_prod_list)
- #
- amvs = get_amvs(amv_files, ts)
- amvs_list.append(amvs)
- bfs = run_best_fit_gfs(amvs, gfs_dir+gfs_files[k], amv_lat_idx=0, amv_lon_idx=1, amv_prs_idx=4, amv_spd_idx=5, amv_dir_idx=6)
- bf_list.append(bfs)
+ m_d = match_amvs_to_raobs(raob_dct, ts, amv_files=amv_files)
+ amvs_list, bf_list, raob_match_list, bf_gfs_list = run_best_fit(m_d, raob_dct, gfs_dir+gfs_files[k],
+ amvs_list=amvs_list, bf_list=bf_list, raob_match_list=raob_match_list, bf_gfs_list=bf_gfs_list)
+ prd_dct = get_product_at_locs(m_d, ts, prd_files, amv_prod_list=amv_prod_list)
+ analyze2(amvs_list, bf_list, raob_match_list, bf_gfs_list, amv_prod_list)
+
+ # amvs = get_amvs(amv_files, ts)
+ # amvs_list.append(amvs)
+ # bfs = run_best_fit_gfs(amvs, gfs_dir+gfs_files[k], amv_lat_idx=0, amv_lon_idx=1, amv_prs_idx=4, amv_spd_idx=5, amv_dir_idx=6)
+ # bf_list.append(bfs)
# bin_ranges, bin_pres, bin_spd, bin_dir = analyze2(amvs_list, bf_list, raob_match_list, bf_gfs_list, amv_prod_list)
- bin_ranges, bin_pres, bin_spd, bin_dir = compare_amvs_bestfit(amvs_list, bf_list, bin_size=100)
+ bin_ranges, bin_pres, bin_spd, bin_dir = compare_amvs_bestfit_all(amvs_list, bf_list, bin_size=100)
return bin_ranges, bin_pres, bin_spd, bin_dir
@@ -1030,11 +1031,19 @@ def analyze2(amvs_list, bf_list, raob_match_list, bf_gfs_list, amv_prod_list):
return bin_ranges, bin_pres, bin_spd, bin_dir
-def compare_amvs_bestfit(amvs_list, bfs_list, bin_size=200):
+def compare_amvs_bestfit_all(amvs_list, bfs_list, bin_size=200):
amvs = np.concatenate(amvs_list, axis=1)
amvs = np.transpose(amvs, axes=[1, 0])
bfs = np.stack(bfs_list, axis=0)
+ return compare_amvs_bestfit(amvs, bfs, bin_size=bin_size)
+
+
+def compare_amvs_bestfit(amvs, bfs, bin_size=200):
+ # amvs = np.concatenate(amvs_list, axis=1)
+ # amvs = np.transpose(amvs, axes=[1, 0])
+ # bfs = np.stack(bfs_list, axis=0)
+
good_amvs = amvs
num_good = good_amvs.shape[0]
didx = 6
@@ -1150,11 +1159,11 @@ def compare_amvs_bestfit(amvs_list, bfs_list, bin_size=200):
def make_plot():
# f = open('/Users/tomrink/amv_raob.pkl', 'rb')
- f = open('/Users/tomrink/bf_bf_gfs.pkl', 'rb')
+ f = open('/Users/tomrink/amv_bf_gfs_all_linear.pkl', 'rb')
tup_r = pickle.load(f)
f.close()
- f = open('/Users/tomrink/amv_gfs.pkl', 'rb')
+ f = open('/Users/tomrink/amv_bf_gfs_all_nearest.pkl', 'rb')
tup_g = pickle.load(f)
f.close()
@@ -1199,8 +1208,48 @@ def make_plot():
spd_mad_g.append(np.average(np.abs(bin_spd_g[i])))
spd_bias_g.append(np.average(bin_spd_g[i]))
- #do_plot(x_values, [pres_mad_r], ['LBF'], ['blue'], title='RAOB-BFS', x_axis_label='MAD', y_axis_label='hPa', invert=True, flip=True)
+ do_plot(x_values, [pres_mad_r, pres_mad_g], ['GFS_linear', 'GFS_nearest'], ['blue', 'red'], title='ACHA - BestFit', x_axis_label='MAD', y_axis_label='hPa', invert=True, flip=True)
+ #do_plot(x_values, [pres_mad_r, pres_mad_g], ['RAOB', 'GFS'], ['blue', 'red'], title='ACHA - BestFit', x_axis_label='MAD', y_axis_label='hPa', invert=True, flip=True)
+ #do_plot(x_values, [spd_mad_r, spd_mad_g], ['RAOB', 'GFS'], ['blue', 'red'], title='ACHA - BestFit', x_axis_label='MAE (m/s)', y_axis_label='hPa', invert=True, flip=True)
+ #do_plot(x_values, [pres_bias_r, pres_bias_g], ['RAOB', 'GFS'], ['blue', 'red'], title='ACHA - BestFit', x_axis_label='BIAS', y_axis_label='hPa', invert=True, flip=True)
+ #do_plot(x_values, [spd_bias_r, spd_bias_g], ['RAOB', 'GFS'], ['blue', 'red'], title='ACHA - BestFit', x_axis_label='BIAS (m/s)', y_axis_label='hPa', invert=True, flip=True)
+ #do_plot(x_values, [num_pres_r, num_pres_g], ['RAOB:'+str(num_r), 'GFS:'+str(num_g)], ['blue', 'red'], x_axis_label='Normalized Count', y_axis_label='hPa', invert=True, flip=True)
+
+
+def make_plot2(bin_pres_r, bin_ranges):
+
+ x_values = []
+ num_pres_r = []
+ num_pres_g = []
+ num_spd = []
+ num_dir = []
+ pres_mad_r = []
+ pres_bias_r = []
+ pres_mad_g = []
+ pres_bias_g = []
+ spd_mad_r = []
+ spd_bias_r = []
+ spd_mad_g = []
+ spd_bias_g = []
+
+ num_r = 0
+ for i in range(len(bin_ranges)):
+ num_r += bin_pres_r[i].shape[0]
+
+ for i in range(len(bin_ranges)):
+ x_values.append(np.average(bin_ranges[i]))
+ num_pres_r.append((bin_pres_r[i].shape[0])/num_r)
+
+ pres_mad_r.append(np.average(np.abs(bin_pres_r[i])))
+ pres_bias_r.append(np.average(bin_pres_r[i]))
+
+
+ #spd_mad_r.append(np.average(np.abs(bin_spd_r[i])))
+ #spd_bias_r.append(np.average(bin_spd_r[i]))
+
+
+ do_plot(x_values, [pres_mad_r], ['RAOB'], ['blue'], title='ACHA - BestFit', x_axis_label='MAD', y_axis_label='hPa', invert=True, flip=True)
#do_plot(x_values, [pres_mad_r, pres_mad_g], ['RAOB', 'GFS'], ['blue', 'red'], title='ACHA - BestFit', x_axis_label='MAD', y_axis_label='hPa', invert=True, flip=True)
#do_plot(x_values, [spd_mad_r, spd_mad_g], ['RAOB', 'GFS'], ['blue', 'red'], title='ACHA - BestFit', x_axis_label='MAE (m/s)', y_axis_label='hPa', invert=True, flip=True)
#do_plot(x_values, [pres_bias_r, pres_bias_g], ['RAOB', 'GFS'], ['blue', 'red'], title='ACHA - BestFit', x_axis_label='BIAS', y_axis_label='hPa', invert=True, flip=True)
@@ -1208,6 +1257,42 @@ def make_plot():
#do_plot(x_values, [num_pres_r, num_pres_g], ['RAOB:'+str(num_r), 'GFS:'+str(num_g)], ['blue', 'red'], x_axis_label='Normalized Count', y_axis_label='hPa', invert=True, flip=True)
+def make_plot3(amvs_list, bfs_list):
+ num = len(amvs_list)
+
+ bin_ranges = get_press_bin_ranges(100, 1000, bin_size=100)
+
+ y_values = []
+ x_values = None
+ line_labels = []
+ colors = []
+ for k in range(num):
+ bfs = np.stack(bfs_list[k])
+ amvs = np.concatenate(amvs_list[k], axis=1)
+ amvs = np.transpose(amvs, axes=[1, 0])
+ #amvs = amvs_list[k]
+
+ bin_ranges, bin_pres, bin_spd, bin_dir = compare_amvs_bestfit(amvs, bfs, bin_size=100)
+
+ #diff = amvs[:, 4] - amvs[:, 7]
+ #bin_pres = bin_data_by(diff, amvs[:, 4], bin_ranges)
+
+ pres_mad = []
+ x_values = []
+ for i in range(len(bin_ranges)):
+ x_values.append(np.average(bin_ranges[i]))
+ #num_pres_r.append((bin_pres[i].shape[0]) / num_r)
+
+ pres_mad.append(np.average(np.abs(bin_pres[i])))
+ #pres_bias.append(np.average(bin_pres[i]))
+
+ line_labels.append(None)
+ colors.append(None)
+ y_values.append(pres_mad)
+
+ do_plot(x_values, y_values, line_labels, colors, title='Daily ACHA - BestFit RAOB', x_axis_label='MAD', y_axis_label='hPa', invert=True, flip=True)
+
+
# imports the S4 NOAA output
# filename: full path as a string, '/home/user/filename'
# returns a dict: time -> list of profiles (a profile is a list of levels)
--
GitLab