snapshot

modules/aeolus/aeolus_amv.py

@@ -6,7 +6,8 @@ import xarray as xr
 from netCDF4 import Dataset, Dimension, Variable
 from aeolus.geos_nav import GEOSNavigation
 from util.util import haversine_np
-from amv.intercompare import best_fit, bin_data_by, get_press_bin_ranges, spd_dir_from_uv, uv_from_spd_dir, direction_difference
+from amv.intercompare import best_fit, bin_data_by, get_press_bin_ranges, spd_dir_from_uv, uv_from_spd_dir, \
+    direction_difference, run_best_fit_gfs2
 import math
 from metpy.units import units
 
@@ -26,6 +27,8 @@ class AMVFiles:
 
     def __init__(self, files_path, file_time_span, pattern, band='14'):
         self.flist = glob.glob(files_path + pattern)
+        if len(self.flist) == 0:
+            raise MyGenericException('no matching files found in: ' + files_path)
         self.band = band
         self.ftimes = []
         self.span_seconds = datetime.timedelta(minutes=file_time_span).seconds
@@ -172,6 +175,30 @@ class CarrStereo(AMVFiles):
         return self.meta_dict
 
 
+def get_amvs(amv_files, timestamp, filepath=None):
+    if filepath is None:
+        filepath, ftime, f_idx = amv_files.get_file(timestamp)
+    # amv_params = amv_files.get_parameters()
+    amv_params = ['lon', 'lat', 'pressure', 'wind_speed', 'wind_direction']
+    ds = Dataset(filepath)
+
+    param_s = []
+    vld = None
+    for name in amv_params:
+        param = ds[name][:]
+        if vld is None:
+            vld = np.invert(param.mask)
+        else:
+            vld = np.logical_and(vld, np.invert(param.mask))
+        param_s.append(param.data)
+
+    param_nd = np.vstack(param_s)
+    param_nd = param_nd[:, vld]
+    param_nd = np.transpose(param_nd, axes=[1, 0])
+
+    return param_nd
+
+
 # raob_dict: time -> profiles
 # amv_files_path: directory containing AMVs, '/home/user/amvdir/'
 # return dict: raob -> tuple (amv_lon, amv_lat, amv_pres, amv_spd, amv_dir)
@@ -340,9 +367,17 @@ def create_file2(filename, raob_to_amv_dct, raob_dct, amv_files):
     rootgrp.close()
 
 
-def analyze2(raob_to_amv_dct, raob_dct):
+def analyze2(raob_to_amv_dct, raob_dct, gfs_filename=None):
     keys = list(raob_to_amv_dct.keys())
 
+    do_gfs_best_fit = False
+    if gfs_filename is not None:
+        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
+
     amvs_list = []
     bf_list = []
     raob_match_list = []
@@ -379,6 +414,9 @@ def analyze2(raob_to_amv_dct, raob_dct):
                 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 = run_best_fit_gfs2(xr_dataset, gfs_press, amv_lat, amv_lon, amv_prs, amv_spd, amv_dir)
+                bf_gfs_list.append(bf)
 
     amvs = np.concatenate(amvs_list, axis=1)
     amvs = np.transpose(amvs, axes=[1, 0])

modules/amv/intercompare.py

@@ -3,7 +3,8 @@ import xarray as xr
 from util.setup import home_dir
 from util.lon_lat_grid import earth_to_indexs
 from util.util import haversine_np
-from deeplearning.amv_raob import get_profile_multi
+#from deeplearning.amv_raob import get_profile_multi
+from util.gfs_reader import get_vert_profile_s
 import cartopy
 from cartopy import *
 import metpy.calc as mpcalc
@@ -575,14 +576,15 @@ def run_best_fit_gfs_all(amvs, gfs_filename):
     return bfs_s
 
 
-def run_best_fit_gfs(amvs, gfs_filename):
+def run_best_fit_gfs(amvs, gfs_filename, amv_lat_idx=amv_lat_idx, amv_lon_idx=amv_lon_idx, amv_prs_idx=amv_prs_idx,
+                     amv_spd_idx=amv_spd_idx, amv_dir_idx=amv_dir_idx):
     xr_dataset = xr.open_dataset(gfs_filename)
     num_amvs = amvs.shape[0]
     bestfits = []
 
     gfs_press = xr_dataset['pressure levels']
 
-    uv_wind = get_profile_multi(xr_dataset, ['u-wind', 'v-wind'], amvs[:, amv_lon_idx], amvs[:, amv_lat_idx])
+    uv_wind = get_vert_profile_s(xr_dataset, ['u-wind', 'v-wind'], amvs[:, amv_lon_idx], amvs[:, amv_lat_idx])
     uv_wind = uv_wind.data
     wspd, wdir = spd_dir_from_uv(uv_wind[0,:,:], uv_wind[1,:,:])
     wspd = wspd.magnitude
@@ -606,6 +608,21 @@ def run_best_fit_gfs(amvs, gfs_filename):
 
     return bestfits
 
+def run_best_fit_gfs2(xr_dataset, gfs_press, amv_lat, amv_lon, amv_prs, amv_spd, amv_dir):
+    uv_wind = get_vert_profile_s(xr_dataset, ['u-wind', 'v-wind'], np.array([amv_lon]), np.array([amv_lat]))
+    uv_wind = uv_wind.data
+    wspd, wdir = spd_dir_from_uv(uv_wind[0,:,:], uv_wind[1,:,:])
+
+    rspd = wspd.magnitude
+    rdir = wdir.magnitude
+
+    rspd = rspd.flatten()
+    rdir = rdir.flatten()
+
+    bf = best_fit(amv_spd, amv_dir, amv_prs, amv_lat, amv_lon, rspd, rdir, gfs_press)
+    bf = (bf[0], bf[1], bf[2], bf[3])
+
+    return bf
 
 def raob_gfs_diff(raob_dct, gfs_filename):
     r_press = 900.0