snapshot...

modules/aeolus/aeolus_amv.py

@@ -167,8 +167,8 @@ def match_amvs_to_raobs(raob_dict,  raob_time, amv_files):
     match_dict = {}
 
     #fname, ftime, f_idx = amv_files.get_file_containing_time(raob_time)
-    # fname = '/Users/tomrink/data/OR_ABI-L2-DMWF-M6C14_G16_s20201190000156_e20201190009464_c20201190023107.nc'
-    fname = '/Users/tomrink/data/OR_ABI-L2-DMWF-M6C14_G16_s20201191200158_e20201191209466_c20201191223041.nc'
+    fname = '/Users/tomrink/data/OR_ABI-L2-DMWF-M6C14_G16_s20201190000156_e20201190009464_c20201190023107.nc'
+    # fname = '/Users/tomrink/data/OR_ABI-L2-DMWF-M6C14_G16_s20201191200158_e20201191209466_c20201191223041.nc'
 
     ds = Dataset(fname)
 
@@ -336,6 +336,7 @@ def analyze2(raob_to_amv_dct, raob_dct):
 
     amvs_list = []
     bf_list = []
+    raob_match_list = []
     for key in keys:
         rlat = key[0]
         rlon = key[1]
@@ -361,9 +362,19 @@ def analyze2(raob_to_amv_dct, raob_dct):
                 bspd, bdir = spd_dir_from_uv(bf[0], bf[1])
                 #print(amv_spd, bspd, amv_dir, bdir)
 
+            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)
+
+
     amvs = np.concatenate(amvs_list, axis=1)
     amvs = np.transpose(amvs, axes=[1, 0])
     bfs = np.stack(bf_list, axis=0)
+    raob_match = np.stack(raob_match_list, axis=0)
 
     good_amvs = amvs
     num_good = good_amvs.shape[0]
@@ -430,6 +441,7 @@ def analyze2(raob_to_amv_dct, raob_dct):
     bf_spd, bf_dir = spd_dir_from_uv(bfs[keep_idxs, 0], bfs[keep_idxs, 1])
 
     diff = amv_spd * units('m/s') - bf_spd
+    diff = diff.magnitude
     spd_mad = np.average(np.abs(diff))
     spd_bias = np.average(diff)
     print('spd, MAD: {0:.2f}'.format(spd_mad))
@@ -454,21 +466,23 @@ def analyze2(raob_to_amv_dct, raob_dct):
     vd_mean = np.mean(vd)
     vd_std = np.std(vd)
     print('VD bias/rms: {0:.2f}  {1:.2f}'.format(vd_mean, np.sqrt(vd_mean**2 + vd_std**2)))
-    print('------------------------------------------')
+    print('******************************************************')
 
     x_values = []
     num_pres = []
     num_spd = []
     num_dir = []
+    print('level   num cases   hgt MAD/bias    spd MAD/bias     dir MAD/bias')
+    print('-------------------------------------------------------------------')
     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
+        print('{0:d}        {1:d}       {2:.2f}/{3:.2f}       {4:.2f}/{5:.2f}        {6:.2f}/{7:.2f}'
+              .format(int(x_values[i]), num_pres[i], np.average(np.abs(bin_pres[i])), np.average(bin_pres[i]),
+                      np.average(np.abs(bin_spd[i])), np.average(bin_spd[i]), np.average(np.abs(bin_dir[i])), np.average(bin_dir[i])))
 
 
 # imports the S4 NOAA output