Bugfix for missing anticyclones

Author: evanmason

eddy_tracker_configuration.yaml

@@ -10,6 +10,7 @@ DIAGNOSTIC_TYPE: 'SLA'
 # Specify the AVISO domain
 DOMAIN:
   THE_DOMAIN: 'Global'
+#   THE_DOMAIN: 'Regional'
   #THE_DOMAIN = 'BlackSea'
   #THE_DOMAIN = 'MedSea' # not yet implemented
   LONMIN: -40.

make_eddy_track_AVISO.py

@@ -145,7 +145,9 @@ def set_initial_indices(self, LONMIN, LONMAX, LATMIN, LATMAX):
         """
         Get indices for desired domain
         """
-        print '--- Setting initial indices to LONMIN, LONMAX, LATMIN, LATMAX'
+        print '--- Setting initial indices to %s domain' % self.THE_DOMAIN
+        print '------ LONMIN = %s, LONMAX = %s, LATMIN = %s, LATMAX = %s' % (
+                                           LONMIN, LONMAX, LATMIN, LATMAX)
         self.i0, junk = self.nearest_point(LONMIN, 
                                            LATMIN + 0.5 * (LATMAX - LATMIN))
         self.i1, junk = self.nearest_point(LONMAX,
@@ -199,7 +201,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
 
@@ -489,13 +491,14 @@ class AvisoGrid (PyEddyTracker):
     Class to satisfy the need of the eddy tracker
     to have a grid class
     """
-    def __init__(self, AVISO_FILE, LONMIN, LONMAX, LATMIN, LATMAX,
+    def __init__(self, AVISO_FILE, THE_DOMAIN, LONMIN, LONMAX, LATMIN, LATMAX,
                  with_pad=True):
         """
         Initialise the grid object
         """
         super(AvisoGrid, self).__init__()
         print '\nInitialising the AVISO_grid'
+        self.THE_DOMAIN = THE_DOMAIN
         self.LONMIN = LONMIN
         self.LONMAX = LONMAX
         self.LATMIN = LATMIN
@@ -766,10 +769,21 @@ def set_interp_coeffs(self, sla, uspd):
     DIAGNOSTIC_TYPE = config['DIAGNOSTIC_TYPE']
 
     config['THE_DOMAIN'] = config['DOMAIN']['THE_DOMAIN']
-    config['LONMIN'] = config['DOMAIN']['LONMIN']
-    config['LONMAX'] = config['DOMAIN']['LONMAX']
-    config['LATMIN'] = config['DOMAIN']['LATMIN']
-    config['LATMAX'] = config['DOMAIN']['LATMAX']
+    
+    # It is not recommended to change values given below
+    # for 'Global', 'BlackSea' or 'MedSea'...
+    if 'Global' in config['THE_DOMAIN']:
+        config['LONMIN'] = -100.
+        config['LONMAX'] = 290.
+        config['LATMIN'] = -80.
+        config['LATMAX'] = 80.
+    
+    elif 'Regional' in config['THE_DOMAIN']:
+        config['LONMIN'] = config['DOMAIN']['LONMIN']
+        config['LONMAX'] = config['DOMAIN']['LONMAX']
+        config['LATMIN'] = config['DOMAIN']['LATMIN']
+        config['LATMAX'] = config['DOMAIN']['LATMAX']
+    
     DATE_STR = config['DATE_STR'] = config['DOMAIN']['DATE_STR']
     DATE_END = config['DATE_END'] = config['DOMAIN']['DATE_END']
 
@@ -865,8 +879,9 @@ def set_interp_coeffs(self, sla, uspd):
         AVISO_FILES = AVISO_FILES[5:-5:np.int(DAYS_BTWN_RECORDS)]
 
     # Set up a grid object using first AVISO file in the list
-    sla_grd = AvisoGrid(AVISO_FILES[0], config['LONMIN'], config['LONMAX'],
-                                        config['LATMIN'], config['LATMAX'])
+    sla_grd = AvisoGrid(AVISO_FILES[0], config['THE_DOMAIN'],
+                        config['LONMIN'], config['LONMAX'],
+                        config['LATMIN'], config['LATMAX'])
 
     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])
@@ -897,9 +912,10 @@ def set_interp_coeffs(self, sla, uspd):
     #kwargs = config
     A_eddy = eddy_tracker.TrackList('AVISO', 'Anticyclonic', A_SAVEFILE,
                             sla_grd, search_ellipse, **config)
+    
     C_eddy = eddy_tracker.TrackList('AVISO', 'Cyclonic', C_SAVEFILE,
                             sla_grd, search_ellipse, **config)
-
+    
     A_eddy.search_ellipse = search_ellipse
     C_eddy.search_ellipse = search_ellipse
 
@@ -920,8 +936,8 @@ def set_interp_coeffs(self, sla, uspd):
     C_eddy.PIXEL_THRESHOLD = [PIXMIN, PIXMAX]
 
     # Create nc files for saving of eddy tracks
-    A_eddy.create_netcdf(DATA_DIR, A_SAVEFILE, 'Anticyclonic')
-    C_eddy.create_netcdf(DATA_DIR, C_SAVEFILE, 'Cyclonic')
+    A_eddy.create_netcdf(DATA_DIR, A_SAVEFILE)
+    C_eddy.create_netcdf(DATA_DIR, C_SAVEFILE)
 
 
     # Loop through the AVISO files...
@@ -1084,12 +1100,11 @@ def set_interp_coeffs(self, sla, uspd):
                 A_CS = ax.contour(sla_grd.lon(),
                                   sla_grd.lat(),
                                   A_eddy.sla, A_eddy.CONTOUR_PARAMETER)
-                # Note that CSc is for the cyclonics,
-                #   CONTOUR_PARAMETER in reverse order
+                # Note that C_CS is in reverse order
                 C_CS = ax.contour(sla_grd.lon(),
                                   sla_grd.lat(),
                                   C_eddy.sla, C_eddy.CONTOUR_PARAMETER)
-            
+                
             else:
                 Exception
 
@@ -1113,10 +1128,6 @@ def set_interp_coeffs(self, sla, uspd):
             C_eddy.swirl = SwirlSpeed(C_CS)
 
             # Now we loop over the CS collection
-            A_eddy.SIGN_TYPE = 'Anticyclonic'
-            C_eddy.SIGN_TYPE = 'Cyclonic'
-            
-            qqq
             if 'Q' in DIAGNOSTIC_TYPE:
                 A_eddy, C_eddy = collection_loop(CS, sla_grd, rtime,
                                    A_list_obj=A_eddy, C_list_obj=C_eddy,

make_eddy_tracker_list_obj.py

@@ -287,7 +287,7 @@ def __init__(self, DATATYPE, SIGN_TYPE, SAVE_DIR, grd, search_ellipse,
 
         self.DIAGNOSTIC_TYPE = kwargs.get('DIAGNOSTIC_TYPE', 'SLA')
 
-        self.THE_DOMAIN = kwargs.get('THE_DOMAIN', 'Global')
+        self.THE_DOMAIN = kwargs.get('THE_DOMAIN', 'Regional')
         self.LONMIN = np.float64(kwargs.get('LONMIN', -40))
         self.LONMAX = np.float64(kwargs.get('LONMAX', -30))
         self.LATMIN = np.float64(kwargs.get('LATMIN', 20))
@@ -304,10 +304,13 @@ def __init__(self, DATATYPE, SIGN_TYPE, SAVE_DIR, grd, search_ellipse,
 
         self.INTERP_METHOD = kwargs.get('INTERP_METHOD', 'RectBivariate')
         self.JDAY_REFERENCE = kwargs.get('JDAY_REFERENCE', 2448623.0)
+        
+        # NOTE: '.copy()' suffix is essential here
         self.CONTOUR_PARAMETER = kwargs.get('CONTOUR_PARAMETER',
-                                   np.linspace(-100., 101, 1))
+                                   np.arange(-100., 101, 1)).copy()
         if 'Cyclonic' in SIGN_TYPE:
-            self.CONTOUR_PARAMETER = self.CONTOUR_PARAMETER[::-1]
+            self.CONTOUR_PARAMETER *= -1
+        
         self.SHAPE_ERROR = kwargs.get('SHAPE_ERROR',
                               np.full(self.CONTOUR_PARAMETER.size, 55.))
 
@@ -514,7 +517,7 @@ def get_inactive_tracks(self, rtime):
 
 
 
-    def create_netcdf(self, directory, savedir, title,
+    def create_netcdf(self, directory, savedir,
                       grd=None, Ymin=None, Ymax=None,
                       Mmin=None, Mmax=None, model=None,
                       sigma_lev=None, rho_ntr=None):
@@ -526,7 +529,7 @@ def create_netcdf(self, directory, savedir, title,
         else:
            self.savedir = savedir.replace('.nc', '_ARGO_enabled.nc')
         nc = Dataset(self.savedir, 'w', format='NETCDF4')
-        nc.title = ''.join((title, ' eddy tracks'))
+        nc.title = ''.join((self.SIGN_TYPE, ' eddy tracks'))
         nc.directory = directory
         nc.DAYS_BTWN_RECORDS = np.float64(self.DAYS_BTWN_RECORDS)
         nc.TRACK_DURATION_MIN = np.float64(self.TRACK_DURATION_MIN)
@@ -1074,7 +1077,7 @@ def __init__(self, THE_DOMAIN, grd, RW_PATH=None):
         self.ZERO_CROSSING = grd.ZERO_CROSSING
         self.RW_PATH = RW_PATH
         self._tree = None
-        if self.THE_DOMAIN in ('Global', 'ROMS'):
+        if self.THE_DOMAIN in ('Global', 'Regional', 'ROMS'):
             assert self.RW_PATH is not None, \
                 'Must supply a path for the Rossby deformation radius data'
             data = np.loadtxt(RW_PATH)
@@ -1094,6 +1097,7 @@ def __init__(self, THE_DOMAIN, grd, RW_PATH=None):
 
     def __getstate__(self):
         """
+        Needed for Pickle
         """
         result = self.__dict__.copy()
         result.pop('_tree')
@@ -1102,15 +1106,18 @@ def __getstate__(self):
 
     def __setstate__(self, thedict):
         """
+        Needed for Pickle
         """
         self.__dict__ = thedict
         self._make_kdtree()
 
 
     def get_rwdistance(self, xpt, ypt, DAYS_BTWN_RECORDS):
         """
+        Return the distance required by SearchEllipse
+        to construct a search ellipse for eddy tracking.
         """
-        if self.THE_DOMAIN in ('Global', 'ROMS'):
+        if self.THE_DOMAIN in ('Global', 'Regional', 'ROMS'):
 	    #print 'xpt, ypt', xpt, ypt
             self.distance[:] = self._get_rlongwave_spd(xpt, ypt)
             self.distance *= 86400.
@@ -1152,6 +1159,8 @@ def get_rwdistance(self, xpt, ypt, DAYS_BTWN_RECORDS):
 
     def _make_subset(self):
         """
+        Make a subset of _defrad data over the domain.
+        If 'Global' is defined then widen the domain.
         """
         pad = 1.5 # degrees
         LONMIN, LONMAX, LATMIN, LATMAX = self.limits
@@ -1168,11 +1177,21 @@ def _make_subset(self):
         self._lon = self._lon[lloi]
         self._lat = self._lat[lloi]
         self._defrad = self._defrad[lloi]
+        
+        if 'Global' in self.THE_DOMAIN:
+            lloi = self._lon > 260.
+            self._lon = np.append(self._lon, self._lon[lloi] - 360.)
+            self._lat = np.append(self._lat, self._lat[lloi])
+            self._defrad = np.append(self._defrad, self._defrad[lloi])
+        
         self.x, self.y = self.M(self._lon, self._lat)
         return self
 
 
     def _make_kdtree(self):
+        """
+        Compute KDE tree for nearest indices.
+        """
         points = np.vstack([self.x, self.y]).T
         self._tree = spatial.cKDTree(points)
         return self
@@ -1207,8 +1226,25 @@ def _get_rlongwave_spd(self, xpt, ypt):
         return self.r_spd_long
 
 
-    
-    
+    def view_grid_subset(self):
+        """
+        Figure to check RossbyWaveSpeed grid after call to
+        self._make_subset()
+        To use, uncomment in SearchEllipse __init__ method
+        """
+        stride = 30
+        plt.figure()
+        ax = plt.subplot()
+        self.M.scatter(self.x, self.y, c='b')
+        self.M.drawcoastlines()
+        self.M.fillcontinents()
+        self.M.drawparallels(np.arange(-90, 90. + stride, stride),
+                             labels=[1, 0, 0, 0], ax=ax)
+        self.M.drawmeridians(np.arange(-360, 360. + stride, stride),
+                             labels=[0, 0, 0, 1], ax=ax)
+        plt.show()
+
+
 class SearchEllipse (object):
     """
     Class to construct a search ellipse/circle around a specified point.
@@ -1238,6 +1274,7 @@ def __init__(self, THE_DOMAIN, grd, DAYS_BTWN_RECORDS, RW_PATH=None):
         self.n_s_minor = self.DAYS_BTWN_RECORDS * 15e4 / 7.
         self.semi_n_s_minor = 0.5 * self.n_s_minor
         self.rwv = RossbyWaveSpeed(THE_DOMAIN, grd, RW_PATH=RW_PATH)
+        #self.rwv.view_grid_subset()
         self.rw_c = np.empty(1)
         self.rw_c_mod = np.empty(1)
         self.rw_c_fac = 1.75
@@ -1306,7 +1343,7 @@ def set_search_ellipse(self, xpt, ypt):
         self.xpt = xpt
         self.ypt = ypt
 
-        if self.THE_DOMAIN in ('Global', 'ROMS'):
+        if self.THE_DOMAIN in ('Global', 'Regional', 'ROMS'):
             self.rw_c[:] = self.rwv.get_rwdistance(xpt, ypt,
                                   self.DAYS_BTWN_RECORDS)
             self.rw_c_mod[:] = 1.75

py_eddy_tracker_classes.py

@@ -421,7 +421,7 @@ def get_uavg(Eddy, CS, collind, centlon_e, centlat_e, poly_eff,
                            poly_eff.vertices[:, 1].copy()
 
     theseglon, theseglat = eddy_tracker.uniform_resample(
-               theseglon, theseglat)
+               theseglon, theseglat, method='akima')
     if 'RectBivariate' in Eddy.INTERP_METHOD:
         uavg = Eddy.uspd_coeffs.ev(theseglat[1:], theseglon[1:]).mean()
 
@@ -480,7 +480,7 @@ def get_uavg(Eddy, CS, collind, centlon_e, centlat_e, poly_eff,
 
                         seglon, seglat = poly_i.vertices[:, 0], poly_i.vertices[:, 1]
                         seglon, seglat = eddy_tracker.uniform_resample(
-                                            seglon, seglat)
+                                            seglon, seglat, method='akima')
 
                         # Interpolate uspd to seglon, seglat, then get mean
                         if 'RectBivariate' in Eddy.INTERP_METHOD:
@@ -681,7 +681,7 @@ def collection_loop(CS, grd, rtime, A_list_obj, C_list_obj,
                                 # circumferential distribution
                                 contlon_e, contlat_e = \
                                     eddy_tracker.uniform_resample(contlon_e,
-                                                     contlat_e)
+                                                     contlat_e, method='akima')
 
                                 if 'Q' in Eddy.DIAGNOSTIC_TYPE:
                                     # Note, eddy amplitude == max(abs(vort/f)) within eddy, KCCMC11
@@ -1299,7 +1299,8 @@ def accounting(Eddy, old_ind, centlon, centlat,
             Eddy.insert_at_index('new_shape_error', Eddy.index, shape_error)
 
         if Eddy.new_list is True: # initialise a new list
-            print '------ starting a new track list for %ss' % Eddy.SIGN_TYPE
+            print ('------ starting a new track list for %s' % 
+                               Eddy.SIGN_TYPE.replace('nic', 'nes'))
             Eddy.new_list = False
 
         args = (centlon, centlat, rtime, uavg, teke,

py_eddy_tracker_property_classes.py

@@ -213,8 +213,8 @@ def _set_local_extrema(self, sign):
         """
         Set count of local SLA maxima/minima within eddy
         """
-        #local_extrema = np.ma.masked_where(self.mask == False, self.sla)
-        local_extrema = self.sla
+        local_extrema = np.ma.masked_where(self.mask == False, self.sla)
+        #local_extrema = self.sla
         local_extrema *= sign
         self._detect_local_minima(local_extrema)
         return self