Skip to content

Commit c5eac2a

Browse files
evanmasonevanmason
committed
Bugfix (incomplete) for Ellipse calculation
1 parent a573f89 commit c5eac2a

File tree

4 files changed

+106
-63
lines changed

4 files changed

+106
-63
lines changed

‎eddy_tracker_configuration.yaml‎

Lines changed: 12 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -13,34 +13,37 @@ DOMAIN:
1313
#THE_DOMAIN = 'BlackSea'
1414
#THE_DOMAIN = 'MedSea' # not yet implemented
1515
LONMIN: -40.
16-
LONMAX: -5.
17-
LATMIN: 20.
18-
LATMAX: 30.
19-
DATE_STR: 20020101
20-
DATE_END: 20020331
16+
LONMAX: -30.
17+
LATMIN: 30.
18+
LATMAX: 40.
19+
DATE_STR: 19980101
20+
DATE_END: 19980331
2121

2222
# Specify the AVISO product
2323
AVISO:
2424
AVISO_DT14: Yes # Use new AVISO DT14 data (e.g., Capet et al, 2014)
25-
AVISO_FILES: 'dt_global_twosat_msla_h_????????_20140106.nc'
25+
AVISO_FILES: 'dt_global_allsat_msla_h_????????_20140106.nc'
2626
AVISO_DT14_SUBSAMP: No # subsample the dqily AVISO DT14
2727
DAYS_BTWN_RECORDS: 7 # sampling rate (days)
2828

2929
PATHS:
3030
# Path to AVISO data
31-
DATA_DIR: '/data/PVA/Externe/global/delayed-time/grids/msla/two-sat-merged/h/2002/'
31+
#DATA_DIR: '/data/PVA/Externe/global/delayed-time/grids/msla/two-sat-merged/h/2002/'
32+
DATA_DIR: '/Users/emason/data/AVISO/global/delayed-time/grids/msla/all-sat-merged/h/'
3233
#DATA_DIR: '/marula/emason/data/altimetry/global/delayed-time/grids/msla/two-sat-merged/h/'
3334
# Path and filename of Chelton et al (1998) Rossby radius data
3435
# Obtain file from:
3536
# http://www-po.coas.oregonstate.edu/research/po/research/rossby_radius/
36-
RW_PATH: '/homelocal/emason/rossrad.dat'
37+
#RW_PATH: '/homelocal/emason/rossrad.dat'
38+
RW_PATH: '/Users/emason/Downloads/rossrad.dat'
3739
#RW_PATH: '/home/emason/data_tmp/chelton_eddies/rossrad.dat'
3840

3941
# Path for saving of outputs
4042
# SAVE_DIR: '/homelocal/emason/outputs/before_swirl/'
4143
# SAVE_DIR: '/homelocal/emason/outputs/after_swirl/'
4244
# SAVE_DIR: '/homelocal/emason/outputs/after_swirl_slice/'
43-
SAVE_DIR: '/homelocal/emason/outputs/'
45+
#SAVE_DIR: '/homelocal/emason/outputs/'
46+
SAVE_DIR: '/Users/emason/toto/'
4447
#SAVE_DIR: '/marula/emason/aviso_eddy_tracking/fd0cdb4b2bd9/'
4548

4649
# Reference Julian day (Julian date at Jan 1, 1992)

‎make_eddy_track_AVISO.py‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -853,8 +853,9 @@ def set_interp_coeffs(self, sla, uspd):
853853

854854

855855
# Instantiate search ellipse object
856-
search_ellipse = eddy_tracker.SearchEllipse(THE_DOMAIN, DAYS_BTWN_RECORDS,
857-
RW_PATH, [LONMIN, LONMAX, LATMIN, LATMAX])
856+
search_ellipse = eddy_tracker.SearchEllipse(THE_DOMAIN, sla_grd,
857+
DAYS_BTWN_RECORDS, sla_grd.ZERO_CROSSING, RW_PATH,
858+
[LONMIN, LONMAX, LATMIN, LATMAX])
858859

859860

860861
if 'Gaussian' in SMOOTHING_TYPE:

‎make_eddy_tracker_list_obj.py‎

Lines changed: 88 additions & 52 deletions
Original file line numberDiff line numberDiff line change
@@ -1203,50 +1203,58 @@ def reshape_mask_eff(self, grd):
12031203

