From 72d9ce45a48cb3a3c6b2927651a9cd2fcdfc28dc Mon Sep 17 00:00:00 2001
From: tomrink <rink@ssec.wisc.edu>
Date: Wed, 6 Jul 2022 13:28:17 -0500
Subject: [PATCH] improvements to best_fit...
---
modules/amv/intercompare.py | 66 ++++++++++++++++++++++---------------
1 file changed, 39 insertions(+), 27 deletions(-)
diff --git a/modules/amv/intercompare.py b/modules/amv/intercompare.py
index e3ae4142..8e96671c 100644
--- a/modules/amv/intercompare.py
+++ b/modules/amv/intercompare.py
@@ -112,6 +112,8 @@ def filter_amvs(amvs, qitype=amv_cqi_idx, qival=50, lat_range=[-61, 61]):
def get_raob_dict_cdf(pathname):
rtgrp = Dataset(pathname, 'r', format='NETCDF3')
+ wmo_stat_id_var = rtgrp['wmoStat']
+
sta_lats_var = rtgrp['staLat']
sta_lats_var.set_always_mask(False)
sta_lats_var.set_auto_mask(False)
@@ -122,12 +124,14 @@ def get_raob_dict_cdf(pathname):
sta_lats = sta_lats_var[:]
sta_lons = sta_lons_var[:]
+ sta_ids = wmo_stat_id_var[:]
sta_msk = np.logical_and(sta_lats > -90.0, sta_lats < 90.0)
# sta_msk = np.invert(sta_lats._mask)
num_sta = np.sum(sta_msk)
sta_lats = sta_lats[sta_msk]
sta_lons = sta_lons[sta_msk]
+ sta_ids = sta_ids[sta_msk]
# mandatory levels
man_levs = rtgrp['prMan']
@@ -171,6 +175,7 @@ def get_raob_dict_cdf(pathname):
for k in range(num_sta):
lon = sta_lons[k]
lat = sta_lats[k]
+ id = sta_ids[k]
vld_man = np.invert(man_levs[k]._mask)
vld_man = np.logical_and(vld_man, np.invert(man_spd[k]._mask))
@@ -976,23 +981,24 @@ def best_fit(amv_spd, amv_dir, amv_prs, amv_lat, amv_lon, fcst_spd, fcst_dir, fc
fcst_num_levels = fcst_prs.shape[0]
- PressDiff = 150.0 # pressure above and below AMV to look for fit
- TopPress = 50.0 # highest level to allow search
+ PressDiff = 150.0 # pressure above and below AMV to look for fit
+ TopPress = 50.0 # highest level to allow search
+ BotPress = fcst_prs[0] # lowest level to allow search
if amv_prs < TopPress:
if verbose:
- print('AMV location lat,lon,prs ({0},{1},{2}) is higher than pressure {3}'.format(amv_lat,amv_lon,amv_prs,TopPress))
+ print('AMV location lat,lon,prs ({0},{1},{2}) is higher than pressure {3}'.format(amv_lat, amv_lon, amv_prs, TopPress))
return bf_data
# Calculate the pressure +/- 150 hPa from the AMV pressure.
- PressMax = amv_prs + PressDiff
+ PressMax = min((amv_prs + PressDiff), BotPress)
PressMin = max((amv_prs - PressDiff), TopPress)
# 1d array of indicies to consider for best fit location
kk = np.where((fcst_prs < PressMax) & (fcst_prs > PressMin))
if len(kk[0]) == 0:
if verbose:
- print('AMV location lat,lon,prs ({0},{1},{2}) failed to find fcst prs around AMV'.format(amv_lat,amv_lon,amv_prs))
+ print('AMV location lat,lon,prs ({0},{1},{2}) failed to find fcst prs around AMV'.format(amv_lat, amv_lon, amv_prs))
return bf_data
# Diagnostic field: Find the model minimum speed and maximum speed within PressDiff of the AMV.
@@ -1015,7 +1021,7 @@ def best_fit(amv_spd, amv_dir, amv_prs, amv_lat, amv_lon, fcst_spd, fcst_dir, fc
imin = np.where(VecDiff == MinVecDiff)
if imin[0].size == 0:
if verbose:
- print('AMV location lat,lon,prs ({0},{1},{2}) failed to find the min vector difference in layers around AMV'.format(amv_lat,amv_lon,amv_prs))
+ print('AMV location lat,lon,prs ({0},{1},{2}) failed to find the min vector difference in layers around AMV'.format(amv_lat, amv_lon, amv_prs))
return bf_data
imin = imin[0][0] # index of min vec diff in absolute coordinates
@@ -1038,44 +1044,50 @@ def best_fit(amv_spd, amv_dir, amv_prs, amv_lat, amv_lon, fcst_spd, fcst_dir, fc
# if bottom model level
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 or bottom'.format(amv_lat,amv_lon,amv_prs))
+ print('AMV location lat,lon,alt ({0},{1},{2}) min vector difference at top or bottom'.format(amv_lat, amv_lon, amv_prs))
else:
p3 = fcst_prs[imin+1]
p1 = fcst_prs[imin-1]
v3 = VecDiff[imin+1]
v1 = VecDiff[imin-1]
- if p3 > TopPress and v1 != v2 and v2 != v3: # check not collinear, above allowed region
+ if p3 >= TopPress and v1 != v2 and v2 != v3: # check not collinear, above allowed region
SatwindBestFitPress = p2 - (0.5 *
((((p2 - p1) * (p2 - p1) * (v2 - v3)) - ((p2 - p3) * (p2 - p3) * (v2 - v1))) /
(((p2 - p1) * (v2 - v3)) - ((p2 - p3) * (v2 - v1)))))
if (SatwindBestFitPress < p3) or (SatwindBestFitPress > p1):
if (verbose):
print('Best Fit not found between two pressure layers')
- print('SatwindBestFitPress {0} p1 {1} p2 {2} p3 {3} imin {4}'.format(SatwindBestFitPress,p1,p2,p3,imin))
+ print('SatwindBestFitPress {0} p1 {1} p2 {2} p3 {3} imin {4}'.format(SatwindBestFitPress, p1, p2, p3, imin))
SatwindBestFitPress = p2
-
-# Find best fit U and V by linear interpolation.
- if p2 != SatwindBestFitPress:
- if p2 < SatwindBestFitPress:
- LevBelow = imin - 1
- LevAbove = imin
- Prop = (SatwindBestFitPress - p1) / (p2 - p1)
else:
- LevBelow = imin
- LevAbove = imin + 1
- Prop = (SatwindBestFitPress - p2) / (p3 - p2)
+ if verbose:
+ print('Top of parabolic fit window above TopPress, or co-linear vector difference profile')
+ print('SatwindBestFitPress {0} p1 {1} p2 {2} p3 {3} imin {4}'.format(SatwindBestFitPress, p1, p2, p3, imin))
+ print('Vector difference profile p1 {0} p2 {1} p3 {2}'.format(v1, v2, v3))
+
+ # Find best fit U and V by linear interpolation.
+ if p2 != SatwindBestFitPress:
+ if p2 < SatwindBestFitPress:
+ LevBelow = imin - 1
+ LevAbove = imin
+ Prop = (SatwindBestFitPress - p1) / (p2 - p1)
+ else:
+ LevBelow = imin
+ LevAbove = imin + 1
+ Prop = (SatwindBestFitPress - p2) / (p3 - p2)
- SatwindBestFitU = fcst_uwind[LevBelow] * (1.0 - Prop) + fcst_uwind[LevAbove] * Prop
- SatwindBestFitV = fcst_vwind[LevBelow] * (1.0 - Prop) + fcst_vwind[LevAbove] * Prop
+ SatwindBestFitU = fcst_uwind[LevBelow] * (1.0 - Prop) + fcst_uwind[LevAbove] * Prop
+ SatwindBestFitV = fcst_vwind[LevBelow] * (1.0 - Prop) + fcst_vwind[LevAbove] * Prop
+ # best fit success and good constraints -----
SatwindGoodConstraint = 1
flag = 0
# Check for overall good agreement
- vdiff = np.sqrt((fcst_uwind - SatwindBestFitU) ** 2 + (fcst_vwind - SatwindBestFitV) ** 2)
- min_vdiff = min(vdiff)
- if min_vdiff > 4.0:
+ vdiff_best = np.sqrt((fcst_uwind - SatwindBestFitU) ** 2 + (fcst_vwind - SatwindBestFitV) ** 2)
+ min_vdiff_best = min(vdiff_best)
+ if min_vdiff_best > 4.0:
SatwindGoodConstraint = 0
flag = 2
@@ -1083,7 +1095,7 @@ 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(vdiff[mm] < (min_vdiff + 2.0)) + np.sum(vdiff[nn] < (min_vdiff + 2.0))) > 0:
+ if (np.sum(vdiff_best[mm] < (min_vdiff_best + 2.0)) + np.sum(vdiff_best[nn] < (min_vdiff_best + 2.0))) > 0:
SatwindGoodConstraint = 0
flag = 1
@@ -1103,9 +1115,9 @@ def best_fit(amv_spd, amv_dir, amv_prs, amv_lat, amv_lon, fcst_spd, fcst_dir, fc
print('Bestfit -> pbest,bfu,bfv,amvu,amvv,bgu,bgv: {0} {1} {2} {3} {4} {5} {6}'.format(
SatwindBestFitPress,SatwindBestFitU,SatwindBestFitV,amv_uwind,amv_vwind,fcst_uwind[imin],fcst_vwind[imin]))
print('Good Constraint: {0}'.format(SatwindGoodConstraint))
- print('Minimum Vector Difference: {0}'.format(VecDiff[imin]))
+ print('Minimum Vector Difference: {0}'.format(min_vdiff_best))
print('Vector Difference Profile: ')
- print(vdiff)
+ print(vdiff_best)
print('Pressure Profile: ')
print(fcst_prs)
--
GitLab