Bugfix: ROMS geostrophic velocity calculation

Author: evanmason

make_eddy_track_AVISO.py

@@ -360,22 +360,16 @@ def make_gridmask(self, with_pad=True, use_maskoceans=False):
         self.Mx, self.My = x, y
         return self
 
-    def getSurfGeostrVel(self, zeta):
+    def get_geostrophic_velocity(self, zeta):
         '''
         Returns u and v geostrophic velocity at
         surface from variables f, zeta, pm, pn...
         Note: output at rho points
-        Adapted from IRD surf_geostr_vel.m function
-        by Evan Mason
-        Changed to gv2, 14 May 07...
         '''
-        #def gv2(zeta, gof, pm, pn, umask, vmask): # Pierrick's version
         self.upad[:] = -self.gof() * self.v2rho_2d(self.vmask() * (zeta.data[1:] - zeta.data[:-1]) \
                                         * 0.5 * (self.pn()[1:] + self.pn()[:-1]))
         self.vpad[:] =  self.gof() * self.u2rho_2d(self.umask() * (zeta.data[:, 1:] - zeta.data[:, :-1]) \
                                         * 0.5 * (self.pm()[:, 1:] + self.pm()[:, :-1]))
-            #return ugv, vgv
-        #return gv2(zeta, self.gof(), self.pm(), self.pn(), self.umask(), self.vmask())
         return self
 
 
@@ -735,12 +729,8 @@ def pcol_2dxy(self, x, y):
     # Path to directory where outputs are to be saved...
     #savedir = directory
     #savedir = '/marula/emason/aviso_eddy_tracking/pablo_exp/'
-    savedir = '/marula/emason/aviso_eddy_tracking/new_AVISO_test/'
+    #savedir = '/marula/emason/aviso_eddy_tracking/new_AVISO_test/'
     savedir = '/marula/emason/aviso_eddy_tracking/junk/'
-    #savedir = '/marula/emason/aviso_eddy_tracking/new_AVISO_test_0.35/'
-    #savedir = '/marula/emason/aviso_eddy_tracking/new_AVISO_test_0.3/'
-    #savedir = '/marula/emason/aviso_eddy_tracking/new_AVISO/CEC/'
-    #savedir = '/shared/emason/eddy_tracks/Barbara/2009/'
     #savedir = '/marula/emason/aviso_eddy_tracking/new_AVISO_test/BlackSea/'
     #savedir = '/path/to/save/your/outputs/'
 
@@ -812,9 +802,9 @@ def pcol_2dxy(self, x, y):
     if the_domain in 'Global':
 
         lonmin = -36.     # Canary
-        lonmax = -7.5
+        lonmax = -5.5
         latmin = 18.
-        latmax = 33.2
+        latmax = 35.5
 
         #lonmin = -30.        # BENCS
         #lonmax =  22.
@@ -869,10 +859,10 @@ def pcol_2dxy(self, x, y):
 
     # Parameters used by Chelton etal and Kurian etal (Sec. 3.2) to ensure the slow evolution
     # of the eddies over time; they use min and max values of 0.25 and 2.5
-    evolve_ammin = 0.01# 0.25 # min change in amplitude
-    evolve_ammax = 10#2.5  # max change in amplitude
-    evolve_armin = 0.01# 0.25 # min change in area
-    evolve_armax = 10  # max change in area
+    evolve_ammin = 0.0005# 0.25 # min change in amplitude
+    evolve_ammax = 500#2.5  # max change in amplitude
+    evolve_armin = 0.0005# 0.25 # min change in area
+    evolve_armax = 500  # max change in area
 
 
     separation_method = 'ellipse' # see CSS11
@@ -1212,7 +1202,7 @@ def pcol_2dxy(self, x, y):
                     #umask, vmask, junk = sla_grd.uvpmask(np.zeros_like(sla))
 
                 # Multiply by 0.01 for m
-                sla_grd.getSurfGeostrVel(sla * 0.01)
+                sla_grd.get_geostrophic_velocity(sla * 0.01)
 
                 # Remove padded boundary
                 sla = sla[sla_grd.jup0:sla_grd.jup1, sla_grd.iup0:sla_grd.iup1]

make_eddy_track_ROMS.py

@@ -134,7 +134,7 @@
         qparameter = np.linspace(0, 5*10**-11, 25)
 
     elif 'sla' in diag_type: # Units, cm; for anticyclones use -50 to 50
-        slaparameter = np.arange(-100., 101., .25)
+        slaparameter = np.arange(-100., 100.25, 0.25)
 
 
     # Apply a filter to the Q parameter
@@ -157,10 +157,10 @@
     subdomain = True
     if the_domain in 'Canary_Islands':
 
-        lonmin = -38.     # Canary
+        lonmin = -36.     # Canary
         lonmax = -5.5
-        latmin = 20.5
-        latmax = 38.5
+        latmin = 18.
+        latmax = 35.5
 
         #lonmin = -25.     # Canary
         #lonmax = -15
@@ -215,6 +215,14 @@
 
     verbose = False
 
+    #Define track_extra_variables to track and save:
+    # - effective contour points
+    # - speed-based contour points
+    # - shape test values
+    # - profiles of swirl velocity from effective contour inwards
+    # Useful for working with ARGO data
+    track_extra_variables = True
+    
 
     # Note this is a global value
     fillval = -9999
@@ -292,8 +300,8 @@
 
     ## Initialise
     # Objects to hold data (qparameter)
-    A_eddy = eddy_tracker.track_list('ROMS', track_duration_min)
-    C_eddy = eddy_tracker.track_list('ROMS', track_duration_min)
+    A_eddy = eddy_tracker.track_list('ROMS', track_duration_min, track_extra_variables)
+    C_eddy = eddy_tracker.track_list('ROMS', track_duration_min, track_extra_variables)
 
     if 'Q' in diag_type:
         A_savefile = "".join([savedir, 'eddy_tracks_anticyclonic.nc'])
@@ -526,7 +534,7 @@
                     #sla *= grd.mask()
                     sla = np.ma.masked_where(grd.mask() == 0, sla)
 
-                    grd.getSurfGeostrVel(sla)#, grd.f(), grd.pm(), grd.pn(),
+                    grd.get_geostrophic_velocity(sla)#, grd.f(), grd.pm(), grd.pn(),
                                                 #grd.u_mask, grd.v_mask)
 
                     sla *= 100. # m to cm
@@ -697,8 +705,8 @@
                     if verbose:
                         print '--- saving to nc'
                     #print 'ddsssssss', rtime
-                    A_eddy.write2chelton_nc(A_savefile, rtime)
-                    C_eddy.write2chelton_nc(C_savefile, rtime)
+                    A_eddy.write2chelton_nc(rtime)
+                    C_eddy.write2chelton_nc(rtime)
                 #print 'Saving the eddies', time.time() - saving_start_time, 'seconds'
 
             # Running time for a single monthly file

py_eddy_tracker_classes.py

@@ -449,11 +449,11 @@ def calc_uavg(points, uspd, lon, lat):
         return np.mean(uavg[np.isfinite(uavg)])
 
     # True for debug figures
-    debug_U = False
-    if debug_U:
-        schem_dic = {}
-        conts_x = np.array([])
-        conts_y = np.array([])
+    #debug_U = False
+    #if debug_U:
+        #schem_dic = {}
+        #conts_x = np.array([])
+        #conts_y = np.array([])
 
     # Unpack indices for convenience
     istr, iend, jstr, jend = Eddy.i0, Eddy.i1, Eddy.j0, Eddy.j1
@@ -520,42 +520,42 @@ def calc_uavg(points, uspd, lon, lat):
                 seglon, seglat = poly_i.vertices[:,0], poly_i.vertices[:,1]
                 seglon, seglat = eddy_tracker.uniform_resample(seglon, seglat)
 
-                if debug_U:
-                    px, py = pcol_2dxy(grd.lon()[jmin:jmax,imin:imax],
-                                       grd.lat()[jmin:jmax,imin:imax])
-                    plt.figure(55)
-                    ax1 = plt.subplot(121)
-                    if start:
-                        pcm = ax1.pcolormesh(px, py, Eddy.Uspd[jmin:jmax,imin:imax], cmap=plt.cm.gist_earth_r)
-                        pcm.set_clim(0, .5)
-                        plt.colorbar(pcm, orientation='horizontal')
-                        plt.scatter(grd.lon()[jmin:jmax,imin:imax],
-                                    grd.lat()[jmin:jmax,imin:imax], s=5, c='k')
-                        start = False
-                    plt.plot(seglon, seglat, '.-r')
-                    plt.plot(Eddy.circlon, Eddy.circlat, 'b')
-                    plt.scatter(centlon_e, centlat_e, s=125, c='r')
-                    plt.axis('image')
-                    ax2 = plt.subplot(122)
-                    ax2.set_title('sum mask: %s' %np.sum(mask_i))
-                    plt.pcolormesh(px, py, mask_i.reshape(px.shape), cmap=plt.cm.Accent, edgecolors='orange')
-                    plt.scatter(grd.lon()[jmin:jmax,imin:imax],
-                                grd.lat()[jmin:jmax,imin:imax], s=5, c='k')
-                    plt.plot(seglon, seglat, '.-r')
-                    conts_x = np.append(conts_x, seglon)
-                    conts_y = np.append(conts_y, seglat)
-                    conts_x = np.append(conts_x, 9999)
-                    conts_y = np.append(conts_y, 9999)
-                    plt.plot(Eddy.circlon, Eddy.circlat, 'b')
-                    plt.scatter(centlon_e, centlat_e, s=125, c='r')
+                #if debug_U:
+                    #px, py = pcol_2dxy(grd.lon()[jmin:jmax,imin:imax],
+                                       #grd.lat()[jmin:jmax,imin:imax])
+                    #plt.figure(55)
+                    #ax1 = plt.subplot(121)
+                    #if start:
+                        #pcm = ax1.pcolormesh(px, py, Eddy.Uspd[jmin:jmax,imin:imax], cmap=plt.cm.gist_earth_r)
+                        #pcm.set_clim(0, .5)
+                        #plt.colorbar(pcm, orientation='horizontal')
+                        #plt.scatter(grd.lon()[jmin:jmax,imin:imax],
+                                    #grd.lat()[jmin:jmax,imin:imax], s=5, c='k')
+                        #start = False
+                    #plt.plot(seglon, seglat, '.-r')
+                    #plt.plot(Eddy.circlon, Eddy.circlat, 'b')
+                    #plt.scatter(centlon_e, centlat_e, s=125, c='r')
+                    #plt.axis('image')
+                    #ax2 = plt.subplot(122)
+                    #ax2.set_title('sum mask: %s' %np.sum(mask_i))
+                    #plt.pcolormesh(px, py, mask_i.reshape(px.shape), cmap=plt.cm.Accent, edgecolors='orange')
+                    #plt.scatter(grd.lon()[jmin:jmax,imin:imax],
+                                #grd.lat()[jmin:jmax,imin:imax], s=5, c='k')
+                    #plt.plot(seglon, seglat, '.-r')
+                    #conts_x = np.append(conts_x, seglon)
+                    #conts_y = np.append(conts_y, seglat)
+                    #conts_x = np.append(conts_x, 9999)
+                    #conts_y = np.append(conts_y, 9999)
+                    #plt.plot(Eddy.circlon, Eddy.circlat, 'b')
+                    #plt.scatter(centlon_e, centlat_e, s=125, c='r')
 
                 # Interpolate Uspd to seglon, seglat, then get mean
                 Uavgseg = calc_uavg(points, Eddy.Uspd[jmin:jmax,imin:imax], seglon[:-1], seglat[:-1])
 
                 if save_all_uavg:
                     all_uavg.append(Uavgseg)
 
-                if Uavgseg >= Uavg:
+                if Uavgseg >= Uavg and np.sum(mask_i) >= Eddy.pixel_threshold[0]:
                     Uavg = Uavgseg.copy()
                     theseglon, theseglat = seglon.copy(), seglat.copy()
 
@@ -569,42 +569,42 @@ def calc_uavg(points, uspd, lon, lat):
         # Speed based eddy radius (eddy_radius_s)
         centx_s, centy_s, eddy_radius_s, junk, junk = fit_circle(cx, cy)
         centlon_s, centlat_s = Eddy.M.projtran(centx_s, centy_s, inverse=True)
-        if debug_U:
-            ax1.set_title('Speed-based radius: %s km' %np.int(eddy_radius_s/1e3))
-            ax1.plot(poly_e.vertices[:,0], poly_e.vertices[:,1], 'om')
-            ax1.plot(theseglon, theseglat, 'g', lw=3)
-            ax1.plot(inner_seglon, inner_seglat, 'k')
-            ax1.scatter(centlon_s, centlat_s, s=50, c='g')
-            ax1.axis('image')
-            ax2.plot(theseglon, theseglat, 'g', lw=3)
-            ax1.plot(inner_seglon, inner_seglat, 'k')
-            ax2.axis('image')
-            plt.show()
-            schem_dic['lon'] = grd.lon()[jmin:jmax,imin:imax]
-            schem_dic['lat'] = grd.lat()[jmin:jmax,imin:imax]
-            schem_dic['spd'] = Eddy.Uspd[jmin:jmax,imin:imax]
-            schem_dic['inner_seg'] = np.array([inner_seglon, inner_seglat]).T
-            schem_dic['effective_seg'] = np.array([poly_e.vertices[:,0], poly_e.vertices[:,1]]).T
-            schem_dic['effective_circ'] = np.array([Eddy.circlon, Eddy.circlat]).T
-            schem_dic['speed_seg'] = np.array([theseglon, theseglat]).T
-            schem_dic['all_segs'] = np.array([conts_x, conts_y]).T
-            schem_dic['Ls_centroid'] = np.array([centlon_s, centlat_s])
-            schem_dic['Leff_centroid'] = np.array([centlon_e, centlat_e])
-            io.savemat('schematic_fig_%s' %np.round(eddy_radius_s/1e3, 3), schem_dic)
+        #if debug_U:
+            #ax1.set_title('Speed-based radius: %s km' %np.int(eddy_radius_s/1e3))
+            #ax1.plot(poly_e.vertices[:,0], poly_e.vertices[:,1], 'om')
+            #ax1.plot(theseglon, theseglat, 'g', lw=3)
+            #ax1.plot(inner_seglon, inner_seglat, 'k')
+            #ax1.scatter(centlon_s, centlat_s, s=50, c='g')
+            #ax1.axis('image')
+            #ax2.plot(theseglon, theseglat, 'g', lw=3)
+            #ax1.plot(inner_seglon, inner_seglat, 'k')
+            #ax2.axis('image')
+            #plt.show()
+            #schem_dic['lon'] = grd.lon()[jmin:jmax,imin:imax]
+            #schem_dic['lat'] = grd.lat()[jmin:jmax,imin:imax]
+            #schem_dic['spd'] = Eddy.Uspd[jmin:jmax,imin:imax]
+            #schem_dic['inner_seg'] = np.array([inner_seglon, inner_seglat]).T
+            #schem_dic['effective_seg'] = np.array([poly_e.vertices[:,0], poly_e.vertices[:,1]]).T
+            #schem_dic['effective_circ'] = np.array([Eddy.circlon, Eddy.circlat]).T
+            #schem_dic['speed_seg'] = np.array([theseglon, theseglat]).T
+            #schem_dic['all_segs'] = np.array([conts_x, conts_y]).T
+            #schem_dic['Ls_centroid'] = np.array([centlon_s, centlat_s])
+            #schem_dic['Leff_centroid'] = np.array([centlon_e, centlat_e])
+            #io.savemat('schematic_fig_%s' %np.round(eddy_radius_s/1e3, 3), schem_dic)
         if not save_all_uavg:
             return Uavg, centlon_s, centlat_s, eddy_radius_s, theseglon, theseglat, inner_seglon, inner_seglat
         else:  
             return Uavg, centlon_s, centlat_s, eddy_radius_s, theseglon, theseglat, inner_seglon, inner_seglat, all_uavg
 
     except Exception: # If no interior contours found, use eddy_radius_e
 
-        if debug_U and 0:
-            plt.title('Effective radius: %s km' %np.int(eddy_radius_e/1e3))
-            plt.plot(poly_e.vertices[:,0], poly_e.vertices[:,1], 'om')
-            plt.plot(theseglon, theseglat, 'g')
-            plt.scatter(centlon_e, centlat_e, s=125, c='r')
-            plt.axis('image')
-            plt.show()
+        #if debug_U and 0:
+            #plt.title('Effective radius: %s km' %np.int(eddy_radius_e/1e3))
+            #plt.plot(poly_e.vertices[:,0], poly_e.vertices[:,1], 'om')
+            #plt.plot(theseglon, theseglat, 'g')
+            #plt.scatter(centlon_e, centlat_e, s=125, c='r')
+            #plt.axis('image')
+            #plt.show()
         if not save_all_uavg:
             return Uavg, centlon_e, centlat_e, eddy_radius_e, theseglon, theseglat, theseglon, theseglat
         else: