Updates to SwirlSpeed object

Author: evanmason

eddy_tracker_configuration.yaml

@@ -12,12 +12,12 @@ DOMAIN:
   THE_DOMAIN: 'Global'
   #THE_DOMAIN = 'BlackSea'
   #THE_DOMAIN = 'MedSea' # not yet implemented
-  LONMIN: -70.
+  LONMIN: -40.
   LONMAX: -5.
-  LATMIN: 5.
-  LATMAX: 55.
+  LATMIN: 20.
+  LATMAX: 30.
   DATE_STR: 20020101
-  DATE_END: 20020630
+  DATE_END: 20020331
 
 # Specify the AVISO product
 AVISO:
@@ -37,6 +37,9 @@ PATHS:
   #RW_PATH: '/home/emason/data_tmp/chelton_eddies/rossrad.dat'
 
   # Path for saving of outputs
+#   SAVE_DIR: '/homelocal/emason/outputs/before_swirl/'
+#   SAVE_DIR: '/homelocal/emason/outputs/after_swirl/'
+#   SAVE_DIR: '/homelocal/emason/outputs/after_swirl_slice/'
   SAVE_DIR: '/homelocal/emason/outputs/'
   #SAVE_DIR: '/marula/emason/aviso_eddy_tracking/fd0cdb4b2bd9/'
 

make_eddy_track_AVISO.py

@@ -1168,7 +1168,7 @@ def set_interp_coeffs(self, sla, uspd):
             # 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'
             C_eddy.sign_type = 'Cyclonic'

make_eddy_tracker_list_obj.py

@@ -344,6 +344,7 @@ def find_nearest_contour(contcoll, x, y):
         dss = lc[:,1] - y
         dss **= 2
         ds += dss
+        #print 'ds', ds
         d = ds.min()
         if d < dmin:
             dmin = d

py_eddy_tracker_amplitude.py

@@ -226,70 +226,196 @@ def __init__(self, contour):
         ci : index to contours
         li : index to levels
         """
-        self.x = []
-        self.y = []
-        self.ci = []
-        self.li = []
+        x = []
+        y = []
+        ci = []
+        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])
+                x.append(coll.vertices[:, 0])
+                y.append(coll.vertices[:, 1])
                 thelen = len(coll.vertices[:, 0])
-                self.ci.append([cind] * thelen)
-                self.li.append([lind] * thelen)
+                ci.append(thelen)
+                #ci.append(coll.vertices.shape[1])
+            li.append(len(cont.get_paths()))
 
-        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.x = np.array([val for sublist in x for val in sublist])
+        self.y = np.array([val for sublist in y for val in sublist])
+        self.nb_pt_per_c = np.array(ci)
+        self.ci = self.nb_pt_per_c.cumsum() - self.nb_pt_per_c
+        self.nb_c_per_l = np.array(li)
+        self.li = self.nb_c_per_l.cumsum() - self.nb_c_per_l
 
         self.nearesti = None
-        self.dist = None
-        self.lbool = None
+        self.level_slice = None
 
 
-    def _set_level_bool(self, thelevel):
+    def _set_level_slice(self, thelevel):
         """
-        Set a boolean to select only values to a discrete
-        contour level in self.li
+        Set slices
         """
-        self.lbool = self.li == thelevel
+        if self.nb_c_per_l[thelevel] == 0:
+            self.level_slice = None
+        else:
+            self.level_view_of_contour = slice(self.li[thelevel],
+                self.li[thelevel] + self.nb_c_per_l[thelevel])
+            nb_pts = self.nb_pt_per_c[self.level_view_of_contour].sum()
+            self.level_slice = slice(self.ci[self.li[thelevel]],
+                self.ci[self.li[thelevel]] + nb_pts)
         return self
 
 
     def set_dist_array_size(self, thelevel):
         """
         """
-        self._set_level_bool(thelevel)
-        self.dist = np.empty_like(self.li[self.lbool])
+        self._set_level_slice(thelevel)
         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
+        #print self.x[self.level_slice]
+        #print self.y[self.level_slice]
+        self.dist = (self.x[self.level_slice] - xpt)**2
+        self.dist += (self.y[self.level_slice] - ypt)**2
+        #print 'self.dist',self.dist
         try:
             self.nearesti = self.dist.argmin()
         except:
+            #print '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]]
+        if self.nearesti is not None:
+            #print self.level_view_of_contour
+            indices_of_first_pts = self.ci[self.level_view_of_contour]
+            for i, index_of_first_pt in enumerate(indices_of_first_pts):
+                if (index_of_first_pt - indices_of_first_pts[0]) > self.nearesti:
+                    return i - 1
+            #print '/////'
+            return i
         else:
+            #print '-----'
             return False
 
 
+    #def get_uavg(self, Eddy, CS, collind, centlon_e, centlat_e, poly_eff,
+                #grd, eddy_radius_e, properties, save_all_uavg=False):
+        #"""
+        #Calculate geostrophic speed around successive contours
+        #Returns the average
+        
+        #If save_all_uavg == True we want uavg for every contour
+        #"""
+        ## Unpack indices for convenience
+        #imin, imax, jmin, jmax = Eddy.imin, Eddy.imax, Eddy.jmin, Eddy.jmax
+        
+        #points = np.array([grd.lon()[jmin:jmax, imin:imax].ravel(),
+                           #grd.lat()[jmin:jmax, imin:imax].ravel()]).T
+        
+        ## First contour is the outer one (effective)
+        #theseglon, theseglat = poly_eff.vertices[:, 0].copy(), \
+                            #poly_eff.vertices[:, 1].copy()
+        
+        #theseglon, theseglat = eddy_tracker.uniform_resample(
+            #theseglon, theseglat, method='akima')
+        
+        #uavg = Eddy.uspd_coeffs.ev(theseglat[:-1], theseglon[:-1]).mean()
 
+        #if save_all_uavg:
+            #all_uavg = [uavg]
+            #pixel_min = 1 # iterate until 1 pixel
 
-        
+        #else:
+            ## Iterate until PIXEL_THRESHOLD[0] number of pixels
+            #pixel_min = Eddy.PIXEL_THRESHOLD[0]
+        
+        ##start = True
+        #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)
+            
+            
+            ## 
+            #self.swirl.set_dist_array_size(citer.index)
+            #self.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]
+            
+                ###print poly_ii is poly_i
+            
+            ##if poly_i is not None:
+                
+                
+                ## 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])]):
+                #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]):
+                    
+                            #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()
+        
+        #try: # Assuming interior contours have been found
+            
+            #cx, cy = Eddy.M(theseglon, theseglat)
+            ## Speed based eddy radius (eddy_radius_s)
+            #centx_s, centy_s, eddy_radius_s, junk = fit_circle(cx, cy)
+            #centlon_s, centlat_s = Eddy.M.projtran(centx_s, centy_s, inverse=True)
+            #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 not save_all_uavg:
+                #return (uavg, centlon_e, centlat_e, eddy_radius_e,
+                        #theseglon, theseglat, theseglon, theseglat)
+            #else:
+                #return (uavg, centlon_e, centlat_e, eddy_radius_e,
+                        #theseglon, theseglat, theseglon, theseglat, all_uavg)
+
+            
+