Skip to content

Commit 3643eae

Browse files
committed
Radius calculations based on pyproj Azimuthal Equidistant projection
1 parent dda01ce commit 3643eae

File tree

3 files changed

+42
-15
lines changed

3 files changed

+42
-15
lines changed

make_eddy_track_AVISO.py

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -499,7 +499,10 @@ def __init__(self, AVISO_FILE, THE_DOMAIN, PRODUCT,
499499
super(AvisoGrid, self).__init__()
500500
print '\nInitialising the *AVISO_grid*'
501501
self.THE_DOMAIN = THE_DOMAIN
502-
self.PRODUCT = PRODUCT
502+
if 'AVISO' in PRODUCT:
503+
self.PRODUCT = 'AVISO'
504+
else:
505+
self.PRODUCT = PRODUCT
503506
self.LONMIN = LONMIN
504507
self.LONMAX = LONMAX
505508
self.LATMIN = LATMIN
@@ -531,7 +534,7 @@ def __init__(self, AVISO_FILE, THE_DOMAIN, PRODUCT,
531534
self.sla_coeffs = None
532535
self.uspd_coeffs = None
533536

534-
self.set_initial_indices(LONMIN, LONMAX, LATMIN, LATMAX)
537+
self.set_initial_indices()
535538
self.set_index_padding()
536539
self.set_basemap(with_pad=with_pad)
537540
self.get_AVISO_f_pm_pn()

make_eddy_track_ROMS.py

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -62,8 +62,7 @@
6262
#SAVE_DIR = DIRECTORY
6363
SAVE_DIR = '/marula/emason/runs2009/na_2009_7pt5km/eddy_tracking_exps/sla_1/'
6464

65-
# True to save outputs in same style as CSS11
66-
chelton_style_nc = True
65+
6766

6867
print 'SAVE_DIR', SAVE_DIR
6968

py_eddy_tracker_classes.py

Lines changed: 36 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -39,12 +39,12 @@
3939
import matplotlib.dates as dt
4040
import matplotlib.path as path
4141
import matplotlib.patches as patch
42-
42+
import mpl_toolkits.basemap.pyproj as pyproj
4343
import make_eddy_tracker_list_obj as eddy_tracker
4444
from py_eddy_tracker_property_classes import Amplitude, EddyProperty, interpolate
45-
#import haversine_distmat as hav # needs compiling with f2py
4645
from haversine import haversine # needs compiling with f2py
4746

47+
4848
def datestr2datetime(datestr):
4949
"""
5050
Take strings with format YYYYMMDD and convert to datetime instance
@@ -414,6 +414,7 @@ def get_uavg(Eddy, CS, collind, centlon_e, centlat_e, poly_eff,
414414

415415
theseglon, theseglat = eddy_tracker.uniform_resample(
416416
theseglon, theseglat, method='akima')
417+
417418
if 'RectBivariate' in Eddy.INTERP_METHOD:
418419
uavg = Eddy.uspd_coeffs.ev(theseglat[1:], theseglon[1:]).mean()
419420

@@ -499,9 +500,14 @@ def get_uavg(Eddy, CS, collind, centlon_e, centlat_e, poly_eff,
499500
citer.iternext()
500501

501502
if any_inner_contours: # set speed based contour parameters
502-
cx, cy = Eddy.M(theseglon, theseglat)
503+
proj = pyproj.Proj('+proj=aeqd +lat_0=%s +lon_0=%s'
504+
% (theseglat.mean(), theseglon.mean()))
505+
cx, cy = proj(theseglon, theseglat)
506+
#cx, cy = Eddy.M(theseglon, theseglat)
503507
centx_s, centy_s, eddy_radius_s, junk = fit_circle(cx, cy)
504-
centlon_s, centlat_s = Eddy.M.projtran(centx_s, centy_s, inverse=True)
508+
centlon_s, centlat_s = proj(centx_s, centy_s, inverse=True)
509+
#print 'fffffffffff', centlon_s, centlat_s
510+
#centlon_s, centlat_s = Eddy.M.projtran(centx_s, centy_s, inverse=True)
505511

506512
else: # use the effective contour
507513
centlon_s, centlat_s = centlon_e, centlat_e
@@ -565,18 +571,35 @@ def collection_loop(CS, grd, rtime, A_list_obj, C_list_obj,
565571
properties = EddyProperty()
566572

567573
# Prepare for shape test and get eddy_radius_e
568-
cx, cy = Eddy.M(contlon_e, contlat_e)
574+
#cx, cy = Eddy.M(contlon_e, contlat_e)
575+
576+
##################0
577+
#proj = pyproj.Proj(proj='aeqd', lat_0=contlat_e.mean(), lon_0=contlon_e.mean())#, width=5.e2, height=5.e2)
578+
# http://www.geo.hunter.cuny.edu/~jochen/gtech201/lectures/lec6concepts/map%20coordinate%20systems/how%20to%20choose%20a%20projection.htm
579+
proj = pyproj.Proj('+proj=aeqd +lat_0=%s +lon_0=%s'
580+
% (contlat_e.mean(), contlon_e.mean()))
581+
cx, cy = proj(contlon_e, contlat_e)
569582
centlon_e, centlat_e, eddy_radius_e, aerr = fit_circle(cx, cy)
570-
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+
571593
aerr = np.atleast_1d(aerr)
572594

573595
# Filter for shape: >35% (>55%) is not an eddy for Q (SLA)
574596
if aerr >= 0. and aerr <= Eddy.SHAPE_ERROR[collind]:
575597

576598
# Get centroid in lon lat
577-
centlon_e, centlat_e = Eddy.M.projtran(centlon_e,
578-
centlat_e,
579-
inverse=True)
599+
centlon_e, centlat_e = proj(centlon_e, centlat_e, inverse=True)
600+
#centlon_e, centlat_e = Eddy.M.projtran(centlon_e,
601+
#centlat_e,
602+
#inverse=True)
580603

581604
# For some reason centlat_e is transformed
582605
# by projtran to 'float'...
@@ -648,8 +671,10 @@ def collection_loop(CS, grd, rtime, A_list_obj, C_list_obj,
648671
centx, centy = Eddy.M(centlon_e, centlat_e)
649672
circlon, circlat = get_circle(centx, centy,
650673
eddy_radius_e, 180)
651-
circlon[:], circlat[:] = Eddy.M.projtran(
652-
circlon, circlat, inverse=True)
674+
#circlon[:], circlat[:] = Eddy.M.projtran(
675+
#circlon, circlat, inverse=True)
676+
circlon[:], circlat[:] = proj(circlon, circlat,
677+
inverse=True)
653678
Eddy.circlon, Eddy.circlat = circlon, circlat
654679

655680
# Set masked points within bounding box around eddy

0 commit comments

Comments
 (0)