Skip to content

Commit e666035

Browse files
evanmasonevanmason
committed
Rewrite of get_circle function
1 parent 3643eae commit e666035

File tree

1 file changed

+55
-41
lines changed

1 file changed

+55
-41
lines changed

‎py_eddy_tracker_classes.py‎

Lines changed: 55 additions & 41 deletions
Original file line numberDiff line numberDiff line change
@@ -95,8 +95,8 @@ def do_basemap(M, ax):
9595
labels=[1, 0, 0, 0], linewidth=0.25, size=8, ax=ax)
9696
M.drawmeridians(np.arange(-360, 360 + stride, stride),
9797
labels=[0, 0, 0, 1], linewidth=0.25, size=8, ax=ax)
98-
M.fillcontinents('k', ax=ax)
99-
M.drawcoastlines(linewidth=0.5, ax=ax)
98+
M.fillcontinents('k', ax=ax, zorder=20)
99+
M.drawcoastlines(linewidth=0.5, ax=ax, zorder=20)
100100
return
101101

102102

@@ -474,7 +474,7 @@ def get_uavg(Eddy, CS, collind, centlon_e, centlat_e, poly_eff,
474474
seglon, seglat = (poly_i.vertices[:, 0],
475475
poly_i.vertices[:, 1])
476476
seglon, seglat = eddy_tracker.uniform_resample(
477-
seglon, seglat, method='akima')
477+
seglon, seglat, method='interp1d')
478478

