snapshot...

modules/amv/intercompare.py

@@ -69,6 +69,8 @@ amv_qi_idx = amv_hdr_list.index('qi')
 amv_cqi_idx = amv_hdr_list.index('cqi')
 amv_qif_idx = amv_hdr_list.index('qif')
 
+amv_centers_list = ['EUM', 'BRZ', 'JMA', 'KMA', 'NOA', 'NWC']
+
 # AMV inter-compare stuff --------------------------------------------------------
 
 
@@ -731,12 +733,12 @@ def direction_difference(dir_a, dir_b):
 
 def get_amv_winds_match(dist_threshold=150, qitype=amv_cqi_idx, qival=50, lat_range=[-60, 60]):
     #----- match all to EUM
-    amvs_eum = get_amv_nd('/Users/rink/data/amv_intercompare/EUM321_csv.csv')
-    amvs_brz = get_amv_nd('/Users/rink/data/amv_intercompare/BRZCPTECfin_121_csv.csv')
-    amvs_jma = get_amv_nd('/Users/rink/data/amv_intercompare/JMA421_csv.csv')
-    amvs_kma = get_amv_nd('/Users/rink/data/amv_intercompare/KMA521NI_csv.csv')
-    amvs_noa = get_amv_nd('/Users/rink/data/amv_intercompare/NOA621_csv.csv')
-    amvs_nwc = get_amv_nd('/Users/rink/data/amv_intercompare/NWC721_csv.csv')
+    amvs_eum = get_amv_nd('/Users/tomrink/data/amv_intercompare/EUM321_csv.csv')
+    amvs_brz = get_amv_nd('/Users/tomrink/data/amv_intercompare/BRZCPTECfin_121_csv.csv')
+    amvs_jma = get_amv_nd('/Users/tomrink/data/amv_intercompare/JMA421_csv.csv')
+    amvs_kma = get_amv_nd('/Users/tomrink/data/amv_intercompare/KMA521NI_csv.csv')
+    amvs_noa = get_amv_nd('/Users/tomrink/data/amv_intercompare/NOA621_csv.csv')
+    amvs_nwc = get_amv_nd('/Users/tomrink/data/amv_intercompare/NWC721_csv.csv')
 
     amvs_eum = filter_amvs(amvs_eum, qitype=qitype, qival=qival, lat_range=lat_range)
     amvs_brz = filter_amvs(amvs_brz, qitype=qitype, qival=qival, lat_range=lat_range)
@@ -803,18 +805,36 @@ def get_amv_winds_match(dist_threshold=150, qitype=amv_cqi_idx, qival=50, lat_ra
     noa_idxs = np.array(noa_idxs)
     nwc_idxs = np.array(nwc_idxs)
 
-
-    amvs_all = []
-
-    amvs_all.append(amvs_eum[eum_idxs, :])
-    amvs_all.append(amvs_brz[brz_idxs, :])
-    amvs_all.append(amvs_kma[kma_idxs, :])
-    amvs_all.append(amvs_jma[jma_idxs, :])
-    amvs_all.append(amvs_noa[noa_idxs, :])
-    amvs_all.append(amvs_nwc[nwc_idxs, :])
-
-    return amvs_all
-
+    # amvs_all = []
+    # amvs_all.append(amvs_eum[eum_idxs, :])
+    # amvs_all.append(amvs_brz[brz_idxs, :])
+    # amvs_all.append(amvs_kma[kma_idxs, :])
+    # amvs_all.append(amvs_jma[jma_idxs, :])
+    # amvs_all.append(amvs_noa[noa_idxs, :])
+    # amvs_all.append(amvs_nwc[nwc_idxs, :])
+
+    amvs_dict = {}
+    amvs_dict['EUM'] = amvs_eum[eum_idxs, :]
+    amvs_dict['BRZ'] = amvs_brz[brz_idxs, :]
+    amvs_dict['KMA'] = amvs_kma[kma_idxs, :]
+    amvs_dict['JMA'] = amvs_jma[jma_idxs, :]
+    amvs_dict['NOA'] = amvs_noa[noa_idxs, :]
+    amvs_dict['NWC'] = amvs_nwc[nwc_idxs, :]
+
+    return amvs_dict
+
+def run_best_fit_all(amvs_dict, raobs_pathname, dist_threshold=200, min_num_levs=20):
+    raob_dct = get_raob_dict(raobs_pathname)
+    print('got raobs')
+
+    bfs_dct = {}
+    for key in amv_centers_list:
+        amvs = amvs_dict.get(key)
+        bfs = run_best_fit(amvs, raob_dct, dist_threshold=dist_threshold, min_num_levs=min_num_levs)
+        print('done: ', key)
+        bfs_dct[key] = bfs
+
+    return bfs_dct
 
 def best_fit(amv_spd, amv_dir, amv_prs, amv_lat, amv_lon, fcst_spd, fcst_dir, fcst_prs, verbose=False):
 
@@ -1143,9 +1163,7 @@ def get_press_bin_ranges(lop, hip, bin_size=100):
     return bin_ranges
 
 
-def analyze_best_fit_all(amvs, bfs, bin_size=200):
-    num_centers = len(amvs)
-
+def analyze_best_fit_all(amvs_dct, bfs_dct, bin_size=200):
     pres_mad = []
     pres_bias = []
     spd_mad = []
@@ -1153,8 +1171,10 @@ def analyze_best_fit_all(amvs, bfs, bin_size=200):
     dir_mad = []
     dir_bias = []
 
-    for k in range(num_centers):
-        x_values, bin_pres, num_pres, bin_spd, num_spd, bin_dir, num_dir = analyze_best_fit_single(amvs[k], bfs[k], bin_size=bin_size)
+    for key in amv_centers_list:
+        amvs = amvs_dct.get(key)
+        bfs = bfs_dct.get(key)
+        x_values, bin_pres, num_pres, bin_spd, num_spd, bin_dir, num_dir = analyze_best_fit_single(amvs, bfs, bin_size=bin_size)
 
         pres_mad_c = []
         pres_bias_c = []
@@ -1188,6 +1208,7 @@ def analyze_best_fit_single(amvs, bfs, bin_size=200):
     # keep_idxs = bfs[keep, 0]
     # keep_idxs = keep_idxs.astype(np.int32)
     keep_idxs = bfs[:, 0]
+    keep_idxs = keep_idxs.astype(np.int32)
 
     amv_p = amvs[keep_idxs, amv_pres_idx]
     bf_p = bfs[:, 3]