Version 1.1.0

Improved AVISO file iteration
Chosen start and end dates no longer have
to correspond exactly to AVISO file dates;
this sometimes causes problems with weekly files.

Author: evanmason

haversine_distmat.f90

@@ -24,7 +24,7 @@
 !    f2py -m haversine_distmat -h haversine_distmat.pyf haversine_distmat.f90  --overwrite-signature
 !    f2py -c --fcompiler=gfortran haversine_distmat.pyf haversine_distmat.f90
 !
-!   Version 1.0.1
+!   Version 1.1.0
 !
 !===========================================================================
 

make_eddy_track_AVISO.py

@@ -24,7 +24,7 @@
 
 make_eddy_track_AVISO.py
 
-Version 1.0.1
+Version 1.1.0
 
 
 Scroll down to line ~630 to get started
@@ -33,7 +33,9 @@
 
 
 import glob as glob
+#import matplotlib.pyplot as plt
 from py_eddy_tracker_classes import *
+from make_eddy_tracker_list_obj import *
 
 
 class py_eddy_tracker (object):
@@ -634,7 +636,6 @@ def pcol_2dxy(self, x, y):
     #diag_type = 'Q' <<< not implemented in 1.0.1
     diag_type = 'sla'
 
-    days_btwn_recs = 7.
 
     # Path to directory where outputs are to be saved...
     #savedir = directory
@@ -646,20 +647,16 @@ def pcol_2dxy(self, x, y):
 
     print '\nOutputs saved to', savedir
 
+    days_btwn_recs = 7.
 
-
-    
-    
-    
-    
-
 
 
     # Reference Julian day (Julian date at Jan 1, 1992)
     jday_ref = 2448623
 
     # Min and max permitted eddy radii [m]
     if 'Q' in diag_type:
+        Exception # Not implemented
         radmin = 15000.
         radmax = 250000.
         ampmin = 0.02
@@ -673,25 +670,22 @@ def pcol_2dxy(self, x, y):
 
 
     # Obtain this file from:
-    #  http://www-po.coas.oregonstate.edu/research/po/research/rossby_radius/
+    # http://www-po.coas.oregonstate.edu/research/po/research/rossby_radius/
     rw_path = '/home/emason/data_tmp/chelton_eddies/rossrad.dat'
-        
-    # Option to track using speed-based radius or effective radius
-    #track_with_speed = True
+
 
     # Make figures for animations
     anim_figs = True
 
 
     # Define contours
     if 'Q' in diag_type:
-        # Set Q contours
+        # Set Q contour spacing 
         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., 100.5, 0.5)
-        slaparameter = np.arange(-100., 101., 1.0)
-        #slaparameter = np.arange(-50., 51., 1.)
+        # Set SLA contour spacing
+        slaparameter = np.arange(-100., 101., 1.0) # cm
+
 
 
     # Apply a filter to the Q parameter
@@ -737,7 +731,7 @@ def pcol_2dxy(self, x, y):
 
 
     # Typical parameters
-    dist0 = 25000. # m separation distance after ~7 days (see Chelton2011 fig 22)
+    dist0 = 25000. # m separation distance after ~7 days (see CSS11 fig 22)
     if 'Q' in diag_type:
         amp0 = 0.02 # vort/f
     elif 'sla' in diag_type:
@@ -754,22 +748,20 @@ def pcol_2dxy(self, x, y):
     evolve_armax = 5  # max change in area
 
 
-    separation_method = 'ellipse' # see Chelton etal (2011)
+    separation_method = 'ellipse' # see CSS11
     #separation_method = 'sum_radii' # see Kurian etal (2011)
 
-    if 'sum_radii' in separation_method:
-        # Separation distance factor. Adjust according to number of days between records
-        # For 7 days, Chelton uses 150 km search ellipse
-        # So, given typical eddy radius of r=50 km, 1.5 * (r1 + r2) = 150 km.
-        #sep_dist_fac = 1.0
-        sep_dist_fac = 1.15 # Seems ok for AVISO 7-day
-        #sep_dist_fac = 1.5 # Causes tracks to jump for AVISO 7-day
+    #if 'sum_radii' in separation_method:
+        ## Separation distance factor. Adjust according to number of days between records
+        ## For 7 days, Chelton uses 150 km search ellipse
+        ## So, given typical eddy radius of r=50 km, 1.5 * (r1 + r2) = 150 km.
+        ##sep_dist_fac = 1.0
+        #sep_dist_fac = 1.15 # Seems ok for AVISO 7-day
+        ##sep_dist_fac = 1.5 # Causes tracks to jump for AVISO 7-day
 
 
-    #cmap = plt.cm.Spectral_r
-    #cmap = plt.cm.jet
-    cmap = plt.cm.RdBu
 
+    cmap = plt.cm.RdBu
 
     verbose = False
 
@@ -780,6 +772,9 @@ def pcol_2dxy(self, x, y):
     assert date_str < date_end, 'date_end must be larger than date_str'
     assert diag_type in ('Q','sla'), 'diag_type not properly defined'
 
