Introduction of SwirlSpeed object

Author: emason

eddy_tracker_configuration.yaml

@@ -12,34 +12,33 @@ DOMAIN:
   THE_DOMAIN: 'Global'
   #THE_DOMAIN = 'BlackSea'
   #THE_DOMAIN = 'MedSea' # not yet implemented
-  LONMIN: -65. # -38.
-  LONMAX: -5.5 # -25.
-  LATMIN: 11.5 # 30.
-  LATMAX: 38.5 # 38.
-  DATE_STR: 19930101
-  DATE_END: 20131231
+  LONMIN: -70.
+  LONMAX: -5.
+  LATMIN: 5.
+  LATMAX: 55.
+  DATE_STR: 20020101
+  DATE_END: 20020630
 
 # Specify the AVISO product
 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_DT14_SUBSAMP: Yes # subsample the dqily AVISO DT14
+  AVISO_DT14_SUBSAMP: No # subsample the dqily AVISO DT14
   DAYS_BTWN_RECORDS: 7  # sampling rate (days)
 
 PATHS:
   # Path to AVISO data
-  #DATA_DIR: '/data/PVA/Externe/global/delayed-time/grids/msla/two-sat-merged/h/2002/'
-  DATA_DIR: '/marula/emason/data/altimetry/global/delayed-time/grids/msla/two-sat-merged/h/'
+  DATA_DIR: '/data/PVA/Externe/global/delayed-time/grids/msla/two-sat-merged/h/2002/'
+  #DATA_DIR: '/marula/emason/data/altimetry/global/delayed-time/grids/msla/two-sat-merged/h/'
   # Path and filename of Chelton et al (1998) Rossby radius data
   # Obtain file from:
   # http://www-po.coas.oregonstate.edu/research/po/research/rossby_radius/
-  #RW_PATH: '/homelocal/emason/rossrad.dat'
-  RW_PATH: '/home/emason/data_tmp/chelton_eddies/rossrad.dat'
+  RW_PATH: '/homelocal/emason/rossrad.dat'
+  #RW_PATH: '/home/emason/data_tmp/chelton_eddies/rossrad.dat'
 
   # Path for saving of outputs
-  #SAVE_DIR: '/homelocal/emason/outputs/'
-#   SAVE_DIR: '/marula/emason/aviso_eddy_tracking/fd0cdb4b2bd9/'
-  SAVE_DIR: '/marula/emason/aviso_eddy_tracking/fd0cdb4b2bd9/'
+  SAVE_DIR: '/homelocal/emason/outputs/'
+  #SAVE_DIR: '/marula/emason/aviso_eddy_tracking/fd0cdb4b2bd9/'
 
   # Reference Julian day (Julian date at Jan 1, 1992)
 JDAY_REFERENCE: 2448623.

make_eddy_track_AVISO.py

@@ -36,6 +36,7 @@
                                     datestr2datetime, gaussian_resolution, \
                                     get_cax, collection_loop, track_eddies, \
                                     anim_figure, interpolate
+from py_eddy_tracker_amplitude import SwirlSpeed
 import make_eddy_tracker_list_obj as eddy_tracker
 from dateutil import parser
 from mpl_toolkits.basemap import Basemap
