merge

Author: adelepoulle

.hgignore

@@ -0,0 +1,10 @@
+syntax: glob
+*.pyc
+*.so
+*.c
+*.egg-*
+*.eclipse.*
+syntax: regexp
+^build$
+syntax: regexp
+^dist$

README.md

@@ -6,7 +6,7 @@ To avoid problems with installation, use of the virtualenv Python virtual enviro
 
 Then use pip to install all dependencies (numpy, scipy, matplotlib, netCDF4, cython, pyproj, Shapely, ...), e.g.:
 
-* ** pip install cython **
+* ** pip install cython numpy matplotlib scipy netCDF4 shapely pyproj**
 
 Then run the following to install the eddy tracker:
 

setup.py

@@ -23,12 +23,16 @@
         # 'src/scripts/EddyTrackingFull'
         ],
     zip_safe=False,
-    cmdclass={
-        'build_ext': cython_build_ext,
-    },
+    cmdclass=dict(
+        build_ext=cython_build_ext,
+        ),
+    entry_points=dict(console_scripts=[
+        'EddyUpdate = py_eddy_tracker.update.__init__:eddy_update',
+        ]),
     ext_modules=[Extension("py_eddy_tracker.tools",
                            ["src/py_eddy_tracker/tools.pyx"],
                            include_dirs=[numpy.get_include()])],
