Skip to content

Commit a56f2fc

Browse files
evanmasonevanmason
committed
Added option for RectBivariate or griddata
1 parent 2a5ce4a commit a56f2fc

File tree

5 files changed

+220
-213
lines changed

5 files changed

+220
-213
lines changed

‎eddy_tracker_configuration.yaml‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -118,7 +118,8 @@ SAVE_FIGURES: Yes
118118
# Useful for working with ARGO data
119119
TRACK_EXTRA_VARIABLES: No
120120

121-
121+
# INTERP_METHOD: 'RectBivariate'
122+
INTERP_METHOD: 'griddata'
122123

123124

124125

‎make_eddy_track_AVISO.py‎

Lines changed: 57 additions & 141 deletions
Original file line numberDiff line numberDiff line change
@@ -504,15 +504,15 @@ def __init__(self, AVISO_FILE, LONMIN, LONMAX, LATMIN, LATMAX,
504504
try: # new AVISO (2014)
505505
self._lon = self.read_nc(AVISO_FILE, 'lon')
506506
self._lat = self.read_nc(AVISO_FILE, 'lat')
507-
self.fillval = self.read_nc_att(AVISO_FILE, 'sla', '_FillValue')
507+
self.FILLVAL = self.read_nc_att(AVISO_FILE, 'sla', '_FillValue')
508508
base_date = self.read_nc_att(AVISO_FILE, 'time', 'units')
509509
self.base_date = dt.date2num(
510510
parser.parse(base_date.split(' ')[2:4][0]))
511511

512512
except Exception: # old AVISO
513513
self._lon = self.read_nc(AVISO_FILE, 'NbLongitudes')
514514
self._lat = self.read_nc(AVISO_FILE, 'NbLatitudes')
515-
self.fillval = self.read_nc_att(AVISO_FILE,
515+
self.FILLVAL = self.read_nc_att(AVISO_FILE,
516516
'Grid_0001', '_FillValue')
517517

518518
if LONMIN < 0 and LONMAX <= 0:
@@ -765,13 +765,13 @@ def set_interp_coeffs(self, sla, uspd):
765765

766766
DIAGNOSTIC_TYPE = config['DIAGNOSTIC_TYPE']
767767

768-
THE_DOMAIN = config['DOMAIN']['THE_DOMAIN']
769-
LONMIN = config['DOMAIN']['LONMIN']
770-
LONMAX = config['DOMAIN']['LONMAX']
771-
LATMIN = config['DOMAIN']['LATMIN']
772-
LATMAX = config['DOMAIN']['LATMAX']
773-
DATE_STR = config['DOMAIN']['DATE_STR']
774-
DATE_END = config['DOMAIN']['DATE_END']
768+
config['THE_DOMAIN'] = config['DOMAIN']['THE_DOMAIN']
769+
config['LONMIN'] = config['DOMAIN']['LONMIN']
770+
config['LONMAX'] = config['DOMAIN']['LONMAX']
771+
config['LATMIN'] = config['DOMAIN']['LATMIN']
772+
config['LATMAX'] = config['DOMAIN']['LATMAX']
773+
DATE_STR = config['DATE_STR'] = config['DOMAIN']['DATE_STR']
774+
DATE_END = config['DATE_END'] = config['DOMAIN']['DATE_END']
775775

776776
AVISO_DT14 = config['AVISO']['AVISO_DT14']
777777
AVISO_FILES = config['AVISO']['AVISO_FILES']
@@ -784,31 +784,35 @@ def set_interp_coeffs(self, sla, uspd):
784784
else:
785785
DAYS_BTWN_RECORDS = 7. # old seven day AVISO
786786

787-
TRACK_DURATION_MIN = config['TRACK_DURATION_MIN']
787+
#TRACK_DURATION_MIN = config['TRACK_DURATION_MIN']
788788

789789
if 'SLA' in DIAGNOSTIC_TYPE:
790790
MAX_SLA = config['CONTOUR_PARAMETER']['CONTOUR_PARAMETER_SLA']['MAX_SLA']
791791
INTERVAL = config['CONTOUR_PARAMETER']['CONTOUR_PARAMETER_SLA']['INTERVAL']
792-
CONTOUR_PARAMETER = np.arange(-MAX_SLA, MAX_SLA + INTERVAL, INTERVAL)
793-
SHAPE_ERROR = config['SHAPE_ERROR'] * np.ones(CONTOUR_PARAMETER.size)
792+
config['CONTOUR_PARAMETER'] = np.arange(-MAX_SLA, MAX_SLA + INTERVAL, INTERVAL)
793+
config['SHAPE_ERROR'] = np.full(config['CONTOUR_PARAMETER'].size,
794+
config['SHAPE_ERROR'])
795+
794796
elif 'Q' in DIAGNOSTIC_TYPE:
795797
MAX_Q = config['CONTOUR_PARAMETER']['CONTOUR_PARAMETER_Q']['MAX_Q']
796798
NUM_LEVS = config['CONTOUR_PARAMETER']['CONTOUR_PARAMETER_Q']['NUM_LEVS']
797-
CONTOUR_PARAMETER = np.linspace(0, MAX_Q, NUM_LEVS)[::-1]
799+
config['CONTOUR_PARAMETER'] = np.linspace(0, MAX_Q, NUM_LEVS)[::-1]
798800
else: Exception
799801

800-
JDAY_REFERENCE = config['JDAY_REFERENCE']
802+
#JDAY_REFERENCE = config['JDAY_REFERENCE']
801803

802-
RADMIN = config['RADMIN']
803-
RADMAX = config['RADMAX']
804+
#RADMIN = config['RADMIN']
805+
#RADMAX = config['RADMAX']
804806

805807
if 'SLA' in DIAGNOSTIC_TYPE:
806-
AMPMIN = config['AMPMIN']
807-
AMPMAX = config['AMPMAX']
808+
#AMPMIN = config['AMPMIN']
809+
#AMPMAX = config['AMPMAX']
810+
pass
808811
elif 'Q' in DIAGNOSTIC_TYPE:
809812
AMPMIN = 0.02 # max(abs(xi/f)) within the eddy
810813
AMPMAX = 100.
811-
else: Exception
814+
else:
815+
Exception
812816

813817
SAVE_FIGURES = config['SAVE_FIGURES']
814818

@@ -823,27 +827,22 @@ def set_interp_coeffs(self, sla, uspd):
823827
SMOOTHING_TYPE = config['SMOOTHING_SLA']['TYPE']
824828
else: Exception
825829

826-
DIST0 = config['DIST0']
827-
AREA0 = np.pi * config['RAD0']**2
828830
if 'Q' in DIAGNOSTIC_TYPE:
829831
AMP0 = 0.02 # vort/f
830832
elif 'SLA' in DIAGNOSTIC_TYPE:
831833
AMP0 = config['AMP0']
832834
TEMP0 = config['TEMP0']
833835
SALT0 = config['SALT0']
834836

835-
EVOLVE_AMP_MIN = config['EVOLVE_AMP_MIN']
836-
EVOLVE_AMP_MAX = config['EVOLVE_AMP_MAX']
837-
EVOLVE_AREA_MIN = config['EVOLVE_AREA_MIN']
838-
EVOLVE_AREA_MAX = config['EVOLVE_AREA_MAX']
837+
#EVOLVE_AMP_MIN = config['EVOLVE_AMP_MIN']
838+
#EVOLVE_AMP_MAX = config['EVOLVE_AMP_MAX']
839+
#EVOLVE_AREA_MIN = config['EVOLVE_AREA_MIN']
840+
#EVOLVE_AREA_MAX = config['EVOLVE_AREA_MAX']
841+
839842

840-
SEPARATION_METHOD = config['SEPARATION_METHOD']
841843

842-
MAX_LOCAL_EXTREMA = config['MAX_LOCAL_EXTREMA']
843844

844-
TRACK_EXTRA_VARIABLES = config['TRACK_EXTRA_VARIABLES']
845845

846-
VERBOSE = config['VERBOSE']
847846

848847
CMAP = plt.cm.RdBu
849848

@@ -866,23 +865,15 @@ def set_interp_coeffs(self, sla, uspd):
866865
AVISO_FILES = AVISO_FILES[5:-5:np.int(DAYS_BTWN_RECORDS)]
867866

868867
# Set up a grid object using first AVISO file in the list
869-
sla_grd = AvisoGrid(AVISO_FILES[0], LONMIN, LONMAX, LATMIN, LATMAX)
870-
868+
sla_grd = AvisoGrid(AVISO_FILES[0], config['LONMIN'], config['LONMAX'],
869+
config['LATMIN'], config['LATMAX'])
871870

872-
873-
874-
875871
Mx, My = (sla_grd.Mx[sla_grd.jup0:sla_grd.jup1, sla_grd.iup0:sla_grd.iup1],
876872
sla_grd.My[sla_grd.jup0:sla_grd.jup1, sla_grd.iup0:sla_grd.iup1])
877-
#pMx, pMy = sla_grd.pcol_2dxy(Mx, My)
878-
879-
880-
881-
882873

883874
# Instantiate search ellipse object
884-
search_ellipse = eddy_tracker.SearchEllipse(THE_DOMAIN, sla_grd,
885-
DAYS_BTWN_RECORDS, RW_PATH)
875+
search_ellipse = eddy_tracker.SearchEllipse(config['THE_DOMAIN'],
876+
sla_grd, DAYS_BTWN_RECORDS, RW_PATH)
886877

887878

888879
if 'Gaussian' in SMOOTHING_TYPE:
@@ -892,120 +883,42 @@ def set_interp_coeffs(self, sla, uspd):
892883

893884
fig = plt.figure(1)
894885

895-
# Initialise two eddy objects to hold data
896-
A_eddy = eddy_tracker.TrackList('AVISO', TRACK_DURATION_MIN,
897-
TRACK_EXTRA_VARIABLES)
898-
C_eddy = eddy_tracker.TrackList('AVISO', TRACK_DURATION_MIN,
899-
TRACK_EXTRA_VARIABLES)
900-
901886
if 'Q' in DIAGNOSTIC_TYPE:
902887
A_SAVEFILE = "".join([SAVE_DIR, 'eddy_tracks_Q_AVISO_anticyclonic.nc'])
903-
A_eddy.qparameter = qparameter
904-
C_eddy.qparameter = qparameter
905888
C_SAVEFILE = "".join([SAVE_DIR, 'eddy_tracks_Q_AVISO_cyclonic.nc'])
906-
A_eddy.SHAPE_ERROR = SHAPE_ERROR
907-
C_eddy.SHAPE_ERROR = SHAPE_ERROR
908889

909890
elif 'SLA' in DIAGNOSTIC_TYPE:
910891
A_SAVEFILE = "".join([SAVE_DIR, 'eddy_tracks_SLA_AVISO_anticyclonic.nc'])
911-
A_eddy.CONTOUR_PARAMETER = CONTOUR_PARAMETER
912-
A_eddy.SHAPE_ERROR = SHAPE_ERROR
913892
C_SAVEFILE = "".join([SAVE_DIR, 'eddy_tracks_SLA_AVISO_cyclonic.nc'])
914-
C_eddy.CONTOUR_PARAMETER = CONTOUR_PARAMETER[::-1]
915-
C_eddy.SHAPE_ERROR = SHAPE_ERROR[::-1]
916-
917-
A_eddy.JDAY_REFERENCE = JDAY_REFERENCE
918-
C_eddy.JDAY_REFERENCE = JDAY_REFERENCE
919-
920-
A_eddy.INTERANNUAL = True
921-
C_eddy.INTERANNUAL = True
922-
923-
A_eddy.DIAGNOSTIC_TYPE = DIAGNOSTIC_TYPE
924-
C_eddy.DIAGNOSTIC_TYPE = DIAGNOSTIC_TYPE
925-
926-
A_eddy.SMOOTHING = SMOOTHING
927-
C_eddy.SMOOTHING = SMOOTHING
928893

929-
A_eddy.MAX_LOCAL_EXTREMA = MAX_LOCAL_EXTREMA
930-
C_eddy.MAX_LOCAL_EXTREMA = MAX_LOCAL_EXTREMA
931894

932-
A_eddy.M = sla_grd.M
933-
C_eddy.M = sla_grd.M
934-
935-
#A_eddy.rwv = rwv
936-
#C_eddy.rwv = rwv
937895

896+
# Initialise two eddy objects to hold data
897+
#kwargs = config
898+
A_eddy = eddy_tracker.TrackList('AVISO', 'Anticyclonic', A_SAVEFILE,
899+
sla_grd, search_ellipse, **config)
900+
C_eddy = eddy_tracker.TrackList('AVISO', 'Cyclonic', C_SAVEFILE,
901+
sla_grd, search_ellipse, **config)
902+
938903
A_eddy.search_ellipse = search_ellipse
939904
C_eddy.search_ellipse = search_ellipse
940905

941-
942-
A_eddy.SEPARATION_METHOD = SEPARATION_METHOD
943-
C_eddy.SEPARATION_METHOD = SEPARATION_METHOD
944-
945-
if 'sum_radii' in SEPARATION_METHOD:
906+
if 'sum_radii' in config['SEPARATION_METHOD']:
946907
A_eddy.SEP_DIST_FAC = SEP_DIST_FACTOR
947908
C_eddy.SEP_DIST_FACTOR = SEP_DIST_FACTOR
948909

949-
950-
A_eddy.points = np.array([sla_grd.lon().ravel(),
951-
sla_grd.lat().ravel()]).T
952-
C_eddy.points = np.array([sla_grd.lon().ravel(),
953-
sla_grd.lat().ravel()]).T
954-
955-
A_eddy.EVOLVE_AMP_MIN = np.float64(EVOLVE_AMP_MIN)
956-
A_eddy.EVOLVE_AMP_MAX = np.float64(EVOLVE_AMP_MAX)
957-
A_eddy.EVOLVE_AREA_MIN = np.float64(EVOLVE_AREA_MIN)
958-
A_eddy.EVOLVE_AREA_MAX = np.float64(EVOLVE_AREA_MAX)
959-
960-
C_eddy.EVOLVE_AMP_MIN = np.float64(EVOLVE_AMP_MIN)
961-
C_eddy.EVOLVE_AMP_MAX = np.float64(EVOLVE_AMP_MAX)
962-
C_eddy.EVOLVE_AREA_MIN = np.float64(EVOLVE_AREA_MIN)
963-
C_eddy.EVOLVE_AREA_MAX = np.float64(EVOLVE_AREA_MAX)
964-
965-
A_eddy.i0, A_eddy.i1 = sla_grd.i0, sla_grd.i1
966-
A_eddy.j0, A_eddy.j1 = sla_grd.j0, sla_grd.j1
967-
C_eddy.i0, C_eddy.i1 = sla_grd.i0, sla_grd.i1
968-
C_eddy.j0, C_eddy.j1 = sla_grd.j0, sla_grd.j1
969-
970-
A_eddy.LONMIN, A_eddy.LONMAX = np.float64(LONMIN), np.float64(LONMAX)
971-
A_eddy.LATMIN, A_eddy.LATMAX = np.float64(LATMIN), np.float64(LATMAX)
972-
C_eddy.LONMIN, C_eddy.LONMAX = np.float64(LONMIN), np.float64(LONMAX)
973-
C_eddy.LATMIN, C_eddy.LATMAX = np.float64(LATMIN), np.float64(LATMAX)
974-
975-
A_eddy.RADMIN = np.float64(RADMIN)
976-
A_eddy.RADMAX = np.float64(RADMAX)
977-
A_eddy.AMPMIN = np.float64(AMPMIN)
978-
A_eddy.AMPMAX = np.float64(AMPMAX)
979-
C_eddy.RADMIN = np.float64(RADMIN)
980-
C_eddy.RADMAX = np.float64(RADMAX)
981-
C_eddy.AMPMIN = np.float64(AMPMIN)
982-
C_eddy.AMPMAX = np.float64(AMPMAX)
983-
984-
A_eddy.fillval = sla_grd.fillval
985-
C_eddy.fillval = sla_grd.fillval
986-
A_eddy.VERBOSE = VERBOSE
987-
C_eddy.VERBOSE = VERBOSE
988-
989-
990910
# See Chelton section B2 (0.4 degree radius)
991911
# These should give 8 and 1000 for 0.25 deg resolution
992-
PIXMIN = np.round((np.pi * RADMIN**2) / sla_grd.get_resolution()**2)
993-
PIXMAX = np.round((np.pi * RADMAX**2) / sla_grd.get_resolution()**2)
994-
print '--- Pixel range = %s-%s' % (np.int(PIXMIN), np.int(PIXMAX))
912+
PIXMIN = np.round((np.pi * config['RADMIN']**2) /
913+
sla_grd.get_resolution()**2)
914+
PIXMAX = np.round((np.pi * config['RADMAX']**2) /
915+
sla_grd.get_resolution()**2)
916+
print '--- Pixel range = %s-%s' % (np.int(PIXMIN),
917+
np.int(PIXMAX))
995918

996919
A_eddy.PIXEL_THRESHOLD = [PIXMIN, PIXMAX]
997920
C_eddy.PIXEL_THRESHOLD = [PIXMIN, PIXMAX]
998921

999-
A_eddy.AREA0 = np.float64(AREA0)
1000-
C_eddy.AREA0 = np.float64(AREA0)
1001-
A_eddy.AMP0 = np.float64(AMP0)
1002-
C_eddy.AMP0 = np.float64(AMP0)
1003-
A_eddy.DIST0 = np.float64(DIST0)
1004-
C_eddy.DIST0 = np.float64(DIST0)
1005-
1006-
A_eddy.DAYS_BTWN_RECORDS = DAYS_BTWN_RECORDS
1007-
C_eddy.DAYS_BTWN_RECORDS = DAYS_BTWN_RECORDS
1008-
1009922
# Create nc files for saving of eddy tracks
1010923
A_eddy.create_netcdf(DATA_DIR, A_SAVEFILE, 'Anticyclonic')
1011924
C_eddy.create_netcdf(DATA_DIR, C_SAVEFILE, 'Cyclonic')
@@ -1064,7 +977,7 @@ def set_interp_coeffs(self, sla, uspd):
1064977
print '------ applying Gaussian high-pass filter'
1065978
# Set landpoints to zero
1066979
np.place(sla, sla_grd.mask == False, 0.)
1067-
np.place(sla, sla.data == sla_grd.fillval, 0.)
980+
np.place(sla, sla.data == sla_grd.FILLVAL, 0.)
1068981
# High pass filter, see
1069982
# http://stackoverflow.com/questions/6094957/high-pass-filter-for-image-processing-in-python-by-using-scipy-numpy
1070983
sla -= ndimage.gaussian_filter(sla, [mres, zres])
@@ -1169,12 +1082,13 @@ def set_interp_coeffs(self, sla, uspd):
11691082

11701083
elif 'SLA' in DIAGNOSTIC_TYPE:
11711084
A_CS = ax.contour(sla_grd.lon(),
1172-
sla_grd.lat(), A_eddy.sla, CONTOUR_PARAMETER)
1085+
sla_grd.lat(),
1086+
A_eddy.sla, A_eddy.CONTOUR_PARAMETER)
11731087
# Note that CSc is for the cyclonics,
11741088
# CONTOUR_PARAMETER in reverse order
11751089
C_CS = ax.contour(sla_grd.lon(),
11761090
sla_grd.lat(),
1177-
C_eddy.sla, CONTOUR_PARAMETER[::-1])
1091+
C_eddy.sla, C_eddy.CONTOUR_PARAMETER)
11781092

