make_eddy_track_CLS.py fixed

evanmason · evanmason · commit 7ae56151dcfa · 2015-03-15T13:47:29.000+01:00
diff --git a/eddy_tracker_configuration.yaml b/eddy_tracker_configuration.yaml
@@ -24,6 +24,7 @@ DOMAIN:
 AVISO:
   AVISO_DT14: Yes # Use new AVISO DT14 data (e.g., Capet et al, 2014)
   AVISO_FILES: 'dt_global_twosat_msla_h_????????_20140106.nc'
+  #AVISO_FILES: 'dt_ref_global_merged_msla_h_qd_????????_????????_????????.nc'
   AVISO_DT14_SUBSAMP: No # subsample the daily AVISO DT14
   DAYS_BTWN_RECORDS: 7  # sampling rate (days)
 
@@ -76,7 +77,7 @@ SMOOTHING_Q:
   SMOOTH_FAC: 5 # number of passes
 
 # Min and max permitted eddy radii [degrees] and amplitudes (cm)
-RADMIN: 0.35 # CSS11 recommend ~50 km minimum
+RADMIN: 0.35 # 0.4 for DT10, 0.35 for DT14
 RADMAX: 4.461
 AMPMIN: 1.
 AMPMAX: 150.
@@ -108,7 +109,7 @@ EVOLVE_AREA_MAX: 500 # max change in area
 
 
 # Make figures for animations
-SAVE_FIGURES: Yes
+SAVE_FIGURES: No
 
 
 # Define track_extra_variables to track and save:
@@ -118,9 +119,10 @@ SAVE_FIGURES: Yes
 # - profiles of swirl velocity from effective contour inwards
 # Useful for working with ARGO data    
 TRACK_EXTRA_VARIABLES: No
-    
-# INTERP_METHOD: 'RectBivariate'
-INTERP_METHOD: 'griddata'
+
+# 'RectBivariate' is faster
+INTERP_METHOD: 'RectBivariate'
+#INTERP_METHOD: 'griddata'
     
     
     
diff --git a/make_eddy_track_AVISO.py b/make_eddy_track_AVISO.py
@@ -30,8 +30,8 @@
 Scroll down to line ~640 to get started
 ===============================================================================
 """
-from matplotlib import use as mpl_use
-mpl_use('Agg')
+#from matplotlib import use as mpl_use
+#mpl_use('Agg')
 import sys
 import glob as glob
 from py_eddy_tracker_classes import plt, np, dt, Dataset, time, \
@@ -520,9 +520,10 @@ def __init__(self, AVISO_FILE, THE_DOMAIN, PRODUCT,
 
     def __getstate__(self):
         """
-        Needed for Pickle
+        Remove refernces to unwanted attributes in self.
+        This reduces the size of saved cPickle objects.
         """
-        print '--- removing unwanted attributes'
+        #print '--- removing unwanted attributes'
         pops = ('Mx', 'My', '_f', '_angle', '_dx', '_dy', '_gof', '_lon',
                 '_lat', '_pm', '_pn', '_umask', '_vmask', 'eke', 'u', 'v',
                 'mask')
@@ -1232,7 +1233,7 @@ def set_interp_coeffs(self, sla, uspd):
 
             # Track the eddies
             A_eddy = track_eddies(A_eddy, first_record)
-            C_eddy = track_eddies(C_eddy, first_record)
+            #C_eddy = track_eddies(C_eddy, first_record)
 
             if SAVE_FIGURES:  # Make figures for animations
 
@@ -1259,7 +1260,9 @@ def set_interp_coeffs(self, sla, uspd):
                     print '+++'
 
                 A_eddy.write2netcdf(rtime)
-                C_eddy.write2netcdf(rtime)
+                #print '______________________eddy.index', A_eddy.index
+                #print '______________________eddy.ncind', A_eddy.ncind
+                #C_eddy.write2netcdf(rtime)
 
         if str(DATE_END) in AVISO_FILE:
             active = False
diff --git a/make_eddy_track_CLS.py b/make_eddy_track_CLS.py
@@ -41,83 +41,94 @@
 
 if __name__ == '__main__':
 
-    DATA_DIR = '/data/OCE_ETU/MSA/emason/Global_DT10/'
+    TRACK_DURATION_MIN = 28
+    
+    #DATA_DIR = '/data/OCE_ETU/MSA/emason/Global_DT10/'
+    DATA_DIR = '/Users/emason/mercurial_projects/py-eddy-tracker-cls/py-eddy-tracker-cls/outputs/'
 
     A_PKL_FILES = 'A_eddy_????????.pkl'
-    C_PKL_FILES = 'C_eddy_????????.pkl'
+    #C_PKL_FILES = 'C_eddy_????????.pkl'
 
-    A_PKL_FILES = sorted(glob.glob(DATA_DIR + A_PKL_FILES))
-    C_PKL_FILES = sorted(glob.glob(DATA_DIR + C_PKL_FILES))
+    PKL_FILES = glob.glob(DATA_DIR + A_PKL_FILES)
+    PKL_FILES.sort()
+    #C_PKL_FILES = glob.glob(DATA_DIR + C_PKL_FILES)
 
-    for PKL_FILES in (A_PKL_FILES, C_PKL_FILES):
+    for active, PKL_FILE in enumerate(PKL_FILES):
 
-        for active, PKL_FILE in enumerate(PKL_FILES):
+        print PKL_FILE
+        
+        try: del eddy
+        except: pass
+        
+        # Unpickle
+        with open(PKL_FILE, 'rb') as the_pickle:
+            eddy = pickle.load(the_pickle)
+            print '--- loaded %s' % PKL_FILE
 
-            print PKL_FILE
+        #print 'eddy.index', eddy.index
+        eddy.savedir = DATA_DIR + eddy.savedir.rpartition('/')[-1]
 
-            # Unpickle
-            with open(PKL_FILE, 'rb') as the_pickle:
-                eddy = pickle.load(the_pickle)
-                print '--- loaded %s' % PKL_FILE
+        eddy.TRACK_DURATION_MIN = TRACK_DURATION_MIN
 
-            #print 'eddy.index', eddy.index
-            eddy.savedir = DATA_DIR + eddy.savedir.rpartition('/')[-1]
+        if active:
 
-            eddy.TRACK_DURATION_MIN = 10
-
-            if active:
-
-                first_record = False
-
-                eddy.new_list = False
-                eddy.tracklist = tracklist.tolist()
+            eddy.new_list = False
+            eddy.tracklist = tracklist#.tolist()
 
+            try:
                 eddy.index = index
                 eddy.ch_index = ch_index
+                #eddy.ncind = ncind
+            except:
+                pass
+            eddy.old_lon = old_lon
+            eddy.old_lat = old_lat
+            eddy.old_amp = old_amp
+            eddy.old_uavg = old_uavg
+            eddy.old_radii_s = old_radii_s
+            eddy.old_radii_e = old_radii_e
+            eddy.old_teke = old_teke
+
+            eddy = track_eddies(eddy, first_record)
+            
+            try:
                 eddy.ncind = ncind
-
-                eddy.old_lon = old_lon
-                eddy.old_lat = old_lat
-                eddy.old_amp = old_amp
-                eddy.old_uavg = old_uavg
-                eddy.old_radii_s = old_radii_s
-                eddy.old_radii_e = old_radii_e
-                eddy.old_teke = old_teke
-                #eddy.old_temp = old_temp
-                #eddy.old_salt = old_salt
-
-                eddy = track_eddies(eddy, first_record)
-                tracklist = np.copy(eddy.tracklist)
-
-            else:
-
-                first_record = True
-
-                eddy.create_netcdf(DATA_DIR, eddy.savedir)
-                eddy.set_old_variables()
-
-                eddy = track_eddies(eddy, first_record)
-                tracklist = np.copy(eddy.tracklist)
-
-            print 'eddy.index', eddy.index
-            print len(eddy.tracklist)
-
-            if not first_record:
-                print eddy.new_time_tmp[0]
-                eddy.write2netcdf(eddy.new_time_tmp[0])
-
-            index = eddy.index
-            ch_index = eddy.ch_index
-            ncind = eddy.ncind
-
-            old_lon = eddy.old_lon
-            old_lat = eddy.old_lat
-            old_amp = eddy.old_amp
-            old_uavg = eddy.old_uavg
-            old_radii_s = eddy.old_radii_s
-            old_radii_e = eddy.old_radii_e
-            old_teke = eddy.old_teke
-            #old_temp = eddy.old_temp
-            #old_salt = eddy.old_salt
-
-            #eddy.reset_holding_variables()
+            except:
+                pass
+
+        else:
+            
+            eddy.create_netcdf(DATA_DIR, eddy.savedir)
+            eddy.set_old_variables()
+            first_record = True
+            eddy = track_eddies(eddy, first_record)
+            first_record = False
+            tracklist = eddy.tracklist
+        
+        if not first_record:
+            
+            #print eddy.new_time_tmp[0]
+            
+            eddy.write2netcdf(eddy.new_time_tmp[0])
+            # tracklist is modified by write2netcdf, so
+            # place update just after
+            tracklist = eddy.tracklist
+            ncind = np.copy(eddy.ncind)
+            ch_index = np.copy(eddy.ch_index)
+            index = np.copy(eddy.index)
+
+        #print '______________________eddy.index', eddy.index
+        #print '______________________eddy.ncind', eddy.ncind
+
+        old_lon = eddy.old_lon
+        #print '______________________eddy.old_lon', eddy.old_lon
+        old_lat = eddy.old_lat
+        old_amp = eddy.old_amp
+        old_uavg = eddy.old_uavg
+        old_radii_s = eddy.old_radii_s
+        old_radii_e = eddy.old_radii_e
+        old_teke = eddy.old_teke
+        #old_temp = eddy.old_temp
+        #old_salt = eddy.old_salt
+
+        eddy.reset_holding_variables()
diff --git a/make_eddy_tracker_list_obj.py b/make_eddy_tracker_list_obj.py
@@ -182,14 +182,13 @@ class Track (object):
         dayzero - True at first appearance of eddy
     """
 
