snapshot...

modules/aeolus/aeolus_amv.py

@@ -6,12 +6,13 @@ 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
+from amv.intercompare import best_fit, bin_data_by, get_press_bin_ranges, spd_dir_from_uv, uv_from_spd_dir, direction_difference
 import math
+from metpy.units import units
 
 
 amv_file_duration = 60  # minutes
-half_width = 20  # search box centered on AEOLUS profile (FGF coordinates)
+half_width = 25  # search box centered on AEOLUS profile (FGF coordinates)
 num_elems = 5424
 num_lines = 5424
 
@@ -87,8 +88,8 @@ class OPS(AMVFiles):
         self.lon_name = 'lon'
         self.lat_name = 'lat'
 
-        self.out_params = {'Lon', 'Lat', 'Element', 'Line', 'pressure', 'wind_speed', 'wind_direction'}
-        self.params = {'pressure', 'wind_speed', 'wind_direction'}
+        self.out_params = ['Lon', 'Lat', 'Element', 'Line', 'pressure', 'wind_speed', 'wind_direction']
+        self.params = ['pressure', 'wind_speed', 'wind_direction']
         self.meta_dict = {'Lon': ('degrees east', 'f4'), 'Lat': ('degrees north', 'f4'), 'Element': (None, 'i4'), 'Line': (None, 'i4'),
                           'pressure': ('hPa', 'f4'), 'wind_speed': ('m s-1', 'f4'), 'wind_direction': ('degrees', 'f4')}
 
@@ -171,8 +172,8 @@ def match_amvs_to_raobs(raob_dict,  raob_time, amv_files):
 
     ds = Dataset(fname)
 
-    amv_lons = ds[amv_files.lon_name][:]
-    amv_lats = ds[amv_files.lat_name][:]
+    amv_lons = ds[amv_files.lon_name][:].data
+    amv_lats = ds[amv_files.lat_name][:].data
     if amv_files.elem_name is None or amv_files.line_name is None:
         cc, ll = nav.earth_to_lc_s(amv_lons, amv_lats)
     else:
@@ -184,12 +185,24 @@ def match_amvs_to_raobs(raob_dict,  raob_time, amv_files):
     param_s.append(amv_lats)
     param_s.append(cc)
     param_s.append(ll)
+
+    vld = None
     for param in amv_params:
+        data = ds[param][:]
+        if vld is None:
+            vld = np.invert(data.mask)
+        else:
+            vld = np.logical_and(vld, np.invert(data.mask))
         if param == 'V_3D':
-            param_s.append(ds[param][:, 0])
-            param_s.append(ds[param][:, 1])
+            param_s.append(data.data[0])
+            param_s.append(data.data[1])
         else:
-            param_s.append(ds[param][:])
+            param_s.append(data.data)
+
+    param_nd = np.vstack(param_s)
+    param_nd = param_nd[:, vld]
+    cc = param_nd[2, :]
+    ll = param_nd[3, :]
 
     ds.close()
 
@@ -215,10 +228,7 @@ def match_amvs_to_raobs(raob_dict,  raob_time, amv_files):
         num_amvs = np.sum(in_box)
         if num_amvs == 0:
             continue
-        # dist = haversine_np(lon, lat, amv_lons[in_box], amv_lats[in_box])
-        param_nd = np.vstack(param_s)
-        param_nd = param_nd[:, in_box]
-        match_dict[key] = param_nd
+        match_dict[key] = param_nd[:, in_box]
 
     return match_dict
 
@@ -317,9 +327,15 @@ def create_file2(filename, raob_to_amv_dct, raob_dct, amv_files):
     rootgrp.close()
 
 
+def bulk_stats(filename):
+    pass
+
+
 def analyze2(raob_to_amv_dct, raob_dct):
     keys = list(raob_to_amv_dct.keys())
 
+    amvs_list = []
+    bf_list = []
     for key in keys:
         rlat = key[0]
         rlon = key[1]
@@ -330,6 +346,7 @@ def analyze2(raob_to_amv_dct, raob_dct):
         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]
@@ -339,9 +356,118 @@ def analyze2(raob_to_amv_dct, raob_dct):
             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)
             if bf[3] == 0:
-                print(nlevs, amv_prs, bf[2])
-
+                bspd, bdir = spd_dir_from_uv(bf[0], bf[1])
+                #print(amv_spd, bspd, amv_dir, bdir)
+
+    amvs = np.concatenate(amvs_list, axis=1)
+    amvs = np.transpose(amvs, axes=[1, 0])
+    bfs = np.stack(bf_list, axis=0)
+
+    good_amvs = amvs
+    num_good = good_amvs.shape[0]
+    didx = 6
+    sidx = 5
+    pidx = 4
+
+    spd_min = good_amvs[:, sidx].min()
+    spd_max = good_amvs[:, sidx].max()
+    print('spd min/max/mean: ', spd_min, spd_max, np.average(good_amvs[:, sidx]))
+
+    p_min = good_amvs[:, pidx].min()
+    p_max = good_amvs[:, pidx].max()
+    print('pres min/max/mean: ', p_min, p_max, np.average(good_amvs[:, pidx]))
+
+    low = good_amvs[:, pidx] >= 700
+    mid = np.logical_and(good_amvs[:, pidx] < 700, good_amvs[:, pidx] > 400)
+    hgh = good_amvs[:, pidx] <= 400
+
+    n_low = np.sum(low)
+    n_mid = np.sum(mid)
+    n_hgh = np.sum(hgh)
+
+    print('% low: ', 100.0*(n_low/num_good))
+    print('% mid: ', 100.0*(n_mid/num_good))
+    print('% hgh: ', 100.0*(n_hgh/num_good))
+    print('---------------------------')
+
+    print('Low Spd min/max/mean: ', good_amvs[low, sidx].min(), good_amvs[low, sidx].max(), good_amvs[low,sidx].mean())
+    print('Low Press min/max/mean: ', good_amvs[low, pidx].min(), good_amvs[low, pidx].max(), good_amvs[low, pidx].mean())
+    print('---------------------------')
+
+    print('Mid Spd min/max/mean: ', good_amvs[mid, sidx].min(), good_amvs[mid, sidx].max(), good_amvs[mid, sidx].mean())
+    print('Mid Press min/max/mean: ', good_amvs[mid, pidx].min(), good_amvs[mid, pidx].max(), good_amvs[mid, pidx].mean())
+    print('---------------------------')
+
+    print('Hgh Spd min/max/mean: ', good_amvs[hgh, sidx].min(), good_amvs[hgh, sidx].max(), good_amvs[hgh, sidx].mean())
+    print('Hgh Press min/max/mean: ', good_amvs[hgh, pidx].min(), good_amvs[hgh, pidx].max(), good_amvs[hgh, pidx].mean())
+
+    bin_size = 200.0
+    vld_bf = bfs[:, 3] == 0
+    keep_idxs = vld_bf
+
+    amv_p = good_amvs[keep_idxs, pidx]
+    bf_p = bfs[keep_idxs, 2]
+    diff = amv_p - bf_p
+    mad = np.average(np.abs(diff))
+    bias = np.average(diff)
+    print('********************************************************')
+    print('num of best fits: ', bf_p.shape[0])
+    print('press, MAD: ', mad)
+    print('press, bias: ', bias)
+    pd_std = np.std(diff)
+    pd_mean = np.mean(diff)
+    print('press bias/rms ', pd_mean, np.sqrt(pd_mean**2 + pd_std**2))
+    print('------------------------------------------')
+
+    bin_ranges = get_press_bin_ranges(50, 1050, bin_size=bin_size)
+    bin_pres = bin_data_by(diff, amv_p, bin_ranges)
+
+    amv_spd = good_amvs[keep_idxs, sidx]
+    amv_dir = good_amvs[keep_idxs, didx]
+    bf_spd, bf_dir = spd_dir_from_uv(bfs[keep_idxs, 0], bfs[keep_idxs, 1])
+
+    diff = amv_spd * units('m/s') - bf_spd
+    spd_mad = np.average(np.abs(diff))
+    spd_bias = np.average(diff)
+    print('spd, MAD: ', spd_mad)
+    print('spd, bias: ', spd_bias)
+    spd_mean = np.mean(diff)
+    spd_std = np.std(diff)
+    print('spd MAD/bias/rms: ', np.average(np.abs(diff)), spd_mean, np.sqrt(spd_mean**2 + spd_std**2))
+    print('-----------------')
+    bin_spd = bin_data_by(diff, amv_p, bin_ranges)
+
+    dir_diff = direction_difference(amv_dir, bf_dir.magnitude)
+    print('dir, MAD: ', np.average(np.abs(dir_diff)))
+    print('dir bias: ', np.average(dir_diff))
+    print('-------------------------------------')
+    bin_dir = bin_data_by(dir_diff, amv_p, bin_ranges)
+
+    amv_u, amv_v = uv_from_spd_dir(good_amvs[keep_idxs, sidx], good_amvs[keep_idxs, didx])
+    u_diffs = amv_u - (bfs[keep_idxs, 0] * units('m/s'))
+    v_diffs = amv_v - (bfs[keep_idxs, 1] * units('m/s'))
+
+    vd = np.sqrt(u_diffs**2 + v_diffs**2)
+    vd_mean = np.mean(vd)
+    vd_std = np.std(vd)
+    print('VD bias/rms: ', vd_mean, np.sqrt(vd_mean**2 + vd_std**2))
+    print('------------------------------------------')
+
+    x_values = []
+    num_pres = []
+    num_spd = []
+    num_dir = []
+    for i in range(len(bin_ranges)):
+        x_values.append(np.average(bin_ranges[i]))
+        num_pres.append(bin_pres[i].shape[0])
+        num_spd.append(bin_spd[i].shape[0])
+        num_dir.append(bin_dir[i].shape[0])
+
+    #return x_values, bin_pres, num_pres, bin_spd, num_spd, bin_dir, num_dir
+
+    return amvs, bfs
 
 
 # imports the S4 NOAA output