snapshot...

modules/aeolus/aeolus_amv.py

@@ -11,6 +11,7 @@ from amv.intercompare import best_fit, bin_data_by, get_press_bin_ranges, spd_di
 import math
 from metpy.units import units
 from util.gfs_reader import get_vert_profile_s
+from amv.intercompare import get_raob_dict_cdf
 
 
 amv_file_duration = 60  # minutes
@@ -387,6 +388,101 @@ def create_file2(filename, raob_to_amv_dct, raob_dct, amv_files):
     rootgrp.close()
 
 
+def run_best_fit_all():
+    amv_files = Framework('/Users/tomrink/data/amv_bust/', 10)
+
+    raob_dir = '/Users/tomrink/data/amv_bust/raob/'
+
+    raob_files = ['raob_soundings20191117_0000.cdf',
+                  'raob_soundings20191118_0000.cdf',
+                  'raob_soundings20191119_0000.cdf',
+                  'raob_soundings20191120_0000.cdf',
+                  'raob_soundings20191121_0000.cdf',
+                  'raob_soundings20191122_0000.cdf']
+
+    raob_dtsts = ['2019-11-17_00:00',
+                  '2019-11-18_00:00',
+                  '2019-11-19_00:00',
+                  '2019-11-20_00:00',
+                  '2019-11-21_00:00',
+                  '2019-11-22_00:00']
+
+    gfs_dir = None
+    gfs_files = ['gfs.19111612_F012.h5',
+                 'gfs.19111712_F012.h5',
+                 'gfs.19111812_F012.h5',
+                 'gfs.19111912_F012.h5',
+                 'gfs.19112012_F012.h5',
+                 'gfs.19112112_F012.h5']
+
+    for k, file in enumerate(raob_files):
+        raob_dct, ts = get_raob_dict_cdf(raob_dir+file, raob_dtsts[k])
+        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])
+
+
+def run_best_fit(raob_to_amv_dct, raob_dct, gfs_filename=None):
+    keys = list(raob_to_amv_dct.keys())
+
+    do_gfs_best_fit = False
+    gfs_press = None
+    gfs_spd = None
+    gfs_dir = None
+    bf_gfs_list = None
+    if gfs_filename is not None:
+        locs = np.array(keys)
+        do_gfs_best_fit = True
+        xr_dataset = xr.open_dataset(gfs_filename)
+        bf_gfs_list = []
+        gfs_press = xr_dataset['pressure levels']
+        gfs_press = gfs_press.values
+        gfs_press = gfs_press[::-1]
+
+        uv_wind = get_vert_profile_s(xr_dataset, ['u-wind', 'v-wind'], locs[:, 1], locs[:, 0], 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]
+
+    amvs_list = []
+    bf_list = []
+    raob_match_list = []
+    for key_idx, key in enumerate(keys):
+        raob = raob_dct.get(key)
+        nlevs = raob.shape[0]
+        raob_prs = raob[:, 0]
+        raob_spd = raob[:, 2]
+        raob_dir = raob[:, 3]
+        amvs = raob_to_amv_dct.get(key)
+        amvs_list.append(amvs)
+        num_amvs = amvs.shape[1]
+        for i in range(num_amvs):
+            amv_lon = amvs[0, i]
+            amv_lat = amvs[1, i]
+            amv_prs = amvs[4, i]
+            amv_spd = amvs[5, i]
+            amv_dir = amvs[6, i]
+
+            bf = best_fit(amv_spd, amv_dir, amv_prs, amv_lat, amv_lon, raob_spd, raob_dir, raob_prs)
+            bf_list.append(bf)
+
+            pdiff = amv_prs - raob_prs
+            lev_idx = np.argmin(np.abs(pdiff))
+            if np.abs(pdiff[lev_idx]) > 100.0:
+                tup = (raob_spd[lev_idx], raob_dir[lev_idx], raob_prs[lev_idx], -9)
+            else:
+                tup = (raob_spd[lev_idx], raob_dir[lev_idx], raob_prs[lev_idx], 0)
+            raob_match_list.append(tup)
+
+            if do_gfs_best_fit:
+                bf = best_fit(amv_spd, amv_dir, amv_prs, amv_lat, amv_lon, gfs_spd[key_idx], gfs_dir[key_idx], gfs_press)
+                bf_gfs_list.append(bf)
+
+    return amvs_list, bf_list, raob_match_list, bf_gfs_list
+
+
 def analyze2(raob_to_amv_dct, raob_dct, gfs_filename=None):
     keys = list(raob_to_amv_dct.keys())