snapshot..

modules/amv/intercompare.py

@@ -834,7 +834,7 @@ def get_amv_winds_match(dist_threshold=150, cqi=50):
     return amvs_all
 
 
-def best_fit(amv_spd, amv_dir, amv_prs, amv_lat, amv_lon, fcst_spd, fcst_dir, fcst_prs, test_a=True, test_b=True):
+def best_fit(amv_spd, amv_dir, amv_prs, amv_lat, amv_lon, fcst_spd, fcst_dir, fcst_prs):
 
     """Finds the background model best fit pressure associated with the AMV.
        The model best-fit pressure is the height (in pressure units) where the
@@ -903,11 +903,9 @@ def best_fit(amv_spd, amv_dir, amv_prs, amv_lat, amv_lon, fcst_spd, fcst_dir, fc
         SatwindMaxSpeed = max(fcst_spd[kk])
 
 # Compute U anv V for both AMVs and forecast
-    amv_uwind = -amv_spd * np.sin((np.pi/180.0)*amv_dir)
-    amv_vwind = -amv_spd * np.cos((np.pi/180.0)*amv_dir)
-#   fcst_uwind = -fcst_spd[:] * np.sin(math.radians(fcst_dir[:]))
-#   fcst_vwind = -fcst_spd[:] * np.cos(math.radians(fcst_dir[:]))
-    dr=0.017453
+    dr = 0.0174533  # pi/180
+    amv_uwind = -amv_spd * np.sin(dr*amv_dir)
+    amv_vwind = -amv_spd * np.cos(dr*amv_dir)
     fcst_uwind = -fcst_spd * np.sin(dr*fcst_dir)
     fcst_vwind = -fcst_spd * np.cos(dr*fcst_dir)
 
@@ -916,22 +914,19 @@ def best_fit(amv_spd, amv_dir, amv_prs, amv_lat, amv_lon, fcst_spd, fcst_dir, fc
 
 # Find the model level of best-fit pressure, from the minimum vector difference.
     MinVecDiff = min(VecDiff[kk])
-    imin = -1
-    for i, item in enumerate(VecDiff):
-        if MinVecDiff == VecDiff[i]:
-            if i in kk[0]:
-                imin = i
-
-    if imin == -1:
+    imin = np.where(VecDiff == MinVecDiff)
+    if imin[0].size == 0:
         if verbose:
             print('AMV location lat,lon,prs ({0},{1},{2}) failed to find min vector difference in layers around AMV'.format(amv_lat,amv_lon,amv_prs))
         return bf_data
 
+    imin = imin[0][0]
+    imin = kk[0][imin]
+
     # TDR - Feb 2020
     if imin == fcst_num_levels - 1:
         if verbose:
             print('AMV location lat,lon,prs ({0},{1},{2}) failed to find min vector difference in layers around AMV'.format(amv_lat,amv_lon,amv_prs))
-
         return bf_data
 
 
@@ -948,8 +943,10 @@ def best_fit(amv_spd, amv_dir, amv_prs, amv_lat, amv_lon, fcst_spd, fcst_dir, fc
 
 # assumes fcst data level 0 at surface and (fcst_num_levels-1) at model top
 # if bottom model level
-    if imin == 0:
+    if imin == 0 or imin == fcst_num_levels - 1:
         SatwindBestFitPress = p2
+        if verbose:
+            print('AMV location lat,lon,alt ({0},{1},{2}) min vector difference at top/bottom'.format(amv_lat,amv_lon,amv_prs))
     else:
         p3 = fcst_prs[imin+1]
         p1 = fcst_prs[imin-1]
@@ -992,22 +989,20 @@ def best_fit(amv_spd, amv_dir, amv_prs, amv_lat, amv_lon, fcst_spd, fcst_dir, fc
 
 
 # Check to see if the best fit pressure is constrained.
-    SatwindGoodConstraint = 0
-    flag = 2
+    SatwindGoodConstraint = 1
+    vdiff = np.sqrt((amv_uwind - SatwindBestFitU) ** 2 + (amv_vwind - SatwindBestFitV) ** 2)
+    if vdiff > 4.0:
+        SatwindGoodConstraint = 0
+        flag = 2
 
-    if not test_a:
-        SatwindGoodConstraint = 1
-    elif MinVecDiff <= 4.0:
-        SatwindGoodConstraint = 1
-    else:
+    # Check for a substantial secondary local minimum or if local minimum is too broad
+    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:
+        SatwindGoodConstraint = 0
         flag = 1
 
-    if test_a and test_b and SatwindGoodConstraint == 1:
-        for ilev in range(fcst_num_levels):
-            if fcst_prs[ilev] >= TopPress:
-                if ((fcst_prs[ilev] < (SatwindBestFitPress - 100.)) or (fcst_prs[ilev] > (SatwindBestFitPress + 100.))) and (VecDiff[ilev] <= (MinVecDiff + 2.0)):
-                    SatwindGoodConstraint = 0
-                    flag = 2
 
     if SatwindGoodConstraint == 1:
         bfit_prs = SatwindBestFitPress
@@ -1101,6 +1096,7 @@ def best_fit_altitude(amv_spd, amv_dir, amv_alt, amv_lat, amv_lon, fcst_spd, fcs
     sat_wind_best_fit_alt = a2
     sat_wind_best_fit_u = fcst_uwind[imin]
     sat_wind_best_fit_v = fcst_vwind[imin]
+
     if imin == 0 or imin == (fcst_num_levels - 1):
         if verbose:
             print('AMV location lat,lon,alt ({0},{1},{2}) min vector difference at top/bottom'.format(amv_lat,amv_lon,amv_alt))