-    def __init__(self, eddy_index, DATATYPE, lon, lat, time, uavg, teke,
+    def __init__(self, DATATYPE, lon, lat, time, uavg, teke,
                  radius_s, radius_e, amplitude, temp=None, salt=None,
                  save_extras=False, contour_e=None, contour_s=None,
                  uavg_profile=None, shape_error=None):
 
-        self.eddy_index = eddy_index
+        #self.eddy_index = eddy_index
         self.DATATYPE   = DATATYPE
-        #self.lon        = np.atleast_1d(lon)
         self.lon        = [lon]
         self.lat        = [lat]
         self.ocean_time = [time]
@@ -264,7 +263,7 @@ class TrackList (object):
         qparameter: Q parameter range used for contours
         new_lon, new_lat: new lon/lat centroids
         old_lon, old_lat: old lon/lat centroids
-        eddy_index:   index of eddy in track_list
+        index:   index of eddy in track_list
     """
     #def __init__(self, DATATYPE, TRACK_DURATION_MIN, TRACK_EXTRA_VARIABLES):
     def __init__(self, DATATYPE, SIGN_TYPE, SAVE_DIR, grd, search_ellipse,
@@ -339,28 +338,29 @@ def __init__(self, DATATYPE, SIGN_TYPE, SAVE_DIR, grd, search_ellipse,
         self.sla = None
         self.slacopy = None
 
-        self.new_lon    = [] #np.array([])
-        self.new_lat    = []
+        self.new_lon = []
+        self.new_lat = []
         self.new_radii_s = []
         self.new_radii_e = []
-        self.new_amp    = []
-        self.new_uavg   = []
-        self.new_teke   = []
+        self.new_amp = []
+        self.new_uavg = []
+        self.new_teke = []
         if 'ROMS' in self.DATATYPE:
             self.new_temp = []
             self.new_salt = []
-        self.new_time   = []
-        self.old_lon    = []
-        self.old_lat    = []
-        self.old_radii_s  = []
-        self.old_radii_e  = []
-        self.old_amp    = []
-        self.old_uavg   = []
-        self.old_teke   = []
+        # NOTE check if new_time and old_time are necessary... 
+        self.new_time = []
+        self.old_lon = []
+        self.old_lat = []
+        self.old_radii_s = []
+        self.old_radii_e = []
+        self.old_amp = []
+        self.old_uavg = []
+        self.old_teke = []
         if 'ROMS' in self.DATATYPE:
             self.old_temp = []
             self.old_salt = []
-        self.old_time   = []
+        self.old_time = []
         if self.TRACK_EXTRA_VARIABLES:
             self.new_contour_e = []
             self.new_contour_s = []
@@ -386,7 +386,7 @@ def __getstate__(self):
         """
         Needed for Pickle
         """
-        print '--- removing unwanted attributes'
+        #print '--- removing unwanted attributes'
         pops = ('uspd', 'uspd_coeffs', 'sla_coeffs', 'points',
                 'circlon', 'circlat', 'sla', 'slacopy')
         result = self.__dict__.copy()
@@ -401,7 +401,7 @@ def add_new_track(self, lon, lat, time, uavg, teke,
         """
         Append a new 'track' object to the list
         """
-        self.tracklist.append(Track(self.index, self.DATATYPE,
+        self.tracklist.append(Track(self.DATATYPE,
                                     lon, lat, time, uavg, teke,
                                     radius_s, radius_e, amplitude,
                                     temp, salt, self.TRACK_EXTRA_VARIABLES,
@@ -841,11 +841,14 @@ def write2netcdf(self, rtime):
                             self.ncind += tsize
                             self.ch_index += 1
                             nc.sync()
-
+        
+        #print 'eddy.ncind', self.ncind
+        
         # Get index to first currently active track
         try:
             lasti = self.get_active_tracks(rtime)[0]
-
+            #print '__________________lasti', lasti
+            #print '__________________rtime', rtime
         except Exception:
             lasti = None
 
@@ -859,6 +862,7 @@ def write2netcdf(self, rtime):
 
         # Update old_lon and old_lat...
         self.old_lon = self.new_lon[lasti:]
+        #print '______________________self.old_lon', self.old_lon
         self.old_lat = self.new_lat[lasti:]
         self.old_radii_s = self.new_radii_s[lasti:]
         self.old_radii_e = self.new_radii_e[lasti:]
diff --git a/py_eddy_tracker_classes.py b/py_eddy_tracker_classes.py
@@ -877,7 +877,9 @@ def track_eddies(Eddy, first_record):
 
     haversine.distance_matrix(X_old, X_new, dist_mat)
     dist_mat = np.ascontiguousarray(dist_mat)
-
+    
+    #print '+++++++++++++ X_old.mean()', X_old.mean()
+    
     dist_mat_copy = dist_mat.copy()
 
     # *new_eddy_inds* contains indices to every newly identified eddy
@@ -1148,7 +1150,7 @@ def track_eddies(Eddy, first_record):
     # Now we need to add new eddies defined by new_eddy_inds
     if np.any(new_eddy_inds):
 
-        if True:
+        if False:
             print '------adding %s new eddies' % new_eddy_inds.sum()
 
         for neind, a_new_eddy in enumerate(new_eddy_inds):
