snapshot...

modules/amv/intercompare.py

@@ -419,25 +419,22 @@ def run_best_fit(amvs, raob_dct, dist_threshold=200, min_num_levs=20):
     for j in range(num_amvs):
         tup = raob_nbor_dct[j]
         if tup is not None:
+            cnt += 1
             nbor_keys = tup[0]
             nbor_dsts = tup[1]
+            s_i = np.argsort(nbor_dsts).tolist()
 
-            last_dist = dist_threshold
-            keep_raob = None
-            for k, ll_tup in enumerate(nbor_keys):
-                dist = nbor_dsts[k]
+            for k in range(len(nbor_keys)):
+                dist = nbor_dsts[s_i[k]]
+                ll_tup = nbor_keys[s_i[k]]
                 raob = raob_dct.get(ll_tup)
 
                 if raob.shape[0] < min_num_levs:
                     continue
-                if dist < last_dist:
-                    last_dist = dist
-                    keep_raob = raob
 
-            if keep_raob is not None:
-                rspd = keep_raob[:, raob_spd_idx]
-                rdir = keep_raob[:, raob_dir_idx]
-                rprs = keep_raob[:, raob_pres_idx]
+                rspd = raob[:, raob_spd_idx]
+                rdir = raob[:, raob_dir_idx]
+                rprs = raob[:, raob_pres_idx]
 
                 aspd = amvs[j, amv_spd_idx]
                 adir = amvs[j, amv_dir_idx]
@@ -445,9 +442,11 @@ def run_best_fit(amvs, raob_dct, dist_threshold=200, min_num_levs=20):
                 alat = amvs[j, amv_lat_idx]
                 alon = amvs[j, amv_lon_idx]
 
-                bf = best_fit(aspd, adir, aprs, alat, alon, rspd, rdir, rprs)
+                bf = best_fit(aspd, adir, aprs, alat, alon, rspd, rdir, rprs, verbose=True)
                 bf = (j, bf[0], bf[1], bf[2], bf[3])
-                bestfits.append(bf)
+                if bf[4] == 0:
+                    bestfits.append(bf)
+                    break
 
     return bestfits
 
@@ -954,8 +953,10 @@ def best_fit(amv_spd, amv_dir, amv_prs, amv_lat, amv_lon, fcst_spd, fcst_dir, fc
     SatwindGoodConstraint = 1
     flag = 0
     # Check for overall good agreement
-    vdiff = np.sqrt((amv_uwind - SatwindBestFitU) ** 2 + (amv_vwind - SatwindBestFitV) ** 2)
-    if vdiff > 4.0:
+    #vdiff = np.sqrt((amv_uwind - SatwindBestFitU) ** 2 + (amv_vwind - SatwindBestFitV) ** 2)
+    vdiff = np.sqrt((fcst_uwind - SatwindBestFitU) ** 2 + (fcst_vwind - SatwindBestFitV) ** 2)
+    min_vdiff = min(vdiff)
+    if min_vdiff > 4.0:
         SatwindGoodConstraint = 0
         flag = 2
 
@@ -963,7 +964,8 @@ def best_fit(amv_spd, amv_dir, amv_prs, amv_lat, amv_lon, fcst_spd, fcst_dir, fc
     if SatwindGoodConstraint == 1:
         mm = np.where(fcst_prs > (SatwindBestFitPress + 100))[0]
         nn = np.where(fcst_prs < (SatwindBestFitPress - 100))[0]
-        if (np.sum(VecDiff[mm] < (vdiff + 2.0)) + np.sum(VecDiff[nn] < (vdiff + 2.0))) > 0:
+        #if (np.sum(VecDiff[mm] < (vdiff + 2.0)) + np.sum(VecDiff[nn] < (vdiff + 2.0))) > 0:
+        if (np.sum(vdiff[mm] < (min_vdiff + 2.0)) + np.sum(vdiff[nn] < (min_vdiff + 2.0))) > 0:
             SatwindGoodConstraint = 0
             flag = 1
 
@@ -977,7 +979,7 @@ def best_fit(amv_spd, amv_dir, amv_prs, amv_lat, amv_lon, fcst_spd, fcst_dir, fc
         bfit_v = undef
 
     if verbose:
-        print('*** AMV best-fit ***')
+        print('*** AMV best-fit *********************************')
         print('AMV -> p/minspd/maxspd: {0} {1} {2}'.format(amv_prs,SatwindMinSpeed,SatwindMaxSpeed))
         # print('Bestfit -> p1,p2,p3,v1,v2,v3: {0} {1} {2} {3} {4} {5}'.format(p1,p2,p3,v1,v2,v3))
         print('Bestfit -> pbest,bfu,bfv,amvu,amvv,bgu,bgv: {0} {1} {2} {3} {4} {5} {6}'.format(
@@ -985,7 +987,7 @@ def best_fit(amv_spd, amv_dir, amv_prs, amv_lat, amv_lon, fcst_spd, fcst_dir, fc
         print('Good Constraint: {0}'.format(SatwindGoodConstraint))
         print('Minimum Vector Difference: {0}'.format(VecDiff[imin]))
         print('Vector Difference Profile: ')
-        print(VecDiff)
+        print(vdiff)
         print('Pressure Profile: ')
         print(fcst_prs)
 
@@ -1090,8 +1092,9 @@ def best_fit_altitude(amv_spd, amv_dir, amv_alt, amv_lat, amv_lon, fcst_spd, fcs
     good_constraint = 1
     flag = 0
     # Check for overall good agreement: want vdiff less than 4 m/s
-    vdiff = np.sqrt((amv_uwind - sat_wind_best_fit_u) ** 2 + (amv_vwind - sat_wind_best_fit_v) ** 2)
-    if vdiff > 4.0:
+    vdiff = np.sqrt((fcst_uwind - sat_wind_best_fit_u) ** 2 + (fcst_vwind - sat_wind_best_fit_v) ** 2)
+    min_vdiff = min(vdiff)
+    if min_vdiff > 4.0:
         good_constraint = 0
         flag = 2
 
@@ -1099,7 +1102,7 @@ def best_fit_altitude(amv_spd, amv_dir, amv_alt, amv_lat, amv_lon, fcst_spd, fcs
     if good_constraint == 1:
         mm = np.where(fcst_alt > (sat_wind_best_fit_alt + constraint_half_width))[0]
         nn = np.where(fcst_alt < (sat_wind_best_fit_alt - constraint_half_width))[0]
-        if (np.sum(vec_diff[mm] < (vdiff + 2.0)) + np.sum(vec_diff[nn] < (vdiff + 2.0))) > 0:
+        if (np.sum(vdiff[mm] < (min_vdiff + 2.0)) + np.sum(vdiff[nn] < (min_vdiff + 2.0))) > 0:
             good_constraint = 0
             flag = 1
 
@@ -1117,7 +1120,7 @@ def best_fit_altitude(amv_spd, amv_dir, amv_alt, amv_lat, amv_lon, fcst_spd, fcs
         print('Good Constraint: {0}'.format(good_constraint))
         print('Minimum Vector Difference: {0}'.format(vec_diff[imin]))
         print('Vector Difference Profile: ')
-        print(vec_diff)
+        print(vdiff)
         print('Altitude Profile: ')
         print(fcst_alt)