Skip to content

Commit 29fe893

Browse files
evanmasonevanmason
committed
Fix applied for local min/max problem
1 parent 620bc55 commit 29fe893

File tree

4 files changed

+188
-76
lines changed

4 files changed

+188
-76
lines changed

‎make_eddy_track_AVISO.py‎

Lines changed: 54 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -145,7 +145,7 @@ def set_initial_indices(self, LONMIN, LONMAX, LATMIN, LATMAX):
145145
"""
146146
Get indices for desired domain
147147
"""
148-
print '--- Setting initial indices to %s domain' % self.THE_DOMAIN
148+
print '--- Setting initial indices to *%s* domain' % self.THE_DOMAIN
149149
print '------ LONMIN = %s, LONMAX = %s, LATMIN = %s, LATMAX = %s' % (
150150
LONMIN, LONMAX, LATMIN, LATMAX)
151151
self.i0, junk = self.nearest_point(LONMIN,
@@ -408,27 +408,26 @@ def uvmask(self):
408408
return self
409409

410410

411-
def make_gridmask(self, with_pad=True):
411+
def set_basemap(self, with_pad=True):
412412
"""
413413
Use Basemap to make a landmask
414414
Format is 1 == ocean, 0 == land
415415
"""
416416
print '--- Computing Basemap'
417417
# Create Basemap instance for Mercator projection.
418-
self.M = Basemap(projection='merc', llcrnrlon = self.LONMIN - 1,
419-
urcrnrlon = self.LONMAX + 1,
420-
llcrnrlat = self.LATMIN - 1,
421-
urcrnrlat = self.LATMAX + 1,
422-
lat_ts = 0.5 * (self.LATMIN +
423-
self.LATMAX),
424-
resolution = 'h')
418+
self.M = Basemap(projection='merc',
419+
llcrnrlon = self.LONMIN - 1, urcrnrlon = self.LONMAX + 1,
420+
llcrnrlat = self.LATMIN - 1, urcrnrlat = self.LATMAX + 1,
421+
lat_ts = 0.5 * (self.LATMIN + self.LATMAX), resolution = 'h')
422+
425423
if with_pad:
426424
x, y = self.M(self.lonpad(), self.latpad())
427425
else:
428426
x, y = self.M(self.lon(), self.lat())
429427
self.Mx, self.My = x, y
430428
return self
431429

430+
432431

433432
#@timeit
434433
def set_geostrophic_velocity(self, zeta):
@@ -532,7 +531,7 @@ def __init__(self, AVISO_FILE, THE_DOMAIN, LONMIN, LONMAX, LATMIN, LATMAX,
532531

533532
self.set_initial_indices(LONMIN, LONMAX, LATMIN, LATMAX)
534533
self.set_index_padding()
535-
self.make_gridmask(with_pad)
534+
self.set_basemap(with_pad=with_pad)
536535
self.get_AVISO_f_pm_pn()
537536
self.set_u_v_eke()
538537
pad2 = 2 * self.pad
@@ -587,12 +586,38 @@ def set_mask(self, sla):
587586
"""
588587
if sla.mask.size == 1: # all sea points
589588
self.mask = np.ones_like(sla.data)
589+
590590
else:
591-
self.mask = sla.mask.astype(np.int) - 1
592-
self.mask = np.abs(self.mask)
591+
self.mask = np.logical_not(sla.mask).astype(np.int)
592+
593+
if 'Global' in self.THE_DOMAIN:
594+
595+
# Close Drake Passage
596+
minus70 = np.argmin(np.abs(self.lonpad()[0] + 70))
597+
self.mask[:125, minus70] = 0
598+
# Mask all unwanted regions (Caspian Sea, etc)
599+
labels = ndimage.label(self.mask)[0]
600+
plus200 = np.argmin(np.abs(self.lonpad()[0] - 200))
601+
plus10 = np.argmin(np.abs(self.latpad()[:, 0] - 10))
602+
# Set to known sea point
603+
good_lab = labels[plus10, plus200]
604+
self.mask[labels != good_lab] = 0
605+
593606
return self
594607

595608

609+
610+
def set_global_mask(self):
611+
"""
612+
"""
613+
if 'Global' in self.THE_DOMAIN:
614+
labels = ndimage.label(np.logical_not(self.sla.mask))[0]
615+
616+
else:
617+
return
618+
return
619+
620+
596621
def fillmask(self, x, mask):
597622
"""
598623
Fill missing values in an array with an average of nearest
@@ -992,23 +1017,24 @@ def set_interp_coeffs(self, sla, uspd):
9921017
if 'first_record' not in locals():
9931018
print '------ applying Gaussian high-pass filter'
9941019
# Set landpoints to zero
995-
np.place(sla, sla_grd.mask == False, 0.)
1020+
np.place(sla, sla_grd.mask == 0, 0.)
9961021
np.place(sla, sla.data == sla_grd.FILLVAL, 0.)
9971022
# High pass filter, see
9981023
# http://stackoverflow.com/questions/6094957/high-pass-filter-for-image-processing-in-python-by-using-scipy-numpy
9991024
sla -= ndimage.gaussian_filter(sla, [mres, zres])
10001025

10011026
elif 'Hanning' in SMOOTHING_TYPE:
10021027

1003-
print '------ applying %s passes of Hanning filter' \
1028+
if 'first_record' not in locals():
1029+
print '------ applying %s passes of Hanning filter' \
10041030
% SMOOTH_FAC
10051031
# Do SMOOTH_FAC passes of 2d Hanning filter
10061032
sla = func_hann2d_fast(sla, SMOOTH_FAC)
10071033

10081034
else: Exception
10091035

1010-
# Expand the landmask
1011-
sla = np.ma.masked_where(sla_grd.mask == False, sla)
1036+
# Apply the landmask
1037+
sla = np.ma.masked_where(sla_grd.mask == 0, sla)
10121038

10131039
# Get timing
10141040
try:
@@ -1045,7 +1071,7 @@ def set_interp_coeffs(self, sla, uspd):
10451071
sla_grd.iup0:sla_grd.iup1]
10461072
xi = np.ma.masked_where(sla_grd.mask[
10471073
sla_grd.jup0:sla_grd.jup1,
1048-
sla_grd.iup0:sla_grd.iup1] == False,
1074+
sla_grd.iup0:sla_grd.iup1] == 0,
10491075
xi / sla_grd.f()[sla_grd.jup0:sla_grd.jup1,
10501076
sla_grd.iup0:sla_grd.iup1])
10511077

@@ -1062,7 +1088,7 @@ def set_interp_coeffs(self, sla, uspd):
10621088
uspd = np.sqrt(sla_grd.u**2 + sla_grd.v**2)
10631089
uspd = np.ma.masked_where(
10641090
sla_grd.mask[sla_grd.jup0:sla_grd.jup1,
1065-
sla_grd.iup0:sla_grd.iup1] == False,
1091+
sla_grd.iup0:sla_grd.iup1] == 0,
10661092
uspd)
10671093
A_eddy.uspd = uspd.copy()
10681094
C_eddy.uspd = uspd.copy()
@@ -1144,10 +1170,15 @@ def set_interp_coeffs(self, sla, uspd):
11441170

11451171

11461172

1173+
1174+
ymd_str = ''.join((str(yr), str(mo).zfill(2), str(da).zfill(2)))
11471175

11481176

1149-
## Test pickling
1150-
#with open('C_eddy.pkl', 'wb') as save_pickle:
1177+
# Test pickling
1178+
#with open("".join((SAVE_DIR, 'A_eddy_%s.pkl' % ymd_str)), 'wb') as save_pickle:
1179+
#pickle.dump(A_eddy, save_pickle)
1180+
1181+
#with open("".join((SAVE_DIR, 'C_eddy_%s.pkl' % ymd_str)), 'wb') as save_pickle:
11511182
#pickle.dump(C_eddy, save_pickle)
11521183

11531184
## Unpickle
@@ -1213,17 +1244,18 @@ def set_interp_coeffs(self, sla, uspd):
12131244

12141245
if 'Q' in DIAGNOSTIC_TYPE:
12151246
anim_figure(A_eddy, C_eddy, Mx, My, plt.cm.RdBu_r, rtime, DIAGNOSTIC_TYPE,
1216-
SAVE_DIR, 'Q-parameter ' + tit, animax, animax_cbar,
1247+
SAVE_DIR, 'Q-parameter ' + ymd_str, animax, animax_cbar,
12171248
qparam=qparam, qparameter=qparameter, xi=xi, xicopy=xicopy)
12181249

12191250
elif 'SLA' in DIAGNOSTIC_TYPE:
12201251
anim_figure(A_eddy, C_eddy, Mx, My, plt.cm.RdBu_r, rtime, DIAGNOSTIC_TYPE,
1221-
SAVE_DIR, 'SLA ' + tit, animax, animax_cbar)
1252+
SAVE_DIR, 'SLA ' + ymd_str, animax, animax_cbar)
12221253

12231254
# Save inactive eddies to nc file
12241255
# IMPORTANT: this must be done at every time step!!
12251256
#saving_START_TIME = time.time()
12261257
if not first_record:
1258+
12271259
if A_eddy.VERBOSE:
12281260
print '--- saving to nc', A_eddy.SAVE_DIR
12291261
print '--- saving to nc', C_eddy.SAVE_DIR

‎make_eddy_tracker_list_obj.py‎

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -308,6 +308,7 @@ def __init__(self, DATATYPE, SIGN_TYPE, SAVE_DIR, grd, search_ellipse,
308308
# NOTE: '.copy()' suffix is essential here
309309
self.CONTOUR_PARAMETER = kwargs.get('CONTOUR_PARAMETER',
310310
np.arange(-100., 101, 1)).copy()
311+
self.INTERVAL = np.diff(self.CONTOUR_PARAMETER)[0]
311312
if 'Cyclonic' in SIGN_TYPE:
312313
self.CONTOUR_PARAMETER *= -1
313314

@@ -590,6 +591,7 @@ def create_netcdf(self, directory, savedir,
590591
else: # climatological ROMS solution
591592
nc.createVariable('ocean_time', 'f8', ('Nobs'), fill_value=self.FILLVAL)
592593

594+
nc.createVariable('cyc', np.int32, ('Nobs'), fill_value=self.FILLVAL)
593595
nc.createVariable('lon', 'f4', ('Nobs'), fill_value=self.FILLVAL)
594596
nc.createVariable('lat', 'f4', ('Nobs'), fill_value=self.FILLVAL)
595597
nc.createVariable('A', 'f4', ('Nobs'), fill_value=self.FILLVAL)
@@ -631,6 +633,14 @@ def create_netcdf(self, directory, savedir,
631633
else: # climatological ROMS solution
632634
nc.variables['ocean_time'].units = 'ROMS ocean_time (seconds)'
633635

636+
#nc.variables['eddy_duration'].units = 'days'
637+
nc.variables['cyc'].units = 'boolean'
638+
nc.variables['cyc'].min_val = -1
639+
nc.variables['cyc'].max_val = 1
640+
nc.variables['cyc'].long_name = 'cyclonic'
641+
nc.variables['cyc'].description = 'cyclonic -1; anti-cyclonic +1'
642+
643+
634644
#nc.variables['eddy_duration'].units = 'days'
635645
nc.variables['lon'].units = 'deg. longitude'
636646
nc.variables['lon'].min_val = self.LONMIN
@@ -777,13 +787,18 @@ def write2netcdf(self, rtime):
777787
radius_e = np.array([self.tracklist[i].radius_e]) * 1e-3 # to km
778788
n = np.arange(tsize, dtype=np.int32)
779789
track = np.full(tsize, self.ch_index)
790+
if 'Anticyclonic' in self.SIGN_TYPE:
791+
cyc = np.full(tsize, 1)
792+
elif 'Cyclonic' in self.SIGN_TYPE:
793+
cyc = np.full(tsize, -1)
780794
track_max_val = np.array([nc.variables['track'].max_val,
781795
np.int32(self.ch_index)]).max()
782796
#print self.tracklist[i].ocean_time
783797
#exit()
784798
#eddy_duration = np.array([self.tracklist[i].ocean_time]).ptp()
785799

786800
tend = self.ncind + tsize
801+
nc.variables['cyc'][self.ncind:tend] = cyc
787802
nc.variables['lon'][self.ncind:tend] = lon
788803
nc.variables['lat'][self.ncind:tend] = lat
789804
nc.variables['A'][self.ncind:tend] = amp

‎py_eddy_tracker_classes.py‎

Lines changed: 41 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -284,25 +284,25 @@ def psi2rho(var_psi):
284284

285285

286286

287-
def pcol_2dxy(x, y):
288-
"""
289-
Function to shift x, y for subsequent use with pcolor
290-
by Jeroen Molemaker UCLA 2008
291-
"""
292-
Mp, Lp = x.shape
293-
M = Mp - 1
294-
L = Lp - 1
295-
x_pcol = np.zeros((Mp, Lp))
296-
y_pcol = np.zeros_like(x_pcol)
297-
x_tmp = half_interp(x[:, :L], x[:, 1:Lp])
298-
x_pcol[1:Mp, 1:Lp] = half_interp(x_tmp[0:M], x_tmp[1:Mp])
299-
x_pcol[0] = 2. * x_pcol[1] - x_pcol[2]
300-
x_pcol[:, 0] = 2. * x_pcol[:, 1] - x_pcol[:, 2]
301-
y_tmp = half_interp(y[:, :L], y[:, 1:Lp] )
302-
y_pcol[1:Mp, 1:Lp] = half_interp(y_tmp[:M], y_tmp[1:Mp])
303-
y_pcol[0] = 2. * y_pcol[1] - y_pcol[2]
304-
y_pcol[:, 0] = 2. * y_pcol[:, 1] - y_pcol[:, 2]
305-
return x_pcol, y_pcol
287+
#def pcol_2dxy(x, y):
288+
#"""
289+
#Function to shift x, y for subsequent use with pcolor
290+
#by Jeroen Molemaker UCLA 2008
291+
#"""
292+
#Mp, Lp = x.shape
293+
#M = Mp - 1
294+
#L = Lp - 1
295+
#x_pcol = np.zeros((Mp, Lp))
296+
#y_pcol = np.zeros_like(x_pcol)
297+
#x_tmp = half_interp(x[:, :L], x[:, 1:Lp])
298+
#x_pcol[1:Mp, 1:Lp] = half_interp(x_tmp[0:M], x_tmp[1:Mp])
299+
#x_pcol[0] = 2. * x_pcol[1] - x_pcol[2]
300+
#x_pcol[:, 0] = 2. * x_pcol[:, 1] - x_pcol[:, 2]
301+
#y_tmp = half_interp(y[:, :L], y[:, 1:Lp] )
302+
#y_pcol[1:Mp, 1:Lp] = half_interp(y_tmp[:M], y_tmp[1:Mp])
303+
#y_pcol[0] = 2. * y_pcol[1] - y_pcol[2]
304+
#y_pcol[:, 0] = 2. * y_pcol[:, 1] - y_pcol[:, 2]
305+
#return x_pcol, y_pcol
306306

307307

308308

@@ -478,13 +478,15 @@ def get_uavg(Eddy, CS, collind, centlon_e, centlat_e, poly_eff,
478478

479479
any_inner_contours = True
480480

481-
seglon, seglat = poly_i.vertices[:, 0], poly_i.vertices[:, 1]
481+
seglon, seglat = (poly_i.vertices[:, 0],
482+
poly_i.vertices[:, 1])
482483
seglon, seglat = eddy_tracker.uniform_resample(
483484
seglon, seglat, method='akima')
484485

485486
# Interpolate uspd to seglon, seglat, then get mean
486487
if 'RectBivariate' in Eddy.INTERP_METHOD:
487-
uavgseg = Eddy.uspd_coeffs.ev(seglat[1:], seglon[1:]).mean()
488+
uavgseg = Eddy.uspd_coeffs.ev(seglat[1:],
489+
seglon[1:]).mean()
488490

489491
elif 'griddata' in Eddy.INTERP_METHOD:
490492
uavgseg = interpolate.griddata(points, uspd1d,
@@ -568,7 +570,7 @@ def collection_loop(CS, grd, rtime, A_list_obj, C_list_obj,
568570

569571
contlon_e, contlat_e = cont.vertices[:, 0].copy(), \
570572
cont.vertices[:, 1].copy()
571-
573+
572574
# Filter for closed contours
573575
if np.alltrue([contlon_e[0] == contlon_e[-1],
574576
contlat_e[0] == contlat_e[-1],
@@ -696,10 +698,24 @@ def collection_loop(CS, grd, rtime, A_list_obj, C_list_obj,
696698
if 'Anticyclonic' in sign_type:
697699
amp.all_pixels_above_h0()
698700

701+
reset_centroid = amp.all_pixels_above_h0(
702+
CS.levels[collind])
703+
#plt.figure(666)
704+
#plt.pcolormesh(Eddy.mask_eff)
705+
#plt.show()
706+
#amp.debug_figure(grd)
707+
699708
elif 'Cyclonic' in sign_type:
700-
amp.all_pixels_below_h0()
701-
702-
else: Exception
709+
reset_centroid = amp.all_pixels_below_h0(
710+
CS.levels[collind])
711+
else:
712+
Exception
713+
714+
if reset_centroid:
715+
centi = reset_centroid[0]
716+
centj = reset_centroid[1]
717+
centlon_e = grd.lon()[centj, centi]
718+
centlat_e = grd.lat()[centj, centi]
703719

704720
#amp.debug_figure(grd)
705721

0 commit comments

Comments
 (0)