+    thestartdate = dt.date2num(datestr2datetime(str(date_str)))
+    theenddate = dt.date2num(datestr2datetime(str(date_end)))
+    
     # Get complete AVISO file list
     AVISO_files = sorted(glob.glob(directory + AVISO_files))
 
@@ -793,10 +788,8 @@ def pcol_2dxy(self, x, y):
                               qparameter.size), 2)[::-1]
 
         # The shape error can maybe be played with...
-        #shape_err = np.power(np.linspace(200., 40,  qparameter.size), 2) / 100.
         shape_err = np.power(np.linspace(85., 40,  qparameter.size), 2) / 100.
         shape_err[shape_err < 35.] = 35.
-        #shape_err = 35. * np.ones(qparameter.size)
 
     elif 'sla' in diag_type:
         shape_err = 55. * np.ones(slaparameter.size)
@@ -808,7 +801,7 @@ def pcol_2dxy(self, x, y):
 
     Mx, My = sla_grd.Mx[sla_grd.jup0:sla_grd.jup1, sla_grd.iup0:sla_grd.iup1], \
              sla_grd.My[sla_grd.jup0:sla_grd.jup1, sla_grd.iup0:sla_grd.iup1]
-    pMx, pMy = sla_grd.pcol_2dxy(Mx, My)
+    #pMx, pMy = sla_grd.pcol_2dxy(Mx, My)
 
 
     # Get Rossby wave speed data
@@ -932,8 +925,6 @@ def pcol_2dxy(self, x, y):
     A_eddy.dist0 = np.float(dist0)
     C_eddy.dist0 = np.float(dist0)
 
-    start = 0
-    start_time = time.time()
     A_eddy.days_btwn_recs = days_btwn_recs
     C_eddy.days_btwn_recs = days_btwn_recs
 
@@ -943,15 +934,32 @@ def pcol_2dxy(self, x, y):
 
 
     # Loop through the AVISO files...
-    active = False
+    #active = False
+    start = True
+    start_time = time.time()
     print '\nStart tracking'
     for AVISO_file in AVISO_files:
-
-        if str(date_str) in AVISO_file:
+        
+        
+        with netcdf.Dataset(AVISO_file) as nc:
+	    
+            try:
+               thedate = nc.OriginalName
+               thedate = thedate.partition('qd_')[2].partition('_')[0]
+               thedate = datestr2datetime(thedate)
+               thedate = dt.date2num(thedate)
+            except Exception:
+               thedate = nc.variables['time']
+               raise Exception # 2014 AVISO not yet implemented
+        
+        
+        if np.logical_and(thedate >= thestartdate,
+                          thedate <= theenddate):
             active = True
-            file_time = time.time()
-            
-            
+        else:
+            active = False
+        
+        
         if active:
 
             #rec_start_time = time.time()
@@ -1161,13 +1169,14 @@ def pcol_2dxy(self, x, y):
                 plt.close(250)
 
 
-            if str(date_str) not in AVISO_file:
-                first_record = False
-            else:
+            if start:
                 first_record = True
                 # Set old variables equal to new variables
                 A_eddy.set_old_variables()
                 C_eddy.set_old_variables()
+                start = False
+            else:
+                first_record = False
 
 
             # Track the eddies

make_eddy_track_ROMS.py

@@ -24,7 +24,7 @@
 
 make_eddy_track_ROMS.py
 
-Version 1.0.1
+Version 1.1.0
 
 
 
@@ -281,7 +281,7 @@
 
     Mx, My = M(grd.lon()[jstr:jend,istr:iend],
                grd.lat()[jstr:jend,istr:iend])
-    pMx, pMy = pcol_2dxy(Mx, My)
+    #pMx, pMy = pcol_2dxy(Mx, My)
 
     # Get Rossby wave speed data
     rwv = eddy_tracker.rossby_wave_speeds(rw_path)

make_eddy_tracker_list_obj.py

@@ -25,7 +25,7 @@
 
 make_eddy_tracker_list_obj.py
 
-Version 1.0.1
+Version 1.1.0
 
 
 ===========================================================================
@@ -34,8 +34,8 @@
 """
 
 import haversine_distmat as hav # needs compiling with f2py
-from make_eddy_track import *
-
+#from make_eddy_track import *
+from py_eddy_tracker_classes import *
 
 def haversine_distance_vector(lon1, lat1, lon2, lat2):
     '''

py_eddy_tracker_classes.py

@@ -24,7 +24,7 @@
 
 py_eddy_tracker_classes.py
 
-Version 1.0.1
+Version 1.1.0
 ===========================================================================
 
 
@@ -52,6 +52,16 @@
 import haversine_distmat as hav # needs compiling with f2py
 
 
+def datestr2datetime(datestr):
+    '''
+    Take strings with format YYYYMMDD and convert to datetime instance
+    '''
+    assert len(datestr) == 8, "'datestr' must have length 8"
+    assert datestr[0] in ('1', '2'), "first character of 'datestr' should be '1' or '2'"
+    return dt.datetime.datetime(np.int(datestr[:4]),
+                                np.int(datestr[4:6]),
+                                np.int(datestr[6:]))
+
 
 def gaussian_resolution(res, zwl, mwl):
     '''