Version 1.4.0

Added extra variables for use with ARGO

Author: evanmason

haversine_distmat.f90

@@ -27,7 +27,7 @@
 !
 !    If you have ifort on your system, change 'gfortran' to 'intelem'
 !
-!   Version 1.3.0
+!   Version 1.4.0
 !
 !===========================================================================
 

make_eddy_track_AVISO.py

@@ -24,7 +24,7 @@
 
 make_eddy_track_AVISO.py
 
-Version 1.3.0
+Version 1.4.0
 
 
 Scroll down to line ~640 to get started
@@ -785,7 +785,8 @@ def pcol_2dxy(self, x, y):
         qparameter = np.linspace(0, 5*10**-11, 25)
     elif 'sla' in diag_type:
         # Set SLA contour spacing
-        slaparameter = np.arange(-100., 101., 1.0) # cm
+        #slaparameter = np.arange(-100., 101., 1.0) # cm
+        slaparameter = np.arange(-100., 100.25, 0.25) # cm
 
 
 
@@ -886,9 +887,17 @@ def pcol_2dxy(self, x, y):
         ##sep_dist_fac = 1.5 # Causes tracks to jump for AVISO 7-day
 
 
+    #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
+
+
 
     cmap = plt.cm.RdBu
-    
     verbose = False
 
 
@@ -950,8 +959,8 @@ def pcol_2dxy(self, x, y):
     #axis = plt.axes()
 
     # Initialise two eddy objects to hold data
-    A_eddy = eddy_tracker.track_list('AVISO', track_duration_min)
-    C_eddy = eddy_tracker.track_list('AVISO', track_duration_min)
+    A_eddy = eddy_tracker.track_list('AVISO', track_duration_min, track_extra_variables)
+    C_eddy = eddy_tracker.track_list('AVISO', track_duration_min, track_extra_variables)
 
     if 'Q' in diag_type:
         A_savefile = "".join([savedir, 'eddy_tracks_Q_AVISO_anticyclonic.nc'])
@@ -1043,6 +1052,7 @@ def pcol_2dxy(self, x, y):
     A_eddy.verbose = verbose
     C_eddy.verbose = verbose
 
+    
     # See Chelton section B2 (0.4 degree radius)
     # These should give 8 and 1000 for 0.25 deg resolution
     pixmin = np.round((np.pi * radmin**2) / sla_grd.get_resolution()**2)
@@ -1369,12 +1379,12 @@ def pcol_2dxy(self, x, y):
             #saving_start_time = time.time()
             if not first_record:
                 if verbose:
-                    print '--- saving to nc', A_savefile
-                    print '--- saving to nc', C_savefile
+                    print '--- saving to nc', A_eddy.savedir
+                    print '--- saving to nc', C_eddy.savedir
                     print '+++'
                 if chelton_style_nc: # Recommended
-                    A_eddy.write2chelton_nc(A_savefile, rtime)
-                    C_eddy.write2chelton_nc(C_savefile, rtime)
+                    A_eddy.write2chelton_nc(rtime)
+                    C_eddy.write2chelton_nc(rtime)
                 else:
                     A_eddy.write2nc(A_savefile, rtime)
                     C_eddy.write2nc(C_savefile, rtime)

make_eddy_track_ROMS.py

@@ -24,7 +24,7 @@
 
 make_eddy_track_ROMS.py
 
-Version 1.3.0
+Version 1.4.0
 
 
 
@@ -90,7 +90,9 @@
     filetype     = 'avg'
     sigma_lev    = -1  # -1 for surface
 
-    days_btwn_recs = 3.
+    days_btwn_recs = 3. # days
+    # Save only tracks longer than...     
+    track_duration_min = 28 # days
 
     # Set grid
     grdname = 'roms_grd_NA2009_7pt5km.nc'
@@ -132,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., 1.0)
+        slaparameter = np.arange(-100., 101., .25)
 
 
     # Apply a filter to the Q parameter
@@ -149,8 +151,6 @@
     else:
         error
 
-    # Save only tracks longer than...     
-    track_duration_min = 5 # days
 
 
 
@@ -161,6 +161,11 @@
         lonmax = -5.5
         latmin = 20.5
         latmax = 38.5
+        
+        #lonmin = -25.     # Canary
+        #lonmax = -15
+        #latmin = 22.
+        #latmax = 28.
 
     # for canbas4
     #lonmin = -36.
@@ -186,10 +191,10 @@
 
     # Parameters used by CSS11 and KCCMC11 (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.05 #0.25 #0.15 # min change in amplitude
-    evolve_ammax = 5 #2.5 #5.  # max change in amplitude
-    evolve_armin = 0.05 #0.25 #0.15 # min change in area
-    evolve_armax = 5 #2.5 #5.  # max change in area
+    evolve_ammin = 0.0005 #0.25 #0.15 # min change in amplitude
+    evolve_ammax = 500 #2.5 #5.  # max change in amplitude
+    evolve_armin = 0.0005 #0.25 #0.15 # min change in area
+    evolve_armax = 500 #2.5 #5.  # max change in area
 
 
     separation_method = 'ellipse' # see CSS11
@@ -238,6 +243,7 @@
         #shape_err = 35. * np.ones(qparameter.size)
 
     elif 'sla' in diag_type:
+        #shape_err = 55. * np.ones(slaparameter.size)
         shape_err = 55. * np.ones(slaparameter.size)
         #shape_err = np.power(np.linspace(85., 40,  slaparameter.size), 2) / 100.
         #shape_err[shape_err < 50.] = 50.
@@ -387,7 +393,7 @@
     A_eddy.dist0 = np.float(dist0)
     C_eddy.dist0 = np.float(dist0)
 
-    start = 0
+    start = True
     start_time = time.time()
 
     # Loop thru the years...
@@ -402,11 +408,14 @@
             filename = directory + model + '_' + filetype + \
                        '_Y' + str(Yr) + 'M' + str(Mo) + '.nc'
             print 'Opening file:', filename
-            if start == 0:
-                start = 1
+            if start:
+                start = False
                 record_range = get_ROMS_data(filename)
-                A_eddy.days_btwn_recs = np.squeeze(np.diff(record_range[0:2]) / 86400.)
-                C_eddy.days_btwn_recs = np.squeeze(np.diff(record_range[0:2]) / 86400.)
+                
+                check_days_btwn_recs = np.squeeze(np.diff(record_range[0:2]) / 86400.)
+                assert days_btwn_recs == check_days_btwn_recs, '"days_btwn_recs" is incorrectly set'
+                A_eddy.days_btwn_recs = check_days_btwn_recs
+                C_eddy.days_btwn_recs = check_days_btwn_recs
 
                 if 'ip_roms' in model:
                     nc = netcdf.Dataset(filename)
@@ -595,12 +604,11 @@
                                        verbose=verbose)
 
                 elif 'sla' in diag_type:
-		    #print 'rtime 1', rtime
                     A_eddy.sign_type = 'Anticyclonic'
                     A_eddy = collection_loop(A_CS, grd, rtime,
                                              A_list_obj=A_eddy,  C_list_obj=None,
                                              sign_type=A_eddy.sign_type, verbose=verbose)
-                    # Note that CSc is reverse order
+                    # Note that CSc is in reverse order
                     C_eddy.sign_type = 'Cyclonic'
                     C_eddy = collection_loop(C_CS, grd, rtime,
                                              A_list_obj=None, C_list_obj=C_eddy,
@@ -640,11 +648,11 @@
                     plt.close(250)
 
                 if start:
-                    first_record = True
                     # Set old variables equal to new variables
                     A_eddy.set_old_variables()
                     C_eddy.set_old_variables()
-                    start = False
+                    first_record = True
+                    #start = False
                 else:
                     first_record = False