+    package_data={'py_eddy_tracker.featured_tracking': ['*.nc', ]},
     setup_requires=[
         'numpy>=1.8'],
     install_requires=[

share/tracking.yaml

@@ -4,7 +4,6 @@ DIAGNOSTIC_TYPE: 'SLA'
 
 PATHS:
   # Files produces with EddyIdentification
-  FILES_PATTERN: /data/adelepoulle/Test/Test_eddy/20160330_brasilia_comparaison/brasil/Anticyclonic_*.nc
   FILES_PATTERN: /home/emason/toto/Anticyclonic_*.nc
   # Path and filename of Chelton et al (1998) Rossby radius data
   # Obtain file from:

src/py_eddy_tracker/featured_tracking/old_tracker_reference.py

@@ -20,7 +20,13 @@ def mask_function(self, other):
         # Compute Parameter of ellips
         y = self.basic_formula_ellips_major_axis(
             self.obs['lat'],
-            degrees=True)
+            degrees=True,
+            c0=1.05,
+            cmin=1.05,
+            cmax=1.5,
+            lat1=23,
+            lat2=5,
+            )
         # mask from ellips
         mask = self.shifted_ellipsoid_degrees_mask(
             other,

src/py_eddy_tracker/grid/__init__.py

@@ -33,7 +33,7 @@
 from pyproj import Proj
 from numpy import unique, array, unravel_index, r_, floor, interp, arange, \
     sin, cos, deg2rad, arctan2, sqrt, pi, zeros, reciprocal, ma, empty, \
-    concatenate
+    concatenate, bytes_
 import logging
 from ..tracking_objects import nearest
 from re import compile as re_compile
@@ -48,16 +48,35 @@ def browse_dataset_in(data_dir, files_model, date_regexp, date_model,
     full_path = join_path(data_dir, files_model)
     logging.info('Search files : %s', full_path)
 
-    dataset_list = array(glob(full_path),
+    filenames = bytes_(glob(full_path))
+    dataset_list = empty(len(filenames),
                          dtype=[('filename', 'S256'),
                                 ('date', 'datetime64[D]'),
                                 ])
+    dataset_list['filename'] = bytes_(glob(full_path))
 
     logging.info('%s grids available', dataset_list.shape[0])
+    mode_attrs = False
+    if '(' not in date_regexp:
+        logging.debug('Attrs date : %s', date_regexp)
+        mode_attrs = date_regexp.strip().split(':')
+    else:
+        logging.debug('Pattern date : %s', date_regexp)
+
     for item in dataset_list:
-        result = pattern_regexp.match(item['filename'])
-        if result:
-            str_date = result.groups()[0]
+        str_date = None
+        if mode_attrs:
+            with Dataset(item['filename'].decode("utf-8")) as h:
+                if len(mode_attrs) == 1:
+                    str_date = getattr(h, mode_attrs[0])
+                else:
+                    str_date = getattr(h.variables[mode_attrs[0]], mode_attrs[1])
+        else:
+            result = pattern_regexp.match(str(item['filename']))
+            if result:
+                str_date = result.groups()[0]
+            
+        if str_date is not None:
             item['date'] = datetime.strptime(str_date, date_model).date()
 
     dataset_list.sort(order=['date', 'filename'])
@@ -71,6 +90,9 @@ def browse_dataset_in(data_dir, files_model, date_regexp, date_model,
         dataset_list = dataset_list[::sub_sampling_step]
 
     if start_date is not None or end_date is not None:
+        logging.info('Available grid from %s to %s',
+                     dataset_list[0]['date'],
+                     dataset_list[-1]['date'])
         logging.info('Filtering grid by time %s, %s', start_date, end_date)
         mask = (dataset_list['date'] >= start_date) * (
             dataset_list['date'] <= end_date)
@@ -176,14 +198,9 @@ def read_nc(self, varname, indices=slice(None)):
           varname : variable ('temp', 'mask_rho', etc) to read
           indices : slice
         """
-        with Dataset(self.grid_filename) as h_nc:
+        with Dataset(self.grid_filename.decode("utf-8")) as h_nc:
             return h_nc.variables[varname][indices]
 
-    @property
-    def nc_variables(self):
-        with Dataset(self.grid_filename) as h_nc:
-            return h_nc.variables.keys()
-
     @property
     def view(self):
         return (self.slice_j,
@@ -207,7 +224,7 @@ def read_nc_att(self, varname, att):
           varname : variable ('temp', 'mask_rho', etc) to read
           att : string of attribute, eg. 'valid_range'
         """
-        with Dataset(self.grid_filename) as h_nc:
+        with Dataset(self.grid_filename.decode("utf-8")) as h_nc:
             return getattr(h_nc.variables[varname], att)
 
     @property
@@ -416,27 +433,14 @@ def vv2vr(vv_in, m_p, l_p):
         mshp, lshp = vv_in.shape
         return vv2vr(vv_in, mshp + 1, lshp)
 
-    def rho2u_2d(self, rho_in):
-        """
-        Convert a 2D field at rho points to a field at u_val points
-        """
-        assert rho_in.ndim == 2, 'rho_in must be 2d'
-        return ((rho_in[:, :-1] + rho_in[:, 1:]) / 2).squeeze()
-
-    def rho2v_2d(self, rho_in):
-        """
-        Convert a 2D field at rho points to a field at v_val points
-        """
-        assert rho_in.ndim == 2, 'rho_in must be 2d'
-        return ((rho_in[:-1] + rho_in[1:]) / 2).squeeze()
-
     def uvmask(self):
         """
         Get mask at U and V points
         """
         logging.info('--- Computing umask and vmask for padded grid')
-        self._umask = self.mask[:, :-1] * self.mask[:, 1:]
-        self._vmask = self.mask[:-1] * self.mask[1:]
+        if getattr(self, 'mask', None) is not None:
+            self._umask = self.mask[:, :-1] * self.mask[:, 1:]
+            self._vmask = self.mask[:-1] * self.mask[1:]
 
     def set_basemap(self, with_pad=True):
         """
@@ -505,10 +509,11 @@ def uspd(self):
         """Get scalar speed
         """
         uspd = (self.u_val ** 2 + self.v_val ** 2) ** .5
-        if hasattr(uspd, 'mask'):
-            uspd.mask += self.mask[self.view_unpad]
-        else:
-            uspd = ma.array(uspd, mask=self.mask[self.view_unpad])
+        if self.mask is not None:
+            if hasattr(uspd, 'mask'):
+                uspd.mask += self.mask[self.view_unpad]
+            else:
+                uspd = ma.array(uspd, mask=self.mask[self.view_unpad])
         return uspd
 
     def set_interp_coeffs(self, sla, uspd):

src/py_eddy_tracker/grid/aviso.py

@@ -31,7 +31,8 @@
 from scipy import spatial
 from dateutil import parser
 from numpy import meshgrid, zeros, array, where, ma, argmin, vstack, ones, \
-    newaxis, sqrt, diff, r_
+    newaxis, sqrt, diff, r_, arange
+from scipy.interpolate  import interp1d
 import logging
 from netCDF4 import Dataset
 
@@ -60,6 +61,8 @@ class AvisoGrid(BaseData):
         'fillval',
         '_angle',
         'sla_coeffs',
+        'xinterp',
+        'yinterp',
         'uspd_coeffs',
         '__lon',
         '__lat',
@@ -77,22 +80,37 @@ def __init__(self, aviso_file, the_domain,
         super(AvisoGrid, self).__init__()
         logging.info('Initialising the *AVISO_grid*')
         self.grid_filename = aviso_file
-        self.domain = the_domain
-        self.lonmin = float(lonmin)
-        self.lonmax = float(lonmax)
-        self.latmin = float(latmin)
-        self.latmax = float(latmax)
-        self.grid_filename = aviso_file
-        
+
         self.lon_name = lon_name
         self.lat_name = lat_name
         self.grid_name = grid_name
 
         self._lon = self.read_nc(self.lon_name)
         self._lat = self.read_nc(self.lat_name)
+        if the_domain is None:
+            self.domain = 'Automatic Domain'
+            dlon = abs(self._lon[1] - self._lon[0])
+            dlat = abs(self._lat[1] - self._lat[0])
+            self.lonmin = float(self._lon.min()) + dlon * 2
+            self.lonmax = float(self._lon.max()) - dlon * 2
+            self.latmin = float(self._lat.min()) + dlat * 2
+            self.latmax = float(self._lat.max()) - dlat * 2
+            if ((self._lon[-1] + dlon) % 360) == self._lon[0]:
+                self.domain = 'Global'
+                self.lonmin = -100.
+                self.lonmax = 290.
+                self.latmin = -80.
+                self.latmax = 80.
+        else:
+            self.domain = the_domain
+            self.lonmin = float(lonmin)
+            self.lonmax = float(lonmax)
+            self.latmin = float(latmin)
+            self.latmax = float(latmax)
+
         self.fillval = self.read_nc_att(self.grid_name, '_FillValue')
 
-        if lonmin < 0 and lonmax <= 0:
+        if self.lonmin < 0 and self.lonmax <= 0:
             self._lon -= 360.
         self._lon, self._lat = meshgrid(self._lon, self._lat)
         self._angle = zeros(self._lon.shape)
@@ -102,7 +120,7 @@ def __init__(self, aviso_file, the_domain,
 
         # zero_crossing, used for handling a longitude range that
         # crosses zero degree meridian
-        if lonmin < 0 and lonmax >= 0 and 'MedSea' not in self.domain:
+        if self.lonmin < 0 and self.lonmax >= 0 and 'MedSea' not in self.domain:
             if ((self.lonmax < self._lon.max()) and (self.lonmax > self._lon.min()) and (self.lonmin < self._lon.max()) and (self.lonmin > self._lon.min())):
                 pass
             else:
@@ -119,9 +137,15 @@ def __init__(self, aviso_file, the_domain,
         self.get_aviso_f_pm_pn()
         self.set_u_v_eke()
         self.shape = self.lon.shape
-#         pad2 = 2 * self.pad
-#         self.shape = (self.f_coriolis.shape[0] - pad2,
-#                       self.f_coriolis.shape[1] - pad2)
+
+        # self.init_pos_interpolator()
+
+    def init_pos_interpolator(self):
+        self.xinterp = interp1d(self.lon[0].copy(), arange(self.lon.shape[1]), assume_sorted=True, copy=False, fill_value=(0, -1), bounds_error=False, kind='nearest')
+        self.yinterp = interp1d(self.lat[:, 0].copy(), arange(self.lon.shape[0]), assume_sorted=True, copy=False, fill_value=(0, -1), bounds_error=False, kind='nearest')
+
+    def nearest_indice(self, lon, lat):
+        return self.xinterp(lon), self.yinterp(lat)
 
     def set_filename(self, file_name):
         self.grid_filename = file_name
@@ -137,7 +161,7 @@ def get_aviso_data(self, aviso_file, dimensions=None):
         if units not in self.KNOWN_UNITS:
             raise Exception('Unknown units : %s' % units)
 
-        with Dataset(self.grid_filename) as h_nc:
+        with Dataset(self.grid_filename.decode('utf-8')) as h_nc:
             grid_dims = array(h_nc.variables[self.grid_name].dimensions)
             lat_dim = h_nc.variables[self.lat_name].dimensions[0]
             lon_dim = h_nc.variables[self.lon_name].dimensions[0]
@@ -186,44 +210,6 @@ def set_mask(self, sla):
                 sea_label = self.labels[plus9, plus200]
                 self.mask += self.labels != sea_label
 
-    def fillmask(self, data, mask):
-        """
-        Fill missing values in an array with an average of nearest
-        neighbours
-        From http://permalink.gmane.org/gmane.comp.python.scientific.user/19610
-        """
-        raise Exception('Use convolution to fill data')
-        assert data.ndim == 2, 'data must be a 2D array.'
-        fill_value = 9999.99
-        data[mask == 0] = fill_value
-
-        # Create (i, j) point arrays for good and bad data.
-        # Bad data are marked by the fill_value, good data elsewhere.
-        igood = vstack(where(data != fill_value)).T
-        ibad = vstack(where(data == fill_value)).T
-
-        # Create a tree for the bad points, the points to be filled
-        tree = spatial.cKDTree(igood)
-
-        # Get the four closest points to the bad points
-        # here, distance is squared
-        dist, iquery = tree.query(ibad, k=4, p=2)
-
-        # Create a normalised weight, the nearest points are weighted as 1.
-        #   Points greater than one are then set to zero
-        weight = dist / (dist.min(axis=1)[:, newaxis])
-        weight *= ones(dist.shape)
-        weight[weight > 1.] = 0.
-
-        # Multiply the queried good points by the weight, selecting only the
-        # nearest points. Divide by the number of nearest points to get average
-        xfill = weight * data[igood[:, 0][iquery], igood[:, 1][iquery]]
-        xfill = (xfill / weight.sum(axis=1)[:, newaxis]).sum(axis=1)
-
-        # Place average of nearest good points, xfill, into bad point locations
-        data[ibad[:, 0], ibad[:, 1]] = xfill
-        return data
-
     @property
     def lon(self):
         if self.__lon is None: