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

snapshot...

parent 43c82de2
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
modules/aeolus/aeolus_amv.py
+ 101
− 16
View file @ eeb6600b
......@@ -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)
......
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