Skip to content

Commit 0b410dc

Browse files
evanmasonevanmason
committed
Bugfix: ROMS geostrophic velocity calculation
1 parent 96cfa93 commit 0b410dc

File tree

3 files changed

+89
-91
lines changed

3 files changed

+89
-91
lines changed

‎make_eddy_track_AVISO.py‎

Lines changed: 9 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -360,22 +360,16 @@ def make_gridmask(self, with_pad=True, use_maskoceans=False):
360360
self.Mx, self.My = x, y
361361
return self
362362

363-
def getSurfGeostrVel(self, zeta):
363+
def get_geostrophic_velocity(self, zeta):
364364
'''
365365
Returns u and v geostrophic velocity at
366366
surface from variables f, zeta, pm, pn...
367367
Note: output at rho points
368-
Adapted from IRD surf_geostr_vel.m function
369-
by Evan Mason
370-
Changed to gv2, 14 May 07...
371368
'''
372-
#def gv2(zeta, gof, pm, pn, umask, vmask): # Pierrick's version
373369
self.upad[:] = -self.gof() * self.v2rho_2d(self.vmask() * (zeta.data[1:] - zeta.data[:-1]) \
374370
* 0.5 * (self.pn()[1:] + self.pn()[:-1]))
375371
self.vpad[:] = self.gof() * self.u2rho_2d(self.umask() * (zeta.data[:, 1:] - zeta.data[:, :-1]) \
376372
* 0.5 * (self.pm()[:, 1:] + self.pm()[:, :-1]))
377-
#return ugv, vgv
378-
#return gv2(zeta, self.gof(), self.pm(), self.pn(), self.umask(), self.vmask())
379373
return self
380374

381375

@@ -735,12 +729,8 @@ def pcol_2dxy(self, x, y):
735729
# Path to directory where outputs are to be saved...
736730
#savedir = directory
737731
#savedir = '/marula/emason/aviso_eddy_tracking/pablo_exp/'
738-
savedir = '/marula/emason/aviso_eddy_tracking/new_AVISO_test/'
732+
#savedir = '/marula/emason/aviso_eddy_tracking/new_AVISO_test/'
739733
savedir = '/marula/emason/aviso_eddy_tracking/junk/'
740-
#savedir = '/marula/emason/aviso_eddy_tracking/new_AVISO_test_0.35/'
741-
#savedir = '/marula/emason/aviso_eddy_tracking/new_AVISO_test_0.3/'
742-
#savedir = '/marula/emason/aviso_eddy_tracking/new_AVISO/CEC/'
743-
#savedir = '/shared/emason/eddy_tracks/Barbara/2009/'
744734
#savedir = '/marula/emason/aviso_eddy_tracking/new_AVISO_test/BlackSea/'
745735
#savedir = '/path/to/save/your/outputs/'
746736

@@ -812,9 +802,9 @@ def pcol_2dxy(self, x, y):
812802
if the_domain in 'Global':
813803

814804
lonmin = -36. # Canary
815-
lonmax = -7.5
805+
lonmax = -5.5
816806
latmin = 18.
817-
latmax = 33.2
807+
latmax = 35.5
818808

819809
#lonmin = -30. # BENCS
820810
#lonmax = 22.
@@ -869,10 +859,10 @@ def pcol_2dxy(self, x, y):
869859

870860
# Parameters used by Chelton etal and Kurian etal (Sec. 3.2) to ensure the slow evolution
871861
# of the eddies over time; they use min and max values of 0.25 and 2.5
872-
evolve_ammin = 0.01# 0.25 # min change in amplitude
873-
evolve_ammax = 10#2.5 # max change in amplitude
874-
evolve_armin = 0.01# 0.25 # min change in area
875-
evolve_armax = 10 # max change in area
862+
evolve_ammin = 0.0005# 0.25 # min change in amplitude
863+
evolve_ammax = 500#2.5 # max change in amplitude
864+
evolve_armin = 0.0005# 0.25 # min change in area
865+
evolve_armax = 500 # max change in area
876866

877867

878868
separation_method = 'ellipse' # see CSS11
@@ -1212,7 +1202,7 @@ def pcol_2dxy(self, x, y):
12121202
#umask, vmask, junk = sla_grd.uvpmask(np.zeros_like(sla))
12131203

12141204
# Multiply by 0.01 for m
1215-
sla_grd.getSurfGeostrVel(sla * 0.01)
1205+
sla_grd.get_geostrophic_velocity(sla * 0.01)
12161206

12171207
# Remove padded boundary
12181208
sla = sla[sla_grd.jup0:sla_grd.jup1, sla_grd.iup0:sla_grd.iup1]

‎make_eddy_track_ROMS.py‎

Lines changed: 17 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -134,7 +134,7 @@
134134
qparameter = np.linspace(0, 5*10**-11, 25)
135135

136136
elif 'sla' in diag_type: # Units, cm; for anticyclones use -50 to 50
137-
slaparameter = np.arange(-100., 101., .25)
137+
slaparameter = np.arange(-100., 100.25, 0.25)
138138

139139

140140
# Apply a filter to the Q parameter
@@ -157,10 +157,10 @@
157157
subdomain = True
158158
if the_domain in 'Canary_Islands':
159159

160-
lonmin = -38. # Canary
160+
lonmin = -36. # Canary
161161
lonmax = -5.5
162-
latmin = 20.5
163-
latmax = 38.5
162+
latmin = 18.
163+
latmax = 35.5
164164

165165
#lonmin = -25. # Canary
166166
#lonmax = -15
@@ -215,6 +215,14 @@
215215

216216
verbose = False
217217

218+
#Define track_extra_variables to track and save:
219+
# - effective contour points
220+
# - speed-based contour points
221+
# - shape test values
222+
# - profiles of swirl velocity from effective contour inwards
223+
# Useful for working with ARGO data
224+
track_extra_variables = True
225+
218226

219227
# Note this is a global value
220228
fillval = -9999
@@ -292,8 +300,8 @@
292300

293301
## Initialise
294302
# Objects to hold data (qparameter)
295-
A_eddy = eddy_tracker.track_list('ROMS', track_duration_min)
296-
C_eddy = eddy_tracker.track_list('ROMS', track_duration_min)
303+
A_eddy = eddy_tracker.track_list('ROMS', track_duration_min, track_extra_variables)
304+
C_eddy = eddy_tracker.track_list('ROMS', track_duration_min, track_extra_variables)
297305

298306
if 'Q' in diag_type:
299307
A_savefile = "".join([savedir, 'eddy_tracks_anticyclonic.nc'])
@@ -526,7 +534,7 @@
526534
#sla *= grd.mask()
527535
sla = np.ma.masked_where(grd.mask() == 0, sla)
528536

529-
grd.getSurfGeostrVel(sla)#, grd.f(), grd.pm(), grd.pn(),
537+
grd.get_geostrophic_velocity(sla)#, grd.f(), grd.pm(), grd.pn(),
530538
#grd.u_mask, grd.v_mask)
531539

532540
sla *= 100. # m to cm
@@ -697,8 +705,8 @@
697705
if verbose:
698706
print '--- saving to nc'
699707
#print 'ddsssssss', rtime
700-
A_eddy.write2chelton_nc(A_savefile, rtime)
701-
C_eddy.write2chelton_nc(C_savefile, rtime)
708+
A_eddy.write2chelton_nc(rtime)
709+
C_eddy.write2chelton_nc(rtime)
702710
#print 'Saving the eddies', time.time() - saving_start_time, 'seconds'
703711

704712
# Running time for a single monthly file

‎py_eddy_tracker_classes.py‎

Lines changed: 63 additions & 63 deletions
Original file line numberDiff line numberDiff line change
@@ -449,11 +449,11 @@ def calc_uavg(points, uspd, lon, lat):
449449
return np.mean(uavg[np.isfinite(uavg)])
450450

451451
# True for debug figures
452-
debug_U = False
453-
if debug_U:
454-
schem_dic = {}
455-
conts_x = np.array([])
456-
conts_y = np.array([])
452+
#debug_U = False
453+
#if debug_U:
454+
#schem_dic = {}
455+
#conts_x = np.array([])
456+
#conts_y = np.array([])
457457

458458
# Unpack indices for convenience
459459
istr, iend, jstr, jend = Eddy.i0, Eddy.i1, Eddy.j0, Eddy.j1
@@ -520,42 +520,42 @@ def calc_uavg(points, uspd, lon, lat):
520520
seglon, seglat = poly_i.vertices[:,0], poly_i.vertices[:,1]
521521
seglon, seglat = eddy_tracker.uniform_resample(seglon, seglat)
522522

523-
if debug_U:
524-
px, py = pcol_2dxy(grd.lon()[jmin:jmax,imin:imax],
525-
grd.lat()[jmin:jmax,imin:imax])
526-
plt.figure(55)
527-
ax1 = plt.subplot(121)
528-
if start:
529-
pcm = ax1.pcolormesh(px, py, Eddy.Uspd[jmin:jmax,imin:imax], cmap=plt.cm.gist_earth_r)
530-
pcm.set_clim(0, .5)
531-
plt.colorbar(pcm, orientation='horizontal')
532-
plt.scatter(grd.lon()[jmin:jmax,imin:imax],
533-
grd.lat()[jmin:jmax,imin:imax], s=5, c='k')
534-
start = False
535-
plt.plot(seglon, seglat, '.-r')
536-
plt.plot(Eddy.circlon, Eddy.circlat, 'b')
537-
plt.scatter(centlon_e, centlat_e, s=125, c='r')
538-
plt.axis('image')
539-
ax2 = plt.subplot(122)
540-
ax2.set_title('sum mask: %s' %np.sum(mask_i))
541-
plt.pcolormesh(px, py, mask_i.reshape(px.shape), cmap=plt.cm.Accent, edgecolors='orange')
542-
plt.scatter(grd.lon()[jmin:jmax,imin:imax],
543-
grd.lat()[jmin:jmax,imin:imax], s=5, c='k')
544-
plt.plot(seglon, seglat, '.-r')
545-
conts_x = np.append(conts_x, seglon)
546-
conts_y = np.append(conts_y, seglat)
547-
conts_x = np.append(conts_x, 9999)
548-
conts_y = np.append(conts_y, 9999)
549-
plt.plot(Eddy.circlon, Eddy.circlat, 'b')
550-
plt.scatter(centlon_e, centlat_e, s=125, c='r')
523+
#if debug_U:
524+
#px, py = pcol_2dxy(grd.lon()[jmin:jmax,imin:imax],
525+
#grd.lat()[jmin:jmax,imin:imax])
526+
#plt.figure(55)
527+
#ax1 = plt.subplot(121)
528+
#if start:
529+
#pcm = ax1.pcolormesh(px, py, Eddy.Uspd[jmin:jmax,imin:imax], cmap=plt.cm.gist_earth_r)
530+
#pcm.set_clim(0, .5)
531+
#plt.colorbar(pcm, orientation='horizontal')
532+
#plt.scatter(grd.lon()[jmin:jmax,imin:imax],
533+
#grd.lat()[jmin:jmax,imin:imax], s=5, c='k')
534+
#start = False
535+
#plt.plot(seglon, seglat, '.-r')
536+
#plt.plot(Eddy.circlon, Eddy.circlat, 'b')
537+
#plt.scatter(centlon_e, centlat_e, s=125, c='r')
538+
#plt.axis('image')
539+
#ax2 = plt.subplot(122)
540+
#ax2.set_title('sum mask: %s' %np.sum(mask_i))
541+
#plt.pcolormesh(px, py, mask_i.reshape(px.shape), cmap=plt.cm.Accent, edgecolors='orange')
542+
#plt.scatter(grd.lon()[jmin:jmax,imin:imax],
543+
#grd.lat()[jmin:jmax,imin:imax], s=5, c='k')
544+
#plt.plot(seglon, seglat, '.-r')
545+
#conts_x = np.append(conts_x, seglon)
546+
#conts_y = np.append(conts_y, seglat)
547+
#conts_x = np.append(conts_x, 9999)
548+
#conts_y = np.append(conts_y, 9999)
549+
#plt.plot(Eddy.circlon, Eddy.circlat, 'b')
550+
#plt.scatter(centlon_e, centlat_e, s=125, c='r')
551551

552552
# Interpolate Uspd to seglon, seglat, then get mean
553553
Uavgseg = calc_uavg(points, Eddy.Uspd[jmin:jmax,imin:imax], seglon[:-1], seglat[:-1])
554554

555555
if save_all_uavg:
556556
all_uavg.append(Uavgseg)
557557

558-
if Uavgseg >= Uavg:
558+
if Uavgseg >= Uavg and np.sum(mask_i) >= Eddy.pixel_threshold[0]:
559559
Uavg = Uavgseg.copy()
560560
theseglon, theseglat = seglon.copy(), seglat.copy()
561561

@@ -569,42 +569,42 @@ def calc_uavg(points, uspd, lon, lat):
569569
# Speed based eddy radius (eddy_radius_s)
570570
centx_s, centy_s, eddy_radius_s, junk, junk = fit_circle(cx, cy)
571571
centlon_s, centlat_s = Eddy.M.projtran(centx_s, centy_s, inverse=True)
572-
if debug_U:
573-
ax1.set_title('Speed-based radius: %s km' %np.int(eddy_radius_s/1e3))
574-
ax1.plot(poly_e.vertices[:,0], poly_e.vertices[:,1], 'om')
575-
ax1.plot(theseglon, theseglat, 'g', lw=3)
576-
ax1.plot(inner_seglon, inner_seglat, 'k')
577-
ax1.scatter(centlon_s, centlat_s, s=50, c='g')
578-
ax1.axis('image')
579-
ax2.plot(theseglon, theseglat, 'g', lw=3)
580-
ax1.plot(inner_seglon, inner_seglat, 'k')
581-
ax2.axis('image')
582-
plt.show()
583-
schem_dic['lon'] = grd.lon()[jmin:jmax,imin:imax]
584-
schem_dic['lat'] = grd.lat()[jmin:jmax,imin:imax]
585-
schem_dic['spd'] = Eddy.Uspd[jmin:jmax,imin:imax]
586-
schem_dic['inner_seg'] = np.array([inner_seglon, inner_seglat]).T
587-
schem_dic['effective_seg'] = np.array([poly_e.vertices[:,0], poly_e.vertices[:,1]]).T
588-
schem_dic['effective_circ'] = np.array([Eddy.circlon, Eddy.circlat]).T
589-
schem_dic['speed_seg'] = np.array([theseglon, theseglat]).T
590-
schem_dic['all_segs'] = np.array([conts_x, conts_y]).T
591-
schem_dic['Ls_centroid'] = np.array([centlon_s, centlat_s])
592-
schem_dic['Leff_centroid'] = np.array([centlon_e, centlat_e])
593-
io.savemat('schematic_fig_%s' %np.round(eddy_radius_s/1e3, 3), schem_dic)
572+
#if debug_U:
573+
#ax1.set_title('Speed-based radius: %s km' %np.int(eddy_radius_s/1e3))
574+
#ax1.plot(poly_e.vertices[:,0], poly_e.vertices[:,1], 'om')
575+
#ax1.plot(theseglon, theseglat, 'g', lw=3)
576+
#ax1.plot(inner_seglon, inner_seglat, 'k')
577+
#ax1.scatter(centlon_s, centlat_s, s=50, c='g')
578+
#ax1.axis('image')
579+
#ax2.plot(theseglon, theseglat, 'g', lw=3)
580+
#ax1.plot(inner_seglon, inner_seglat, 'k')
581+
#ax2.axis('image')
582+
#plt.show()
583+
#schem_dic['lon'] = grd.lon()[jmin:jmax,imin:imax]
584+
#schem_dic['lat'] = grd.lat()[jmin:jmax,imin:imax]
585+
#schem_dic['spd'] = Eddy.Uspd[jmin:jmax,imin:imax]
586+
#schem_dic['inner_seg'] = np.array([inner_seglon, inner_seglat]).T
587+
#schem_dic['effective_seg'] = np.array([poly_e.vertices[:,0], poly_e.vertices[:,1]]).T
588+
#schem_dic['effective_circ'] = np.array([Eddy.circlon, Eddy.circlat]).T
589+
#schem_dic['speed_seg'] = np.array([theseglon, theseglat]).T
590+
#schem_dic['all_segs'] = np.array([conts_x, conts_y]).T
591+
#schem_dic['Ls_centroid'] = np.array([centlon_s, centlat_s])
592+
#schem_dic['Leff_centroid'] = np.array([centlon_e, centlat_e])
593+
#io.savemat('schematic_fig_%s' %np.round(eddy_radius_s/1e3, 3), schem_dic)
594594
if not save_all_uavg:
595595
return Uavg, centlon_s, centlat_s, eddy_radius_s, theseglon, theseglat, inner_seglon, inner_seglat
596596
else:
597597
return Uavg, centlon_s, centlat_s, eddy_radius_s, theseglon, theseglat, inner_seglon, inner_seglat, all_uavg
598598

599599
except Exception: # If no interior contours found, use eddy_radius_e
600600

601-
if debug_U and 0:
602-
plt.title('Effective radius: %s km' %np.int(eddy_radius_e/1e3))
603-
plt.plot(poly_e.vertices[:,0], poly_e.vertices[:,1], 'om')
604-
plt.plot(theseglon, theseglat, 'g')
605-
plt.scatter(centlon_e, centlat_e, s=125, c='r')
606-
plt.axis('image')
607-
plt.show()
601+
#if debug_U and 0:
602+
#plt.title('Effective radius: %s km' %np.int(eddy_radius_e/1e3))
603+
#plt.plot(poly_e.vertices[:,0], poly_e.vertices[:,1], 'om')
604+
#plt.plot(theseglon, theseglat, 'g')
605+
#plt.scatter(centlon_e, centlat_e, s=125, c='r')
606+
#plt.axis('image')
607+
#plt.show()
608608
if not save_all_uavg:
609609
return Uavg, centlon_e, centlat_e, eddy_radius_e, theseglon, theseglat, theseglon, theseglat
610610
else:

0 commit comments

Comments
 (0)