diff --git a/modules/amv/intercompare.py b/modules/amv/intercompare.py
index ea3c9635dba4cc703fc49eb7b75173640abd14d4..96ed24e5ec5335145a2ef678288c5a4464c88a94 100644
--- a/modules/amv/intercompare.py
+++ b/modules/amv/intercompare.py
@@ -1032,7 +1032,7 @@ def best_fit(amv_spd, amv_dir, amv_prs, amv_lat, amv_lon, fcst_spd, fcst_dir, fc
print('Pressure Profile: ')
print(fcst_prs)
- if (abs(SatwindBestFitU - amv_uwind) > 4.0) or (abs(SatwindBestFitV - amv_vwind) > 4.0):
+ if (abs(SatwindBestFitU - amv_uwind) > 4.0) or (abs(SatwindBestFitV - amv_vwind) > 4.0):
print('U Diff: {0}'.format(abs(SatwindBestFitU - amv_uwind)))
print('V Diff: {0}'.format(abs(SatwindBestFitV - amv_vwind)))
@@ -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,49 +1092,77 @@ 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]
-
- a1 = fcst_alt[imin-1]
- v2 = vec_diff[imin-1]
- a3 = fcst_alt[imin+1]
- v3 = vec_diff[imin+1]
-
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]
+ 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:
+ if a2 < sat_wind_best_fit_alt:
+ lev_below = imin
+ lev_above = imin + 1
+ prop = (sat_wind_best_fit_alt - a2) / (a3 - a2)
+ else:
+ lev_below = imin - 1
+ lev_above = imin
+ prop = (sat_wind_best_fit_alt - a1) / (a2 - a1)
+
+ 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
-
-
-
-
-
- # 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]
+ if good_constraint == 1:
+ bf_tup = (sat_wind_best_fit_u, sat_wind_best_fit_v, sat_wind_best_fit_alt, 0)
else:
- if a2 < sat_wind_best_fit_alt:
- lev_below = imin
- lev_above = imin + 1
- prop = (sat_wind_best_fit_alt - a2) / (a3 - a2)
- else:
- lev_below = imin - 1
- lev_above = imin
- prop = (sat_wind_best_fit_alt - a1) / (a2 - a1)
- pass
+ 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)
- 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
+ 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