From 80ef39b0d9adea211763e1527dcc5dd89c4ebdaa Mon Sep 17 00:00:00 2001
From: tomrink <rink@ssec.wisc.edu>
Date: Sat, 5 Dec 2020 12:59:28 -0600
Subject: [PATCH] snapshot
---
modules/aeolus/aeolus_amv.py | 42 ++++++++++++++++++++++++++++++++++--
modules/amv/intercompare.py | 23 +++++++++++++++++---
2 files changed, 60 insertions(+), 5 deletions(-)
diff --git a/modules/aeolus/aeolus_amv.py b/modules/aeolus/aeolus_amv.py
index 566e59d9..b97220b8 100644
--- a/modules/aeolus/aeolus_amv.py
+++ b/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])
diff --git a/modules/amv/intercompare.py b/modules/amv/intercompare.py
index b89a32d7..57952df3 100644
--- a/modules/amv/intercompare.py
+++ b/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
--
GitLab