Bugfix (incomplete) for Ellipse calculation

Author: evanmason

eddy_tracker_configuration.yaml

@@ -13,34 +13,37 @@ DOMAIN:
   #THE_DOMAIN = 'BlackSea'
   #THE_DOMAIN = 'MedSea' # not yet implemented
   LONMIN: -40.
-  LONMAX: -5.
-  LATMIN: 20.
-  LATMAX: 30.
-  DATE_STR: 20020101
-  DATE_END: 20020331
+  LONMAX: -30.
+  LATMIN: 30.
+  LATMAX: 40.
+  DATE_STR: 19980101
+  DATE_END: 19980331
 
 # 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_FILES: 'dt_global_allsat_msla_h_????????_20140106.nc'
   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: '/data/PVA/Externe/global/delayed-time/grids/msla/two-sat-merged/h/2002/'
+  DATA_DIR: '/Users/emason/data/AVISO/global/delayed-time/grids/msla/all-sat-merged/h/'
   #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: '/homelocal/emason/rossrad.dat'
+  RW_PATH: '/Users/emason/Downloads/rossrad.dat'
   #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: '/homelocal/emason/outputs/'
+  SAVE_DIR: '/Users/emason/toto/'
   #SAVE_DIR: '/marula/emason/aviso_eddy_tracking/fd0cdb4b2bd9/'
 
   # Reference Julian day (Julian date at Jan 1, 1992)

make_eddy_track_AVISO.py

@@ -853,8 +853,9 @@ def set_interp_coeffs(self, sla, uspd):
 
 
     # Instantiate search ellipse object
-    search_ellipse = eddy_tracker.SearchEllipse(THE_DOMAIN, DAYS_BTWN_RECORDS,
-                                    RW_PATH, [LONMIN, LONMAX, LATMIN, LATMAX])
+    search_ellipse = eddy_tracker.SearchEllipse(THE_DOMAIN, sla_grd,
+            DAYS_BTWN_RECORDS, sla_grd.ZERO_CROSSING, RW_PATH,
+            [LONMIN, LONMAX, LATMIN, LATMAX])
 
 
     if 'Gaussian' in SMOOTHING_TYPE:

make_eddy_tracker_list_obj.py

@@ -1203,50 +1203,58 @@ def reshape_mask_eff(self, grd):
 
 class RossbyWaveSpeed (object):
 
-    def __init__(self, domain, limits, rw_path=None):
+    def __init__(self, domain, grd, ZERO_CROSSING, limits, rw_path=None):
         """
         Initialise the RossbyWaveSpeedsobject
         """
         self.domain = domain
+        self.M = grd.M
+        self.ZERO_CROSSING = ZERO_CROSSING
         self.rw_path = rw_path
         if self.domain in ('Global', 'ROMS'):
             assert self.rw_path is not None, \
                 'Must supply a path for the Rossby deformation radius data'
             data = np.loadtxt(rw_path)
             self._lon = data[:, 1] - 360.
             self._lat = data[:, 0]
-            self._defrad = data[:, 3]
+            self._c = data[:, 2]
+            #self._defrad = data[:, 3]
             self.limits = limits
-            self._make_subset()
+            self._make_subset()._make_kdtree()
             self.vartype = 'variable'
         else:
             self.vartype = 'constant'
         self.distance = np.empty(1)
         self.start = True
 
 
-    def get_rwdistance(self, x, y, days_between_records):
+    def get_rwdistance(self, xpt, ypt, DAYS_BTWN_RECORDS):
         """
         """
         if self.domain in ('Global', 'ROMS'):
-	    #print 'x,y', x,y
-            self.distance[:] = self._get_rlongwave_spd(x, y)
+	    #print 'xpt, ypt', xpt, ypt
+            self.distance[:] = self._get_c(xpt, ypt)
             self.distance *= 86400.
             #if self.domain in 'ROMS':
                 #self.distance *= 1.5
         elif 'BlackSea' in self.domain:
             self.distance[:] = 15000. # e.g., Blokhina & Afanasyev, 2003
         elif 'MedSea' in self.domain:
             self.distance[:] = 20000.
-        else: Exception # Unknown domain
+        else:
+            Exception # Unknown domain
 
         if self.start:
             print '--------- setting ellipse for %s domain' %self.domain