12041204
class RossbyWaveSpeed (object):
12051205

1206-
def __init__(self, domain, limits, rw_path=None):
1206+
def __init__(self, domain, grd, ZERO_CROSSING, limits, rw_path=None):
12071207
"""
12081208
Initialise the RossbyWaveSpeedsobject
12091209
"""
12101210
self.domain = domain
1211+
self.M = grd.M
1212+
self.ZERO_CROSSING = ZERO_CROSSING
12111213
self.rw_path = rw_path
12121214
if self.domain in ('Global', 'ROMS'):
12131215
assert self.rw_path is not None, \
12141216
'Must supply a path for the Rossby deformation radius data'
12151217
data = np.loadtxt(rw_path)
12161218
self._lon = data[:, 1] - 360.
12171219
self._lat = data[:, 0]
1218-
self._defrad = data[:, 3]
1220+
self._c = data[:, 2]
1221+
#self._defrad = data[:, 3]
12191222
self.limits = limits
1220-
self._make_subset()
1223+
self._make_subset()._make_kdtree()
12211224
self.vartype = 'variable'
12221225
else:
12231226
self.vartype = 'constant'
12241227
self.distance = np.empty(1)
12251228
self.start = True
12261229

12271230

1228-
def get_rwdistance(self, x, y, days_between_records):
1231+
def get_rwdistance(self, xpt, ypt, DAYS_BTWN_RECORDS):
12291232
"""
12301233
"""
12311234
if self.domain in ('Global', 'ROMS'):
1232-
#print 'x,y', x,y
1233-
self.distance[:] = self._get_rlongwave_spd(x, y)
1235+
#print 'xpt, ypt', xpt, ypt
1236+
self.distance[:] = self._get_c(xpt, ypt)
12341237
self.distance *= 86400.
12351238
#if self.domain in 'ROMS':
12361239
#self.distance *= 1.5
12371240
elif 'BlackSea' in self.domain:
12381241
self.distance[:] = 15000. # e.g., Blokhina & Afanasyev, 2003
12391242
elif 'MedSea' in self.domain:
12401243
self.distance[:] = 20000.
1241-
else: Exception # Unknown domain
1244+
else:
1245+
Exception # Unknown domain
12421246

12431247
if self.start:
12441248
print '--------- setting ellipse for %s domain' %self.domain
1245-
print '--------- using %s Rossby deformation radius of %s m' \
1246-
%(self.vartype, np.abs(self.distance[0]))
1249+
#print '--------- using %s Rossby deformation radius of %s m' \
1250+
#%(self.vartype, np.abs(self.distance[0]))
1251+
print '--------- using %s baroclinic gravity wave phase speed of %s m/s' \
1252+
%(self.vartype, self.distance[0] / 86400.)
1253+
print '--------- ', xpt, ypt
12471254
self.start = False
1248-
self.distance *= days_between_records
1249-
return np.abs(self.distance)
1255+
self.distance *= DAYS_BTWN_RECORDS
1256+
#print 'aaaa', self.distance
1257+
return self.distance
12501258

12511259

12521260
def _make_subset(self):
@@ -1260,77 +1268,91 @@ def _make_subset(self):
12601268
self._lat <= LATMAX + pad)
12611269
self._lon = self._lon[lloi]
12621270
self._lat = self._lat[lloi]
1263-
self._defrad = self._defrad[lloi]
1264-
self._make_kdtree()
1271+
self._c = self._c[lloi]
1272+
self.x, self.y = self.M(self._lon, self._lat)
1273+
return self
12651274

12661275

12671276
def _make_kdtree(self):
1268-
points = np.vstack([self._lon, self._lat]).T
1277+
#points = np.vstack([self._lon, self._lat]).T
1278+
points = np.vstack([self.x, self.y]).T
12691279
self._tree = spatial.cKDTree(points)
12701280

12711281

1272-
def _get_defrad(self, x, y):
1282+
#def _get_defrad(self, x, y):
1283+
#"""
1284+
#Get a point average of the deformation radius
1285+
#at x, y
1286+
#"""
1287+
#weights, i = self._tree.query(np.array([x, y]), k=4, p=2)
1288+
#weights /= weights.sum()
1289+
#self._weights = weights
1290+
#self.i = i
1291+
#return np.average(self._defrad[i], weights=weights)
1292+
1293+
def _get_c(self, xpt, ypt):
12731294
"""
1274-
Get a point average of the deformation radius
1275-
at x, y
1295+
Get a point average of c, the phase speed, at x, y
12761296
"""
1277-
weights, i = self._tree.query(np.array([x, y]), k=4, p=2)
1297+
weights, i = self._tree.query(np.array([xpt, ypt]), k=4, p=2)
12781298
weights /= weights.sum()
12791299
self._weights = weights
12801300
self.i = i
1281-
return np.average(self._defrad[i], weights=weights)
1282-
1301+
return np.average(self._c[i], weights=weights)
12831302

1284-
def _get_rlongwave_spd(self, x, y):
1285-
"""
1286-
Get the longwave phase speed, see Chelton etal (2008) pg 446:
1287-
c = -beta * defrad**2 (this only for extratropical waves...)
1288-
"""
1289-
r_spd_long = self._get_defrad(x, y)
1290-
r_spd_long *= 1000. # km to m
1291-
r_spd_long **= 2
1292-
lat = np.average(self._lat[self.i], weights=self._weights)
1293-
beta = np.cos(np.deg2rad(lat))
1294-
beta *= 1458e-7 # 1458e-7 ~ (2 * 7.29*10**-5)
1295-
beta /= 6371315.0
1296-
r_spd_long *= -beta
1297-
return r_spd_long
1303+
#def _get_rlongwave_spd(self, x, y):
1304+
#"""
1305+
#Get the longwave phase speed, see Chelton etal (2008) pg 446:
1306+
#c = -beta * defrad**2 (this only for extratropical waves...)
1307+
#"""
1308+
##r_spd_long = self._get_c(x, y)
1309+
##r_spd_long *= 1000. # km to m
1310+
##r_spd_long **= 2
1311+
##lat = np.average(self._lat[self.i], weights=self._weights)
1312+
##beta = np.cos(np.deg2rad(lat))
1313+
##beta *= 1458e-7 # 1458e-7 ~ (2 * 7.29*10**-5)
1314+
##beta /= 6371315.0
1315+
##r_spd_long *= -beta
1316+
#return self._get_c(x, y)
12981317

12991318

13001319

13011320

13021321
class SearchEllipse (object):
13031322

1304-
def __init__(self, domain, DAYS_BTWN_RECORDS, rw_path, limits):
1323+
def __init__(self, domain, grd, DAYS_BTWN_RECORDS, ZERO_CROSSING,
1324+
rw_path, limits):
13051325
"""
13061326
Class to construct a search ellipse/circle around a specified point.
13071327
13081328
"""
13091329
self.domain = domain
13101330
self.DAYS_BTWN_RECORDS = DAYS_BTWN_RECORDS
1331+
self.ZERO_CROSSING = ZERO_CROSSING
13111332
self.rw_path = rw_path
13121333
self.limits = limits
13131334
self.rw_d_fac = 1.75
13141335
self.e_w_major = self.DAYS_BTWN_RECORDS * 3e5 / 7.
1315-
self.n_s_minor = self.DAYS_BTWN_RECORDS * 15e4 / 7.
1316-
self.rwv = RossbyWaveSpeed(self.domain,
1317-
self.limits, rw_path=self.rw_path)
1336+
self.n_s_minor = self.DAYS_BTWN_RECORDS * 2 * 15e4 / 7.
1337+
self.semi_n_s_minor = 0.5 * self.n_s_minor
1338+
self.rwv = RossbyWaveSpeed(domain, grd, ZERO_CROSSING,
1339+
limits, rw_path=rw_path)
13181340
self.rw_d = np.empty(1)
13191341
self.rw_d_mod = np.empty(1)
13201342

13211343

