snapshot...

modules/amv/intercompare.py

@@ -1041,17 +1041,16 @@ def best_fit(amv_spd, amv_dir, amv_prs, amv_lat, amv_lon, fcst_spd, fcst_dir, fc
     return bf_data
 
 
-def best_fit_altitude(amv_spd, amv_dir, amv_alt, amv_lat, amv_lon, fcst_spd, fcst_dir, fcst_alt, test_a=True, test_b=True):
+def best_fit_altitude(amv_spd, amv_dir, amv_alt, amv_lat, amv_lon, fcst_spd, fcst_dir, fcst_alt,
+                      alt_top=25000.0, alt_bot=0.0):
     verbose = False
-
     fcst_num_levels = fcst_spd.shape[0]
 
     alt_diff = 1000.0
-    alt_top = 20000.0
-    alt_bot = 0
+    alt_diff_test = 500.0
 
     flag = 3
-    bf_tup = (None, None, None, flag)
+    bf_tup = (np.nan, np.nan, np.nan, flag)
     if amv_alt > alt_top:
         if verbose:
             print('AMV location lat,lon,lat ({0},{1},{2}) is higher than max altitude {3}'.format(amv_lat, amv_lon, amv_alt, alt_top))
@@ -1063,7 +1062,7 @@ def best_fit_altitude(amv_spd, amv_dir, amv_alt, amv_lat, amv_lon, fcst_spd, fcs
 
     # 1d array of indices to consider for best fit height
     kk = np.where((fcst_alt > alt_min) & (fcst_alt < alt_max))
-    if len(kk[0]) == 0:
+    if kk[0].size == 0:
         if verbose:
             print('AMV location lat,lon,alt ({0},{1},{2}) failed to find fcst alt around AMV'.format(amv_lat,amv_lon,amv_alt))
         return bf_tup
@@ -1085,7 +1084,7 @@ def best_fit_altitude(amv_spd, amv_dir, amv_alt, amv_lat, amv_lon, fcst_spd, fcs
     # Find the model level of best-fit pressure, from the minimum vector difference.
     min_vec_diff = min(vec_diff[kk])
     imin = np.where(vec_diff == min_vec_diff)
-    if len(imin[0]) == 0:
+    if imin[0].size == 0:
         if verbose:
             print('AMV location lat,lon,alt ({0},{1},{2}) failed to find min vector difference in layers around AMV'.format(amv_lat,amv_lon,amv_alt))
         return bf_tup
@@ -1093,37 +1092,29 @@ def best_fit_altitude(amv_spd, amv_dir, amv_alt, amv_lat, amv_lon, fcst_spd, fcs
     imin = imin[0][0]
     imin = kk[0][imin]
 
-    if imin == fcst_num_levels - 1:
-        if verbose:
-            print('AMV location lat,lon,alt ({0},{1},{2}) failed to find min vector difference in layers around AMV'.format(amv_lat,amv_lon,amv_alt))
-        return bf_tup
-
-    a2 = fcst_alt[imin]
-    v2 = vec_diff[imin]
-
-    if imin == 0:
-        return bf_tup
-
     a2 = fcst_alt[imin]
     v2 = vec_diff[imin]
 
+    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))
+    else:
         a1 = fcst_alt[imin-1]
-    v2 = vec_diff[imin-1]
+        v1 = vec_diff[imin-1]
         a3 = fcst_alt[imin+1]
         v3 = vec_diff[imin+1]
 
+        # check not co-linear
+        if a1 != a2 and a2 != a3:
+            sat_wind_best_fit_alt = minimize_quadratic(a1, a2, a3, v1, v2, v3)
+            if sat_wind_best_fit_alt < a1 or sat_wind_best_fit_alt > a3:
                 sat_wind_best_fit_alt = a2
 
-
-
-
-
-
         # Find best fit U and V by linear interpolation:
-    if a2 == sat_wind_best_fit_alt:
-        sat_wind_best_fit_u = fcst_uwind[imin]
-        sat_wind_best_fit_v = fcst_vwind[imin]
-    else:
+        if a2 != sat_wind_best_fit_alt:
             if a2 < sat_wind_best_fit_alt:
                 lev_below = imin
                 lev_above = imin + 1
@@ -1132,11 +1123,47 @@ def best_fit_altitude(amv_spd, amv_dir, amv_alt, amv_lat, amv_lon, fcst_spd, fcs
                 lev_below = imin - 1
                 lev_above = imin
                 prop = (sat_wind_best_fit_alt - a1) / (a2 - a1)
-            pass
 
             sat_wind_best_fit_u = fcst_uwind[lev_below] * (1.0 - prop) + fcst_uwind[lev_above] * prop
             sat_wind_best_fit_v = fcst_vwind[lev_below] * (1.0 - prop) + fcst_vwind[lev_above] * prop
 
+    # check contraints
+    good_constraint = 1
+
+    vdiff = np.sqrt((amv_uwind - sat_wind_best_fit_u) ** 2 + (amv_vwind - sat_wind_best_fit_v) ** 2)
+    if vdiff > 4.0:
+        good_constraint = 0
+        flag = 2
+
+    mm = np.where(fcst_alt > (sat_wind_best_fit_alt + alt_diff_test))[0]
+    nn = np.where(fcst_alt < (sat_wind_best_fit_alt - alt_diff_test))[0]
+
+    if (np.sum(vec_diff[mm] < (vec_diff + 2.0)) + np.sum(vec_diff[nn] < (vec_diff + 2.0))) > 0:
+        good_constraint = 0
+        flag = 1
+
+    if good_constraint == 1:
+        bf_tup = (sat_wind_best_fit_u, sat_wind_best_fit_v, sat_wind_best_fit_alt, 0)
+    else:
+        bf_tup = (None, None, None, flag)
+
+    if verbose:
+        print('*** AMV best-fit ***')
+        print('AMV -> p/minspd/maxspd: {0} {1} {2}'.format(amv_alt,sat_wind_min_spd,sat_wind_max_spd))
+        # 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(
+        sat_wind_best_fit_alt,sat_wind_best_fit_u,sat_wind_best_fit_v,amv_uwind,amv_vwind,fcst_uwind[imin],fcst_vwind[imin]))
+        print('Good Constraint: {0}'.format(good_constraint))
+        print('Minimum Vector Difference: {0}'.format(vec_diff[imin]))
+        print('Vector Difference Profile: ')
+        print(vec_diff)
+        print('Altitude Profile: ')
+        print(fcst_alt)
+
+        if (abs(sat_wind_best_fit_u - amv_uwind) > 4.0) or (abs(sat_wind_best_fit_v - amv_vwind) > 4.0):
+            print('U Diff: {0}'.format(abs(sat_wind_best_fit_u - amv_uwind)))
+            print('V Diff: {0}'.format(abs(sat_wind_best_fit_v - amv_vwind)))
+
     return bf_tup