-            print '--------- using %s Rossby deformation radius of %s m' \
-                                %(self.vartype, np.abs(self.distance[0]))
+            #print '--------- using %s Rossby deformation radius of %s m' \
+                                #%(self.vartype, np.abs(self.distance[0]))
+            print '--------- using %s baroclinic gravity wave phase speed of %s m/s' \
+                                %(self.vartype, self.distance[0] / 86400.)
+            print '--------- ', xpt, ypt
             self.start = False
-        self.distance *= days_between_records
-        return np.abs(self.distance)
+        self.distance *= DAYS_BTWN_RECORDS
+        #print 'aaaa', self.distance
+        return self.distance
 
 
     def _make_subset(self):
@@ -1260,77 +1268,91 @@ def _make_subset(self):
                                self._lat <= LATMAX + pad)
         self._lon = self._lon[lloi]
         self._lat = self._lat[lloi]
-        self._defrad = self._defrad[lloi]
-        self._make_kdtree()
+        self._c = self._c[lloi]
+        self.x, self.y = self.M(self._lon, self._lat)
+        return self
 
 
     def _make_kdtree(self):
-        points = np.vstack([self._lon, self._lat]).T
+        #points = np.vstack([self._lon, self._lat]).T
+        points = np.vstack([self.x, self.y]).T
         self._tree = spatial.cKDTree(points)
 
 
-    def _get_defrad(self, x, y):
+    #def _get_defrad(self, x, y):
+        #"""
+        #Get a point average of the deformation radius
+        #at x, y
+        #"""
+        #weights, i = self._tree.query(np.array([x, y]), k=4, p=2)
+        #weights /= weights.sum()
+        #self._weights = weights
+        #self.i = i
+        #return np.average(self._defrad[i], weights=weights)
+    
+    def _get_c(self, xpt, ypt):
         """
-        Get a point average of the deformation radius
-        at x, y
+        Get a point average of c, the phase speed, at x, y
         """
-        weights, i = self._tree.query(np.array([x, y]), k=4, p=2)
+        weights, i = self._tree.query(np.array([xpt, ypt]), k=4, p=2)
         weights /= weights.sum()
         self._weights = weights
         self.i = i
-        return np.average(self._defrad[i], weights=weights)
-    
+        return np.average(self._c[i], weights=weights)
 
-    def _get_rlongwave_spd(self, x, y):
-        """
-        Get the longwave phase speed, see Chelton etal (2008) pg 446:
-          c = -beta * defrad**2 (this only for extratropical waves...)
-        """
-        r_spd_long = self._get_defrad(x, y)
-        r_spd_long *= 1000. # km to m
-        r_spd_long **= 2
-        lat = np.average(self._lat[self.i], weights=self._weights)
-        beta = np.cos(np.deg2rad(lat))
-        beta *= 1458e-7 # 1458e-7 ~ (2 * 7.29*10**-5)
-        beta /= 6371315.0
-        r_spd_long *= -beta
-        return r_spd_long
+    #def _get_rlongwave_spd(self, x, y):
+        #"""
+        #Get the longwave phase speed, see Chelton etal (2008) pg 446:
+          #c = -beta * defrad**2 (this only for extratropical waves...)
+        #"""
+        ##r_spd_long = self._get_c(x, y)
+        ##r_spd_long *= 1000. # km to m
+        ##r_spd_long **= 2
+        ##lat = np.average(self._lat[self.i], weights=self._weights)
+        ##beta = np.cos(np.deg2rad(lat))
+        ##beta *= 1458e-7 # 1458e-7 ~ (2 * 7.29*10**-5)
+        ##beta /= 6371315.0
+        ##r_spd_long *= -beta
+        #return self._get_c(x, y)
 
 
 
 
 class SearchEllipse (object):
 
-    def __init__(self, domain, DAYS_BTWN_RECORDS, rw_path, limits):
+    def __init__(self, domain, grd, DAYS_BTWN_RECORDS, ZERO_CROSSING,
+                 rw_path, limits):
         """
         Class to construct a search ellipse/circle around a specified point.
         
         """
         self.domain = domain
         self.DAYS_BTWN_RECORDS = DAYS_BTWN_RECORDS
+        self.ZERO_CROSSING = ZERO_CROSSING
         self.rw_path = rw_path
         self.limits = limits
         self.rw_d_fac = 1.75
         self.e_w_major = self.DAYS_BTWN_RECORDS * 3e5 / 7.