11791093
else:
11801094
Exception
@@ -1199,8 +1113,10 @@ def set_interp_coeffs(self, sla, uspd):
11991113
C_eddy.swirl = SwirlSpeed(C_CS)
12001114

12011115
# Now we loop over the CS collection
1202-
A_eddy.sign_type = 'Anticyclonic'
1203-
C_eddy.sign_type = 'Cyclonic'
1116+
A_eddy.SIGN_TYPE = 'Anticyclonic'
1117+
C_eddy.SIGN_TYPE = 'Cyclonic'
1118+
1119+
qqq
12041120
if 'Q' in DIAGNOSTIC_TYPE:
12051121
A_eddy, C_eddy = collection_loop(CS, sla_grd, rtime,
12061122
A_list_obj=A_eddy, C_list_obj=C_eddy,
@@ -1209,11 +1125,11 @@ def set_interp_coeffs(self, sla, uspd):
12091125
elif 'SLA' in DIAGNOSTIC_TYPE:
12101126
A_eddy = collection_loop(A_CS, sla_grd, rtime,
12111127
A_list_obj=A_eddy, C_list_obj=None,
1212-
sign_type=A_eddy.sign_type, VERBOSE=VERBOSE)
1128+
sign_type=A_eddy.SIGN_TYPE, VERBOSE=A_eddy.VERBOSE)
12131129
# Note that C_CS is reverse order
12141130
C_eddy = collection_loop(C_CS, sla_grd, rtime,
12151131
A_list_obj=None, C_list_obj=C_eddy,
1216-
sign_type=C_eddy.sign_type, VERBOSE=VERBOSE)
1132+
sign_type=C_eddy.SIGN_TYPE, VERBOSE=C_eddy.VERBOSE)
12171133

12181134

12191135

@@ -1297,7 +1213,7 @@ def set_interp_coeffs(self, sla, uspd):
12971213
# IMPORTANT: this must be done at every time step!!
12981214
#saving_START_TIME = time.time()
12991215
if not first_record:
1300-
if VERBOSE:
1216+
if A_eddy.VERBOSE:
13011217
print '--- saving to nc', A_eddy.SAVE_DIR
13021218
print '--- saving to nc', C_eddy.SAVE_DIR
13031219
print '+++'

0 commit comments

Comments
 (0)