@@ -153,7 +154,7 @@ def set_index_padding(self, pad=2):
         around 2d variables.
         Padded matrices are needed only for geostrophic velocity computation.
         """
-        print '--- Setting padding indices with pad=%s' %pad
+        print '--- Setting padding indices with PAD=%s' %pad
 
         self.pad = pad
 
@@ -1163,7 +1164,10 @@ def set_interp_coeffs(self, sla, uspd):
                 plt.axis('image')
                 plt.show()
 
-
+            # Set contour coordinates and indices for calculation of
+            # speed-based radius
+            A_eddy.swirl = SwirlSpeed(A_CS)
+            C_eddy.swirl = SwirlSpeed(C_CS)
 
             # Now we loop over the CS collection
             A_eddy.sign_type = 'Anticyclonic'

py_eddy_tracker_amplitude.py

@@ -214,4 +214,82 @@ def debug_figure(self):
         lmi_j = lmi_j[0] + self.grd.jmin
         x_i, y_i = self.grd.lon()[lmi_j, lmi_i], self.grd.lat()[lmi_j, lmi_i]
         plt.scatter(x_i, y_i, c='gray')
-        plt.show()
+        plt.show()
+
+        
+        
+class SwirlSpeed(object):
+    """
+    """
+    def __init__(self, contour):
+        """
+        ci : index to contours
+        li : index to levels
+        """
+        self.x = []
+        self.y = []
+        self.ci = []
+        self.li = []
+
+        for lind, cont in enumerate(contour.collections):
+            for cind, coll in enumerate(cont.get_paths()):
+                self.x.append(coll.vertices[:, 0])
+                self.y.append(coll.vertices[:, 1])
+                thelen = len(coll.vertices[:, 0])
+                self.ci.append([cind] * thelen)
+                self.li.append([lind] * thelen)
+    
+        self.x = np.array([val for sublist in self.x for val in sublist])
+        self.y = np.array([val for sublist in self.y for val in sublist])
+        self.ci = np.array([val for sublist in self.ci for val in sublist],
+                                                               dtype=np.int)
+        self.li = np.array([val for sublist in self.li for val in sublist],
+                                                               dtype=np.int)
+        
+        self.nearesti = None
+        self.dist = None
+        self.lbool = None
+        
+    
+    def _set_level_bool(self, thelevel):
+        """
+        Set a boolean to select only values to a discrete
+        contour level in self.li
+        """
+        self.lbool = self.li == thelevel
+        return self
+        
+    
+    def set_dist_array_size(self, thelevel):
+        """
+        """
+        self._set_level_bool(thelevel)
+        self.dist = np.empty_like(self.li[self.lbool])
+        return self
+        
+    
+    def set_nearest_contour_index(self, xpt, ypt):
+        """
+        """
+        self.dist[:] = (self.x[self.lbool] - xpt)**2
+        self.dist += (self.y[self.lbool] - ypt)**2
+        try:
+            self.nearesti = self.dist.argmin()
+        except:
+            self.nearesti = None
+        return self
+    
+    
+    def get_index_nearest_path(self):
+        """
+        """
+        ind = self.nearesti
+        if ind is not None:
+            return self.ci[self.lbool[ind]]
+        else:
+            return False
+    
+    
+        
+        
+        

py_eddy_tracker_classes.py

@@ -522,44 +522,56 @@ def get_uavg(Eddy, CS, collind, centlon_e, centlat_e, poly_eff,
     citer = np.nditer(CS.cvalues, flags=['c_index'])
     #print '************************************************'
     while not citer.finished:
+         
+        ## Get contour around centlon_e, centlat_e at level [collind:][iuavg]
+        #segi, poly_i = eddy_tracker.find_nearest_contour(
+                         #CS.collections[citer.index], centlon_e, centlat_e)
+         
+        
+        # 
+        Eddy.swirl.set_dist_array_size(citer.index)
+        Eddy.swirl.set_nearest_contour_index(centlon_e, centlat_e)
+        segi = Eddy.swirl.get_index_nearest_path()
+        
+        
+        if segi:
+            
+            poly_i = CS.collections[citer.index].get_paths()[segi]
 
-        # Get contour around centlon_e, centlat_e at level [collind:][iuavg]
-        segi, poly_i = eddy_tracker.find_nearest_contour(
-                         CS.collections[citer.index], centlon_e, centlat_e)
+            ##print poly_ii is poly_i
 
-        if poly_i is not None:
+        #if poly_i is not None:
 
-            # NOTE: contains_points requires matplotlib 1.3 +
-            mask_i_sum = poly_i.contains_points(points).sum()
 
             # 1. Ensure polygon_i is within polygon_e
             # 2. Ensure polygon_i contains point centlon_e, centlat_e
             # 3. Respect size range
-            if np.all([poly_eff.contains_path(poly_i),
-                       poly_i.contains_point([centlon_e, centlat_e]),
-                       (mask_i_sum >= pixel_min and
-                        mask_i_sum <= Eddy.PIXEL_THRESHOLD[1])]):
-                
-                seglon, seglat = poly_i.vertices[:, 0], poly_i.vertices[:, 1]
-                seglon, seglat = eddy_tracker.uniform_resample(
-                                      seglon, seglat, method='akima')
-                
-                # Interpolate uspd to seglon, seglat, then get mean
-                uavgseg = Eddy.uspd_coeffs.ev(seglat[:-1], seglon[:-1]).mean()
-                #uavgseg = uavgseg
-                #uavgseg = rbspline.ev(seglat[:-1], seglon[:-1]).mean()
-                #uavgseg =  uavgseg.mean()
-                #print uavgseg
-                #assert uavgseg == uavgsegsss, 'fffffffffffff'
-                
-                if save_all_uavg:
-                    all_uavg.append(uavgseg)
+            #if np.all([poly_eff.contains_path(poly_i),
+                       #poly_i.contains_point([centlon_e, centlat_e]),
+                       #(mask_i_sum >= pixel_min and
+                        #mask_i_sum <= Eddy.PIXEL_THRESHOLD[1])]):
+            if poly_i.contains_point([centlon_e, centlat_e]):
+                if poly_eff.contains_path(poly_i):
+                    # NOTE: contains_points requires matplotlib 1.3 +
+                    mask_i_sum = poly_i.contains_points(points).sum()
+                    if (mask_i_sum >= pixel_min and
+                        mask_i_sum <= Eddy.PIXEL_THRESHOLD[1]):
 
-                if (uavgseg >= uavg) and (mask_i_sum >= pixel_min):
-                    uavg = uavgseg.copy()
-                    theseglon, theseglat = seglon.copy(), seglat.copy()
-                
-                inner_seglon, inner_seglat = seglon.copy(), seglat.copy()
+                        seglon, seglat = poly_i.vertices[:, 0], poly_i.vertices[:, 1]
+                        seglon, seglat = eddy_tracker.uniform_resample(
+                                            seglon, seglat, method='akima')
+                        
+                        # Interpolate uspd to seglon, seglat, then get mean
+                        uavgseg = Eddy.uspd_coeffs.ev(seglat[:-1], seglon[:-1]).mean()
+                        
+                        if save_all_uavg:
+                            all_uavg.append(uavgseg)
+                        
+                        if uavgseg >= uavg:
+                            uavg = uavgseg.copy()
+                            theseglon, theseglat = seglon.copy(), seglat.copy()
+                        
+                        inner_seglon, inner_seglat = seglon.copy(), seglat.copy()
 
         citer.iternext()
 
@@ -747,7 +759,7 @@ def collection_loop(CS, grd, rtime, A_list_obj, C_list_obj,
         if VERBOSE:
             message = '------ doing collection %s, contour value %s'
             print message %  (collind, CS.cvalues[collind])
-                    
+        
         # Loop over individual CS contours (i.e., every eddy in field)
         for cont in coll.get_paths():