13221344
def _set_east_ellipse(self):
13231345
"""
13241346
"""
1325-
self.east_ellipse = patch.Ellipse((self.x, self.y),
1347+
self.east_ellipse = patch.Ellipse((self.xpt, self.ypt),
13261348
self.e_w_major, self.n_s_minor)
13271349
return self
13281350

13291351

13301352
def _set_west_ellipse(self):
13311353
"""
13321354
"""
1333-
self.west_ellipse = patch.Ellipse((self.x, self.y),
1355+
self.west_ellipse = patch.Ellipse((self.xpt, self.ypt),
13341356
self.rw_d_mod, self.n_s_minor)
13351357
return self
13361358

@@ -1362,27 +1384,35 @@ def _set_black_sea_ellipse(self):
13621384
return self
13631385

13641386

1365-
def set_search_ellipse(self, x, y):
1387+
def set_search_ellipse(self, xpt, ypt):
13661388
"""
13671389
"""
1368-
self.x = x
1369-
self.y = y
1390+
self.xpt = xpt
1391+
self.ypt = ypt
13701392

13711393
if self.domain in ('Global', 'ROMS'):
1372-
self.rw_d[:] = self.rwv.get_rwdistance(x, y,
1373-
self.DAYS_BTWN_RECORDS)
1394+
self.rw_d[:] = self.rwv.get_rwdistance(xpt, ypt,
1395+
self.DAYS_BTWN_RECORDS)
1396+
#print self.rw_d
1397+
#exit()
13741398
self.rw_d_mod[:] = 1.75
13751399
self.rw_d_mod *= self.rw_d
1376-
self.rw_d_mod[:] = np.maximum(self.rw_d_mod, self.n_s_minor)
1400+
self.rw_d_mod[:] = np.array([self.rw_d_mod,
1401+
self.semi_n_s_minor]).max()
1402+
1403+
#print self.rw_d_mod, self.semi_n_s_minor
1404+
13771405
self.rw_d_mod *= 2.
1406+
#print self.rw_d_mod,
13781407
self._set_global_ellipse()
13791408
elif 'BlackSea' in self.domain:
13801409
self.rw_d_mod[:] = 1.75
1381-
self.rw_d[:] = self.rwv.get_rwdistance(x, y,
1382-
self.DAYS_BTWN_RECORDS)
1410+
self.rw_d[:] = self.rwv.get_rwdistance(xpt, ypt,
1411+
self.DAYS_BTWN_RECORDS)
13831412
self.rw_d_mod *= self.rw_d
13841413
self._set_black_sea_ellipse()
1385-
else: Exception
1414+
else:
1415+
Exception
13861416

13871417
return self
13881418

@@ -1392,12 +1422,18 @@ def view_search_ellipse(self, Eddy):
13921422
Input A_eddy or C_eddy
13931423
"""
13941424
plt.figure()
1395-
ax = plt.subplot(111)
1396-
ax.set_title('Rossby def. rad %s m' %self.rw_d[0])
1397-
Eddy.M.scatter(self.x, self.y, c='b')
1425+
ax = plt.subplot()
1426+
#ax.set_title('Rossby def. rad %s m' %self.rw_d[0])
1427+
Eddy.M.scatter(self.xpt, self.ypt, c='b')
13981428
Eddy.M.plot(self.ellipse_path.vertices[:, 0],
13991429
self.ellipse_path.vertices[:, 1], 'r')
14001430
Eddy.M.drawcoastlines()
1431+
stride = .5
1432+
Eddy.M.drawparallels(np.arange(-90, 90.+stride, stride),
1433+
labels=[1, 0, 0, 0], ax=ax)
1434+
Eddy.M.drawmeridians(np.arange(-360, 360. + stride, stride),
1435+
labels=[0, 0, 0, 1], ax=ax)
1436+
#plt.axis('image')
14011437
plt.show()
14021438

14031439

‎py_eddy_tracker_classes.py‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -994,6 +994,9 @@ def track_eddies(Eddy, first_record):
994994
if 'ellipse' in Eddy.SEPARATION_METHOD:
995995
if Eddy.search_ellipse.ellipse_path.contains_point(
996996
(new_x[new_ind], new_y[new_ind])):
997+
#if Eddy.search_ellipse.ellipse_path.contains_point(
998+
#(Eddy.new_lon_tmp[new_ind], Eddy.new_lat_tmp[new_ind])):
999+
#Eddy.search_ellipse.view_search_ellipse(Eddy)
9971000
sep_proceed = True
9981001
else:
9991002
sep_proceed = False

0 commit comments

Comments
 (0)