first crack best fit altitude coords

modules/amv/intercompare.py

@@ -1039,6 +1039,79 @@ 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):
+    verbose = False
+
+    fcst_num_levels = fcst_spd.shape[0]
+
+    alt_diff = 1000.0
+    alt_top = 20000.0
+    alt_bot = 0
+
+    flag = 3
+    bf_tup = (None, None, None, 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))
+        return bf_tup
+
+    # Calculate the height +/- alt_diff from the AMV height
+    alt_max = min((amv_alt + alt_diff), alt_top)
+    alt_min = max((amv_lat - alt_diff), alt_bot)
+
+    # 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 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
+
+    # Diagnostic field: Find the model minimum speed and maximum speed within PressDiff of the AMV.
+    if verbose:
+        sat_wind_min_spd = min(fcst_spd[kk])
+        sat_wind_max_spd = 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((np.pi/180.0)*fcst_dir)
+    fcst_vwind = -fcst_spd * np.cos((np.pi/180.0)*fcst_dir)
+
+    # Calculate the vector difference between the AMV and model background at all levels.
+    vec_diff = np.sqrt((amv_uwind - fcst_uwind) ** 2 + (amv_vwind - fcst_vwind) ** 2)
+
+    # 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 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
+
+    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]
+
+    a1 = fcst_alt[imin-1]
+    v2 = vec_diff[imin-1]
+    a3 = fcst_alt[imin+1]
+    v3 = vec_diff[imin+1]
+
+    return bf_tup
+
 def get_press_bin_ranges(lop, hip, bin_size=100):
     bin_ranges = []
     delp = hip - lop