Skip to content

Commit 340d184

Browse files
evanmasonevanmason
committed
Minor changes
1 parent 02d635a commit 340d184

File tree

3 files changed

+105
-72
lines changed

3 files changed

+105
-72
lines changed

‎make_eddy_track_AVISO.py‎

Lines changed: 53 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -31,6 +31,19 @@
3131
===========================================================================
3232
"""
3333

34+
"""
35+
WARNINGS...?
36+
37+
--- AVISO_file: /marula/emason/data/altimetry/MSLA/GLOBAL/DT/REF/dt_ref_global_merged_msla_h_qd_19930915_19930915_20100503.nc
38+
/usr/lib64/python2.7/site-packages/numpy/lib/arraysetops.py:197: RuntimeWarning: invalid value encountered in not_equal
39+
flag = np.concatenate(([True], aux[1:] != aux[:-1]))
40+
--- AVISO_file: /marula/emason/data/altimetry/MSLA/GLOBAL/DT/REF/dt_ref_global_merged_msla_h_qd_19930922_19930922_20100503.nc
41+
"""
42+
43+
44+
45+
46+
3447

3548
import glob as glob
3649
#import matplotlib.pyplot as plt
@@ -114,6 +127,7 @@ def kdt(lon, lat, limits, k=4):
114127
return iind, jind
115128

116129
if 'AvisoGrid' in self.__class__.__name__:
130+
117131
if self.zero_crossing is True:
118132
'''
119133
Used for a zero crossing, e.g., across Agulhas region
@@ -135,6 +149,7 @@ def half_limits(lon, lat):
135149
limits = half_limits(lon, lat)
136150
iind, jind = kdt(self._lon, self._lat, limits)
137151
self.i0 = iind.max()
152+
138153
return self
139154

140155

@@ -395,10 +410,19 @@ def get_geostrophic_velocity(self, zeta):
395410
surface from variables f, zeta, pm, pn...
396411
Note: output at rho points
397412
'''
398-
self.upad[:] = -self.gof() * self.v2rho_2d(self.vmask() * (zeta.data[1:] - zeta.data[:-1]) \
399-
* 0.5 * (self.pn()[1:] + self.pn()[:-1]))
400-
self.vpad[:] = self.gof() * self.u2rho_2d(self.umask() * (zeta.data[:, 1:] - zeta.data[:, :-1]) \
401-
* 0.5 * (self.pm()[:, 1:] + self.pm()[:, :-1]))
413+
gof = self.gof().view()
414+
vmask = self.v2rho_2d(self.vmask().view()
415+
zeta1, zeta2 = zeta.data[1:].view(), zeta.data[:-1].view()
416+
pn1, pn2 = self.pn()[1:].view(), self.pn()[:-1].view()
417+
self.upad[:] = ne.evaluate('-gof * vmask * (zeta1 - zeta2) * 0.5 * (pn1 + pn2)')
418+
#self.upad[:] = -self.gof() * self.v2rho_2d(self.vmask() * (zeta.data[1:] - zeta.data[:-1]) \
419+
#* 0.5 * (self.pn()[1:] + self.pn()[:-1]))
420+
umask = self.u2rho_2d(self.umask().view()
421+
zeta1, zeta2 = zeta.data[:,1:].view(), zeta.data[:,:-1].view()
422+
pm1, pm2 = self.pm()[:,1:].view(), self.pm()[:,:-1].view()
423+
self.vpad[:] = ne.evaluate('gof * umask * (zeta1 - zeta2) * 0.5 * (pm1 + pm2)')
424+
#self.vpad[:] = self.gof() * self.u2rho_2d(self.umask() * (zeta.data[:, 1:] - zeta.data[:, :-1]) \
425+
#* 0.5 * (self.pm()[:, 1:] + self.pm()[:, :-1]))
402426
return self
403427

404428

@@ -424,15 +448,17 @@ def getEKE(self):
424448
'''
425449
self.u[:] = self.upad[self.jup0:self.jup1, self.iup0:self.iup1]
426450
self.v[:] = self.vpad[self.jup0:self.jup1, self.iup0:self.iup1]
427-
np.add(self.u**2, self.v**2, out=self.eke)
428-
self.eke *= 0.5
451+
u, v = self.u.view(), self.v.view()
452+
self.eke[:] = ne.evaluate('0.5 * (u**2 + v**2)')
453+
#np.add(self.u**2, self.v**2, out=self.eke)
454+
#self.eke *= 0.5
429455
return self
430456

431457

432458

433459
class AvisoGrid (PyEddyTracker):
434460
'''
435-
Class to satisfy the need of the ROMS eddy tracker
461+
Class to satisfy the need of the eddy tracker
436462
to have a grid class
437463
'''
438464
def __init__(self, AVISO_file, lonmin, lonmax, latmin, latmax, with_pad=True, use_maskoceans=False):
@@ -705,7 +731,7 @@ def pcol_2dxy(self, x, y):
705731
#the_domain = 'MedSea' # not yet implemented
706732

707733
# Specify use of new AVISO 2014 data
708-
new_AVISO = True
734+
new_AVISO = False
709735

710736
# Specify subsampling new AVISO 2014 data; i.e. you may
711737
# prefer to use every second day rather than every day
@@ -730,8 +756,8 @@ def pcol_2dxy(self, x, y):
730756
directory = '/marula/emason/data/altimetry/global/delayed-time/grids/msla/two-sat-merged/h/'
731757
AVISO_files = 'dt_global_twosat_msla_h_????????_20140106.nc'
732758
else:
733-
directory = '/path/to/your/aviso_data/'
734-
#directory = '/marula/emason/data/altimetry/MSLA/GLOBAL/DT/REF/'
759+
#directory = '/path/to/your/aviso_data/'
760+
directory = '/marula/emason/data/altimetry/MSLA/GLOBAL/DT/REF/'
735761
AVISO_files = 'dt_ref_global_merged_msla_h_qd_????????_*.nc'
736762

737763
elif 'MedSea' in the_domain:
@@ -754,20 +780,21 @@ def pcol_2dxy(self, x, y):
754780
# Set date range (YYYYMMDD)
755781
#date_str, date_end = 19980107, 19991110 #
756782
#date_str, date_end = 20081107, 20100110 #
757-
date_str, date_end = 19930101, 20121231 #
783+
date_str, date_end = 19921014, 20120718 #
758784

759785
# Choose type of diagnostic: either q-parameter ('Q') or sea level anomaly ('sla')
760-
diag_type = 'Q' #<<< not implemented in 1.2.0
761-
#diag_type = 'sla'
786+
#diag_type = 'Q' #<<< not implemented in 1.2.0
787+
diag_type = 'sla'
762788

763789

764790
# Path to directory where outputs are to be saved...
765791
#savedir = directory
766792
#savedir = '/marula/emason/aviso_eddy_tracking/pablo_exp/'
767793
#savedir = '/marula/emason/aviso_eddy_tracking/new_AVISO_test/'
768794
#savedir = '/marula/emason/aviso_eddy_tracking/new_AVISO_SUBSAMP-3days/'
769-
savedir = '/marula/emason/aviso_eddy_tracking/junk/'
795+
#savedir = '/marula/emason/aviso_eddy_tracking/junk/'
770796
#savedir = '/marula/emason/aviso_eddy_tracking/new_AVISO_test/BlackSea/'
797+
savedir = '/marula/emason/aviso_eddy_tracking/Corridor_V3_Dec2014/'
771798
#savedir = '/path/to/save/your/outputs/'
772799

773800

@@ -810,8 +837,8 @@ def pcol_2dxy(self, x, y):
810837

811838
# Set SLA contour spacing
812839
elif 'sla' in diag_type:
813-
#slaparameter = np.arange(-100., 101., 1.0) # cm
814-
slaparameter = np.arange(-100., 100.25, 0.25) # cm
840+
slaparameter = np.arange(-100., 101., 1.0) # cm
841+
#slaparameter = np.arange(-100., 100.25, 0.25) # cm
815842

816843

817844

@@ -844,10 +871,10 @@ def pcol_2dxy(self, x, y):
844871
#latmin = -35.
845872
#latmax = -10.
846873

847-
#lonmin = -65. # Corridor
848-
#lonmax = -5.5
849-
#latmin = 11.5
850-
#latmax = 38.5
874+
lonmin = -65. # Corridor
875+
lonmax = -5.5
876+
latmin = 11.5
877+
latmax = 38.5
851878

852879
#lonmin = -179. # SEP
853880
#lonmax = -65
@@ -938,7 +965,7 @@ def pcol_2dxy(self, x, y):
938965

939966
# Use this for subsampling to get identical list as old_AVISO
940967
#AVISO_files = AVISO_files[5:-5:7]
941-
if new_AVISO_SUBSAMP:
968+
if new_AVISO and new_AVISO_SUBSAMP:
942969
AVISO_files = AVISO_files[5:-5:np.int(days_btwn_recs)]
943970

944971
# Set up a grid object using first AVISO file in the list
@@ -1341,7 +1368,8 @@ def pcol_2dxy(self, x, y):
13411368

13421369
if anim_figs: # Make figures for animations
13431370

1344-
tit = 'Y' + str(yr) + 'M' + str(mo).zfill(2) + 'D' + str(da).zfill(2)
1371+
#tit = 'Y' + str(yr) + 'M' + str(mo).zfill(2) + 'D' + str(da).zfill(2)
1372+
tit = ''.join((str(yr), str(mo).zfill(2), str(da).zfill(2)))
13451373
#tit = str(yr) + str(mo).zfill(2) + str(da).zfill(2)
13461374

13471375
#if 'anim_fig' in locals():
@@ -1380,12 +1408,12 @@ def pcol_2dxy(self, x, y):
13801408
print '--- saving to nc', C_eddy.savedir
13811409
print '+++'
13821410
if chelton_style_nc: # Recommended
1383-
A_eddy.write2chelton_nc(rtime)
1384-
C_eddy.write2chelton_nc(rtime)
1411+
A_eddy.write2netcdf(rtime)
1412+
C_eddy.write2netcdf(rtime)
13851413
else:
13861414
A_eddy.write2nc(A_savefile, rtime)
13871415
C_eddy.write2nc(C_savefile, rtime)
1388-
1416+
raise Exception, 'not supported'
13891417
#print 'Saving the eddies', time.time() - saving_start_time, 'seconds'
13901418

13911419
# Running time for a single monthly file

‎make_eddy_tracker_list_obj.py‎

Lines changed: 15 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -109,6 +109,9 @@ def strcompare(str1, str2):
109109

110110

111111

112+
113+
114+
112115
class track (object):
113116
'''
114117
Class that holds eddy tracks and related info
@@ -371,24 +374,18 @@ def get_inactive_tracks(self, rtime):
371374

372375

373376
def create_netcdf(self, directory, savedir, title,
374-
grd=None,
375-
Ymin=None,
376-
Ymax=None,
377-
Mmin=None,
378-
Mmax=None,
379-
model=None,
380-
sigma_lev=None,
381-
rho_ntr=None):
377+
grd=None, Ymin=None, Ymax=None,
378+
Mmin=None, Mmax=None, model=None,
379+
sigma_lev=None, rho_ntr=None):
382380
'''
383381
Create netcdf file same style as Chelton etal (2011)
384382
'''
385-
#print 'yooooooooooooo'
386383
if not self.track_extra_variables:
387384
self.savedir = savedir
388385
else:
389386
self.savedir = savedir.replace('.nc', '_ARGO_enabled.nc')
390387
nc = netcdf.Dataset(self.savedir, 'w', format='NETCDF4')
391-
nc.title = title + ' eddy tracks'
388+
nc.title = ''.join((title, ' eddy tracks'))
392389
nc.directory = directory
393390
nc.days_between_records = np.float64(self.days_btwn_recs)
394391
nc.track_duration_min = np.float64(self.track_duration_min)
@@ -615,7 +612,7 @@ def _reduce_inactive_tracks(self):
615612
#return np.array([self.tracklist[i].Uavg]) / distances
616613

617614

618-
def write2chelton_nc(self, rtime):
615+
def write2netcdf(self, rtime):
619616
'''
620617
Write inactive tracks to netcdf file.
621618
'ncind' is important because prevents writing of
@@ -790,19 +787,19 @@ def insert_at_index(self, xarr, ind, x):
790787
newsize = tmp.size
791788
else:
792789
newsize = ind + 1
793-
#print 'numpy-------'
794790

795791
elif isinstance(tmp, list):
796792
tmp = list(tmp)
797793
if ind < len(tmp):
798794
newsize = len(tmp)
799795
else:
800796
newsize = ind + 1
801-
#print 'list-------'
797+
802798

803799
else:
804800
Exception
805801

802+
# First, numpy arrays...
806803
if strcompare('new_lon', xarr):
807804
self.new_lon = np.zeros((newsize))
808805
self.new_lon[:tmp.size] = tmp
@@ -852,23 +849,25 @@ def insert_at_index(self, xarr, ind, x):
852849
self.new_bounds[:tmp.shape[0]] = tmp
853850
self.new_bounds[ind] = x
854851

852+
# Second, lists...
855853
elif strcompare('new_contour_e', xarr):
856854
try:
857855
self.new_contour_e[ind] = x
858856
except:
859-
self.new_contour_e += [0] * (ind - len(self.new_contour_e) + 1)
857+
#self.new_contour_e += [0] * (ind - len(self.new_contour_e) + 1)
858+
self.new_contour_e.append([0] * (ind - len(self.new_contour_e) + 1))
860859
self.new_contour_e[ind] = x
861860
elif strcompare('new_contour_s', xarr):
862861
try:
863862
self.new_contour_s[ind] = x
864863
except:
865-
self.new_contour_s += [0] * (ind - len(self.new_contour_s) + 1)
864+
self.new_contour_s.append([0] * (ind - len(self.new_contour_s) + 1))
866865
self.new_contour_s[ind] = x
867866
elif strcompare('new_Uavg_profile', xarr):
868867
try:
869868
self.new_Uavg_profile[ind] = x
870869
except:
871-
self.new_Uavg_profile += [0] * (ind - len(self.new_Uavg_profile) + 1)
870+
self.new_Uavg_profile.append([0] * (ind - len(self.new_Uavg_profile) + 1))
872871
self.new_Uavg_profile[ind] = x
873872
else:
874873
raise Exception

0 commit comments

Comments
 (0)