snapshot...

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