-        self.n_s_minor = self.DAYS_BTWN_RECORDS * 15e4 / 7.
-        self.rwv = RossbyWaveSpeed(self.domain,
-                               self.limits, rw_path=self.rw_path)
+        self.n_s_minor = self.DAYS_BTWN_RECORDS * 2 * 15e4 / 7.
+        self.semi_n_s_minor = 0.5 * self.n_s_minor
+        self.rwv = RossbyWaveSpeed(domain, grd, ZERO_CROSSING,
+                               limits, rw_path=rw_path)
         self.rw_d = np.empty(1)
         self.rw_d_mod = np.empty(1)
 
 
     def _set_east_ellipse(self):
         """
         """
-        self.east_ellipse = patch.Ellipse((self.x, self.y),
+        self.east_ellipse = patch.Ellipse((self.xpt, self.ypt),
                                            self.e_w_major, self.n_s_minor)
         return self
 
 
     def _set_west_ellipse(self):
         """
         """
-        self.west_ellipse = patch.Ellipse((self.x, self.y),
+        self.west_ellipse = patch.Ellipse((self.xpt, self.ypt),
                                            self.rw_d_mod, self.n_s_minor)
         return self
 
@@ -1362,27 +1384,35 @@ def _set_black_sea_ellipse(self):
         return self
 
 
-    def set_search_ellipse(self, x, y):
+    def set_search_ellipse(self, xpt, ypt):
         """
         """
-        self.x = x
-        self.y = y
+        self.xpt = xpt
+        self.ypt = ypt
 
         if self.domain in ('Global', 'ROMS'):
-            self.rw_d[:] = self.rwv.get_rwdistance(x, y,
-                                                   self.DAYS_BTWN_RECORDS)
+            self.rw_d[:] = self.rwv.get_rwdistance(xpt, ypt,
+                                  self.DAYS_BTWN_RECORDS)
+            #print self.rw_d
+            #exit()
             self.rw_d_mod[:] = 1.75
             self.rw_d_mod *= self.rw_d
-            self.rw_d_mod[:] = np.maximum(self.rw_d_mod, self.n_s_minor)
+            self.rw_d_mod[:] = np.array([self.rw_d_mod,
+                                         self.semi_n_s_minor]).max()
+            
+            #print self.rw_d_mod, self.semi_n_s_minor
+            
             self.rw_d_mod *= 2.
+            #print self.rw_d_mod,
             self._set_global_ellipse()
         elif 'BlackSea'  in self.domain:
             self.rw_d_mod[:] = 1.75
-            self.rw_d[:] = self.rwv.get_rwdistance(x, y,
-                                                   self.DAYS_BTWN_RECORDS)
+            self.rw_d[:] = self.rwv.get_rwdistance(xpt, ypt,
+                                    self.DAYS_BTWN_RECORDS)
             self.rw_d_mod *= self.rw_d
             self._set_black_sea_ellipse()
-        else: Exception
+        else:
+            Exception
 
         return self
 
@@ -1392,12 +1422,18 @@ def view_search_ellipse(self, Eddy):
         Input A_eddy or C_eddy
         """
         plt.figure()
-        ax = plt.subplot(111)
-        ax.set_title('Rossby def. rad %s m' %self.rw_d[0])
-        Eddy.M.scatter(self.x, self.y, c='b')
+        ax = plt.subplot()
+        #ax.set_title('Rossby def. rad %s m' %self.rw_d[0])
+        Eddy.M.scatter(self.xpt, self.ypt, c='b')
         Eddy.M.plot(self.ellipse_path.vertices[:, 0],
                     self.ellipse_path.vertices[:, 1], 'r')
         Eddy.M.drawcoastlines()
+        stride = .5
+        Eddy.M.drawparallels(np.arange(-90, 90.+stride, stride),
+                             labels=[1, 0, 0, 0], ax=ax)
+        Eddy.M.drawmeridians(np.arange(-360, 360. + stride, stride),
+                             labels=[0, 0, 0, 1], ax=ax)
+        #plt.axis('image')
         plt.show()
 
 

py_eddy_tracker_classes.py

@@ -994,6 +994,9 @@ def track_eddies(Eddy, first_record):
             if 'ellipse' in Eddy.SEPARATION_METHOD:
                 if Eddy.search_ellipse.ellipse_path.contains_point(
                                           (new_x[new_ind], new_y[new_ind])):
+                #if Eddy.search_ellipse.ellipse_path.contains_point(
+                        #(Eddy.new_lon_tmp[new_ind], Eddy.new_lat_tmp[new_ind])):
+                    #Eddy.search_ellipse.view_search_ellipse(Eddy)
                     sep_proceed = True
                 else:
                     sep_proceed = False