From e2d0bb6331e8518f6e0347f2b734ca397d8ad907 Mon Sep 17 00:00:00 2001 From: tomrink <rink@ssec.wisc.edu> Date: Mon, 4 Jan 2021 12:10:46 -0600 Subject: [PATCH] snapshot... --- modules/aeolus/aeolus_amv.py | 95 ++++++++++++++++++++++++++++++++++++ 1 file changed, 95 insertions(+) diff --git a/modules/aeolus/aeolus_amv.py b/modules/aeolus/aeolus_amv.py index 3b066ead..183c9333 100644 --- a/modules/aeolus/aeolus_amv.py +++ b/modules/aeolus/aeolus_amv.py @@ -489,6 +489,101 @@ def run_best_fit(raob_to_amv_dct, raob_dct, gfs_filename=None): return out_dct +def compare_gfs_raob(raob_dct, gfs_filename): + keys = list(raob_dct.keys()) + + gfs_press = None + gfs_spd = None + gfs_dir = None + + locs = np.array(keys) + do_gfs_best_fit = True + xr_dataset = xr.open_dataset(gfs_filename) + gfs_press = xr_dataset['pressure levels'] + gfs_press = gfs_press.values + gfs_press = gfs_press[::-1] + g_0 = 5 # 850 + g_1 = 8 # 700 + g_2 = 12 # 500 + g_3 = 16 # 300 + g_4 = 20 # 100 + + lon_idx = LatLonTuple._fields.index('lon') + lat_idx = LatLonTuple._fields.index('lat') + + uv_wind = get_vert_profile_s(xr_dataset, ['u-wind', 'v-wind'], locs[:, lon_idx], locs[:, lat_idx], method='nearest') + uv_wind = uv_wind.values + wspd, wdir = spd_dir_from_uv(uv_wind[0, :, :], uv_wind[1, :, :]) + wspd = wspd.magnitude + wdir = wdir.magnitude + gfs_spd = wspd[:, ::-1] + gfs_dir = wdir[:, ::-1] + gfs_u = uv_wind[0, :, ::-1] + gfs_v = uv_wind[1, :, ::-1] + + raob_intrp = [] + + lvl_0 = [] + lvl_1 = [] + lvl_2 = [] + lvl_3 = [] + lvl_4 = [] + + for key_idx, key in enumerate(keys): + raob = raob_dct.get(key) + raob_prs = raob[:, 0] + raob_spd = raob[:, 2] + raob_dir = raob[:, 3] + raob_u, raob_v = uv_from_spd_dir(raob_spd, raob_dir) + raob_u = raob_u.magnitude + raob_v = raob_v.magnitude + + r_0 = np.where(raob_prs == 850)[0] + r_1 = np.where(raob_prs == 700)[0] + r_2 = np.where(raob_prs == 500)[0] + r_3 = np.where(raob_prs == 300)[0] + r_4 = np.where(raob_prs == 100)[0] + + spd_intrp = np.interp(gfs_press, raob_prs, raob_v) + #print(raob_v[r_0] - gfs_v[key_idx, g_0], raob_v[r_1] - gfs_v[key_idx, g_1], raob_v[r_3] - gfs_v[key_idx, g_3]) + + if len(r_0) > 0: + r_0 = r_0[0] + lvl_0.append(raob_v[r_0] - gfs_v[key_idx, g_0]) + if len(r_1) > 0: + r_1 = r_1[0] + lvl_1.append(raob_v[r_1] - gfs_v[key_idx, g_1]) + if len(r_2) > 0: + r_2 = r_2[0] + lvl_2.append(raob_v[r_2] - gfs_v[key_idx, g_2]) + if len(r_3) > 0: + r_3 = r_3[0] + lvl_3.append(raob_v[r_3] - gfs_v[key_idx, g_3]) + if len(r_4) > 0: + r_4 = r_4[0] + lvl_4.append(raob_v[r_4] - gfs_v[key_idx, g_4]) + + raob_intrp.append(spd_intrp) + + raob_v_intrp = np.stack(raob_intrp) + + v_diff = np.abs(raob_v_intrp - gfs_v) + print(np.average(v_diff, axis=0)) + print('---------------') + lvl_0 = np.array(lvl_0) + lvl_1 = np.array(lvl_1) + lvl_2 = np.array(lvl_2) + lvl_3 = np.array(lvl_3) + lvl_4 = np.array(lvl_4) + print(np.average(np.abs(lvl_0))) + print(np.average(np.abs(lvl_1))) + print(np.average(np.abs(lvl_2))) + print(np.average(np.abs(lvl_3))) + print(np.average(np.abs(lvl_4))) + + return None + + def analyze2(amvs, bfs, raob_match, bfs_gfs, amv_prod): good_amvs = amvs -- GitLab