479479
# Interpolate uspd to seglon, seglat, then get mean
480480
if 'RectBivariate' in Eddy.INTERP_METHOD:
@@ -573,23 +573,13 @@ def collection_loop(CS, grd, rtime, A_list_obj, C_list_obj,
573573
# Prepare for shape test and get eddy_radius_e
574574
#cx, cy = Eddy.M(contlon_e, contlat_e)
575575

576-
##################0
577576
#proj = pyproj.Proj(proj='aeqd', lat_0=contlat_e.mean(), lon_0=contlon_e.mean())#, width=5.e2, height=5.e2)
578577
# http://www.geo.hunter.cuny.edu/~jochen/gtech201/lectures/lec6concepts/map%20coordinate%20systems/how%20to%20choose%20a%20projection.htm
579578
proj = pyproj.Proj('+proj=aeqd +lat_0=%s +lon_0=%s'
580579
% (contlat_e.mean(), contlon_e.mean()))
581580
cx, cy = proj(contlon_e, contlat_e)
582581
centlon_e, centlat_e, eddy_radius_e, aerr = fit_circle(cx, cy)
583-
#print centlon_e, centlat_e, eddy_radius_e, aerr
584-
#print cxx, cyy
585-
##################0
586-
587-
#centlon_e, centlat_e, eddy_radius_e, aerr = fit_circle(cx, cy)
588-
#print centlon_e, centlat_e, eddy_radius_e, aerr
589-
#print cx, cy
590-
#print '-----------------------------------------------------'
591-
#aaaaaaaaaaaaaaaaa
592-
582+
593583
aerr = np.atleast_1d(aerr)
594584

595585
# Filter for shape: >35% (>55%) is not an eddy for Q (SLA)
@@ -668,13 +658,15 @@ def collection_loop(CS, grd, rtime, A_list_obj, C_list_obj,
668658
Eddy.jmin, Eddy.jmax)
669659

670660
# Get eddy circumference using eddy_radius_e
671-
centx, centy = Eddy.M(centlon_e, centlat_e)
672-
circlon, circlat = get_circle(centx, centy,
661+
#centx, centy = Eddy.M(centlon_e, centlat_e)
662+
#circlon, circlat = get_circle(centx, centy,
663+
#eddy_radius_e, 180)
664+
circlon, circlat = get_circle(centlon_e, centlat_e,
673665
eddy_radius_e, 180)
674666
#circlon[:], circlat[:] = Eddy.M.projtran(
675667
#circlon, circlat, inverse=True)
676-
circlon[:], circlat[:] = proj(circlon, circlat,
677-
inverse=True)
668+
#circlon[:], circlat[:] = proj(circlon, circlat,
669+
#inverse=True)
678670
Eddy.circlon, Eddy.circlat = circlon, circlat
679671

680672
# Set masked points within bounding box around eddy
@@ -690,7 +682,7 @@ def collection_loop(CS, grd, rtime, A_list_obj, C_list_obj,
690682
# circumferential distribution
691683
contlon_e, contlat_e = \
692684
eddy_tracker.uniform_resample(
693-
contlon_e, contlat_e, method='akima')
685+
contlon_e, contlat_e)
694686

695687
if 'Q' in Eddy.DIAGNOSTIC_TYPE:
696688
# Note, eddy amplitude == max(abs(vort/f)) within eddy, KCCMC11
@@ -963,6 +955,7 @@ def track_eddies(Eddy, first_record):
963955
# Make an ellipse at current old_eddy location
964956
# (See CSS11 sec. B4, pg. 208)
965957
if 'ellipse' in Eddy.SEPARATION_METHOD:
958+
966959
Eddy.search_ellipse.set_search_ellipse(old_x[old_ind],
967960
old_y[old_ind])
968961
#Eddy.search_ellipse.view_search_ellipse(Eddy)
@@ -1388,30 +1381,51 @@ def hann2d_fast(var, ni, nj):
13881381
return var
13891382

13901383

1391-
def get_circle(x0, y0, r, npts):
1392-
"""
1393-
Return points on a circle, with specified (x0, y0) center and radius
1394-
(and optional number of points too).
1384+
#def get_circle_OLD(x0, y0, r, npts):
1385+
#"""
1386+
#Return points on a circle, with specified (x0, y0) center and radius
1387+
#(and optional number of points too).
13951388

1396-
Input : 1 - x0, scalar, center X of circle
1397-
2 - y0, scalar, center Y of circle
1398-
3 - r, scalar, radius
1399-
4 - npts, scalar, number of points (optional)
1389+
#Input : 1 - x0, scalar, center X of circle
1390+
#2 - y0, scalar, center Y of circle
1391+
#3 - r, scalar, radius
1392+
#4 - npts, scalar, number of points (optional)
14001393

1401-
Output : 1 - cx, circle x-points
1402-
2 - cy, circle y-points
1394+
#Output : 1 - cx, circle x-points
1395+
#2 - cy, circle y-points
14031396

1404-
Example : [cx cy] = func_get_circle (5, 5, 3, 256)
1405-
plot (cx, cy, '-')
1397+
#Example : [cx cy] = func_get_circle (5, 5, 3, 256)
1398+
#plot (cx, cy, '-')
14061399

1407-
Written By : UCLA ROMS Team ([email protected])
1408-
Written On : June/05/2008
1409-
Tool : Eddy Tracker
1410-
"""
1411-
theta = np.arange(npts) # NOTE npts is a constant, so
1412-
# *cos(theta)* and *sin(theta)* can be predefined
1400+
#Written By : UCLA ROMS Team ([email protected])
1401+
#Written On : June/05/2008
1402+
#Tool : Eddy Tracker
1403+
#"""
1404+
#theta = np.arange(np.float64(npts)) # NOTE npts is a constant, so
1405+
## *cos(theta)* and *sin(theta)* can be predefined
1406+
## SHOULD BE PART OF CONTOUR OBJECT
1407+
#theta[:] = theta * 2. * (4. * np.arctan(1.)) / npts
1408+
#cx = x0 + r * np.cos(theta)
1409+
#cy = y0 + r * np.sin(theta)
1410+
#return cx, cy
1411+
1412+
1413+
def get_circle(lon, lat, distance, npts):
1414+
'''
1415+
Given the inputs (base lon, base lat, angle [deg], distance [m]) return
1416+
the lon, lat of the new position...
14131417
# SHOULD BE PART OF CONTOUR OBJECT
1414-
theta[:] = theta * 2. * (4. * np.arctan(1.)) / npts
1415-
cx = x0 + r * np.cos(theta)
1416-
cy = y0 + r * np.sin(theta)
1417-
return cx, cy
1418+
Could be speeded up with f2py
1419+
'''
1420+
pio180 = np.pi / 180.
1421+
lon *= pio180
1422+
lat *= pio180
1423+
angle = np.arange(np.float64(npts))
1424+
angle[:] = angle * 2. * (4. * np.arctan(1.)) / npts
1425+
distance /= 6371315.0 # angular distance
1426+
lat1 = (np.arcsin(np.sin(lat) * np.cos(distance) + np.cos(lat) *
1427+
np.sin(distance) * np.cos(angle)))
1428+
lon1 = (lon + np.arctan2(np.sin(angle) * np.sin(distance) * np.cos(lat),
1429+
np.cos(distance) - np.sin(lat) * np.sin(lat1)))
1430+
return lon1 / pio180, lat1 / pio180
1431+

0 commit comments

Comments
 (0)