Change to be ok with python3 and to manage detection without information of area

Author: AntSimi

src/py_eddy_tracker/grid/__init__.py

@@ -5,7 +5,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
@@ -20,16 +20,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'])
@@ -43,6 +62,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)
@@ -148,14 +170,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,
@@ -179,7 +196,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
@@ -388,27 +405,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):
         """
@@ -477,10 +481,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

@@ -4,7 +4,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
 
@@ -33,6 +34,8 @@ class AvisoGrid(BaseData):
         'fillval',
         '_angle',
         'sla_coeffs',
+        'xinterp',
+        'yinterp',
         'uspd_coeffs',
         '__lon',
         '__lat',
@@ -50,22 +53,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)
@@ -75,7 +93,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:
@@ -92,9 +110,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
@@ -110,7 +134,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]
@@ -159,44 +183,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:

src/py_eddy_tracker/observations.py

@@ -110,10 +110,10 @@ def compute_pixel_path(self, x0, y0, x1, y1):
         # Delta index of y
         d_y = i_y1 - i_y0
 
-        d_max = maximum(abs(d_x), abs(d_y))
+        d_max = int_(maximum(abs(d_x), abs(d_y)))
 
         # Compute number of pixel which we go trought
-        nb_value = (abs(d_max) + 1).sum()
+        nb_value = int((abs(d_max) + 1).sum())
         # Create an empty array to store value of pixel across the travel
         # Max Index ~65000
         i_g = empty(nb_value, dtype='u2')
@@ -257,9 +257,12 @@ def merge(self, other):
         """
         nb_obs_self = len(self)
         nb_obs = nb_obs_self + len(other)
-        eddies = self.__class__(size=nb_obs)
-        eddies.obs[:nb_obs_self] = self.obs[:]
-        eddies.obs[nb_obs_self:] = other.obs[:]
+        eddies = self.new_like(self, nb_obs)
+        other_keys = other.obs.dtype.fields.keys()
+        for key in eddies.obs.dtype.fields.keys():
+            eddies.obs[key][:nb_obs_self] = self.obs[key][:]
+            if key in other_keys:
+                eddies.obs[key][nb_obs_self:] = other.obs[key][:]
         eddies.sign_type = self.sign_type
         return eddies
 
@@ -294,11 +297,7 @@ def insert_observations(self, other, index):
             return self
         if index < 0:
             index = self_size + index + 1
-        eddies = self.__class__(new_size,
-            track_extra_variables=self.track_extra_variables,
-            track_array_variables=self.track_array_variables,
-            array_variables=self.array_variables
-            )
+        eddies = self.new_like(self, new_size)
         eddies.obs[:index] = self.obs[:index]
         eddies.obs[index: index + insert_size] = other.obs
         eddies.obs[index + insert_size:] = self.obs[index:]
@@ -323,13 +322,21 @@ def distance(self, other):
             dist_result)
         return dist_result
 
+    @staticmethod
+    def new_like(eddies, new_size):
+        return eddies.__class__(new_size,
+            track_extra_variables=eddies.track_extra_variables,
+            track_array_variables=eddies.track_array_variables,
+            array_variables=eddies.array_variables
+            )
+
     def index(self, index):
         """Return obs from self at the index
         """
         size = 1
         if hasattr(index, '__iter__'):
             size = len(index)
-        eddies = self.__class__(size, self.track_extra_variables)
+        eddies = self.new_like(self, size)
         eddies.obs[:] = self.obs[index]
         return eddies
 
@@ -761,7 +768,10 @@ def create_variable(handler_nc, kwargs_variable, attr_variable,
             zlib=True,
             complevel=1,
             **kwargs_variable)
-        for attr, attr_value in attr_variable.iteritems():
+        attrs = list(attr_variable.keys())
+        attrs.sort()
+        for attr in attrs:
+            attr_value = attr_variable[attr]
             var.setncattr(attr, attr_value)
         if scale_factor is not None:
             var.scale_factor = scale_factor