update tracking

Author: adelepoulle

src/py_eddy_tracker/__init__.py

@@ -9,7 +9,7 @@ class ColoredFormatter(logging.Formatter):
         ERROR="\033[31;47m",
         WARNING="\033[30;47m",
         INFO="\033[36m",
-        DEBUG="\033[34m",
+        DEBUG="\033[34m\t",
         )
 
     def __init__(self, message):

src/py_eddy_tracker/observations.py

@@ -28,7 +28,8 @@
 ===========================================================================
 
 """
-from numpy import zeros, empty, nan, arange, interp, where, unique
+from numpy import zeros, empty, nan, arange, interp, where, unique, \
+    ma
 from netCDF4 import Dataset
 from py_eddy_tracker.tools import distance_matrix
 from . import VAR_DESCR, VAR_DESCR_inv
@@ -82,9 +83,11 @@ def __getitem__(self, attr):
         if attr in self.elements:
             return self.observations[attr]
         raise KeyError('%s unknown' % attr)
-    
+
     @property
     def dtype(self):
+        """Return dtype to build numpy array
+        """
         dtype = list()
         for elt in self.elements:
             dtype.append((elt, VAR_DESCR[elt][
@@ -94,6 +97,8 @@ def dtype(self):
 
     @property
     def elements(self):
+        """Return all variable name
+        """
         elements = [
             'lon',  # 'centlon'
             'lat',  # 'centlat'
@@ -113,9 +118,13 @@ def elements(self):
         return elements
 
     def coherence(self, other):
+        """Check coherence between two dataset
+        """
         return self.track_extra_variables == other.track_extra_variables
-    
+
     def merge(self, other):
+        """Merge two dataset
+        """
         nb_obs_self = len(self)
         nb_obs = nb_obs_self + len(other)
         eddies = self.__class__(size=nb_obs)
@@ -126,6 +135,8 @@ def merge(self, other):
 
     @property
     def obs(self):
+        """returan array observations
+        """
         return self.observations
 
     def __len__(self):
@@ -136,6 +147,8 @@ def __iter__(self):
             yield obs
 
     def insert_observations(self, other, index):
+        """Insert other obs in self at the index
+        """
         if not self.coherence(other):
             raise Exception('Observations with no coherence')
         insert_size = len(other.obs)
@@ -156,6 +169,8 @@ def insert_observations(self, other, index):
         return self
 
     def append(self, other):
+        """Merge
+        """
         return self + other
 
     def __add__(self, other):
@@ -172,13 +187,15 @@ def distance(self, other):
         return dist_result
 
     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.obs[:] = self.obs[index]
         return eddies
-    
+
     @staticmethod
     def load_from_netcdf(filename):
         with Dataset(filename) as h_nc:
@@ -187,22 +204,79 @@ def load_from_netcdf(filename):
             for variable in h_nc.variables:
                 if variable == 'cyc':
                     continue
-                eddies.obs[VAR_DESCR_inv[variable]] = h_nc.variables[variable][:]
+                eddies.obs[VAR_DESCR_inv[variable]
+                           ] = h_nc.variables[variable][:]
             eddies.sign_type = h_nc.variables['cyc'][0]
         return eddies
 
     def tracking(self, other):
-        """Track obs between from self to other
+        """Track obs between self and other
         """
-        dist = self.distance(other)
-        i_self, i_other = where(dist < 20.)
+        cost = self.distance(other)
+        # Links available which respect a maximal cost
+        mask_accept_cost = cost < 75
+        cost = ma.array(cost, mask=-mask_accept_cost, dtype='i2')
+
+        # Count number of link by self obs and other obs
+        self_links = mask_accept_cost.sum(axis=1)
+        other_links = mask_accept_cost.sum(axis=0)
+        max_links = max(self_links.max(), other_links.max())
+        if max_links > 5:
+            logging.warning('One observation have %d links', max_links)
+
+        # If some obs have multiple link, we keep only one link by eddy
+        eddies_separation = 1 < self_links
+        eddies_merge = 1 < other_links
+        test = eddies_separation.any() or eddies_merge.any()
+        if test:
+            # We extract matrix which contains concflict
+            obs_linking_to_self = mask_accept_cost[eddies_separation
+                                                   ].any(axis=0)
+            obs_linking_to_other = mask_accept_cost[:, eddies_merge
+                                                    ].any(axis=1)
+            i_self_keep = where(obs_linking_to_other + eddies_separation)[0]
+            i_other_keep = where(obs_linking_to_self + eddies_merge)[0]
+
+            # Cost to resolve conflict
+            cost_reduce = cost[i_self_keep][:, i_other_keep]
+            shape = cost_reduce.shape
+            logging.debug('Shape conflict matrix : %s', shape)
+
+            matrix_size = shape[0] * shape[1]
+            if (matrix_size) >= 20000:
+                logging.warning('High number of conflict : %d (matrix_size)',
+                                matrix_size)
+
+            links_resolve = 0
+            while False in cost_reduce.mask:
+                i_min_value = cost_reduce.argmin()
+                i, j = i_min_value / shape[1], i_min_value % shape[1]
+                # Set to False all link
+                mask_accept_cost[i_self_keep[i]] = False
+                mask_accept_cost[:, i_other_keep[j]] = False
+                cost_reduce.mask[i] = True
+                cost_reduce.mask[:, j] = True
+                # we active only this link
+                mask_accept_cost[i_self_keep[i], i_other_keep[j]] = True
+                links_resolve += 1
+            logging.debug('%d links resolve', links_resolve)
+
+        i_self, i_other = where(mask_accept_cost)
 
         logging.debug('%d matched with previous', i_self.shape[0])
+
+        # Check
+        if unique(i_other).shape[0] != i_other.shape[0]:
+            raise Exception()
+        if unique(i_self).shape[0] != i_self.shape[0]:
+            raise Exception()
         return i_self, i_other
 
 
 class VirtualEddiesObservations(EddiesObservations):
-    
+    """Class to work with virtual obs
+    """
+
     @property
     def elements(self):
         elements = super(VirtualEddiesObservations, self).elements
@@ -211,7 +285,9 @@ def elements(self):
 
 
 class TrackEddiesObservations(EddiesObservations):
-    
+    """Class to practice Tracking on observations
+    """
+
     def filled_by_interpolation(self, mask):
         """Filled selected values by interpolation
         """
@@ -228,42 +304,47 @@ def filled_by_interpolation(self, mask):
                                          self.obs[var][-mask])
 
     def extract_longer_eddies(self, nb_min, nb_obs, compress_id=True):
-        m = nb_obs >= nb_min
-        nb_obs_select = m.sum()
+        """Select eddies which are longer than nb_min
+        """
+        mask = nb_obs >= nb_min
+        nb_obs_select = mask.sum()
         logging.info('Selection of %d observations', nb_obs_select)
         eddies = TrackEddiesObservations(size=nb_obs_select)
         eddies.sign_type = self.sign_type
         for var, _ in eddies.obs.dtype.descr:
-            eddies.obs[var] = self.obs[var][m]
+            eddies.obs[var] = self.obs[var][mask]
         if compress_id:
             list_id = unique(eddies.obs['track'])
             list_id.sort()
             id_translate = arange(list_id.max() + 1)
             id_translate[list_id] = arange(len(list_id)) + 1
             eddies.obs['track'] = id_translate[eddies.obs['track']]
         return eddies
-    
+
     @property
     def elements(self):
         elements = super(TrackEddiesObservations, self).elements
         elements.extend(['track', 'n', 'virtual'])
         return elements
-    
-    def create_variable(self, handler_nc, kwargs_variable,
+
+    @staticmethod
+    def create_variable(handler_nc, kwargs_variable,
                         attr_variable, data, scale_factor=None):
+        """Create variable
+        """
         var = handler_nc.createVariable(
             zlib=True,
             complevel=1,
             **kwargs_variable)
         for attr, attr_value in attr_variable.iteritems():
             var.setncattr(attr, attr_value)
-            
+
         var[:] = data
-        
+
         #~ var.set_auto_maskandscale(False)
         if scale_factor is not None:
             var.scale_factor = scale_factor
-            
+
         try:
             var.setncattr('min', var[:].min())
             var.setncattr('max', var[:].max())
@@ -291,7 +372,10 @@ def write_netcdf(self):
                          dimensions=VAR_DESCR[name]['nc_dims']),
                     VAR_DESCR[name]['nc_attr'],
                     self.observations[name],
-                    scale_factor=None if 'scale_factor' not in VAR_DESCR[name] else VAR_DESCR[name]['scale_factor'])
+                    scale_factor=None
+                        if 'scale_factor' not in VAR_DESCR[name] else
+                        VAR_DESCR[name]['scale_factor']
+                    )
 
             # Add cyclonic information
             self.create_variable(
@@ -305,7 +389,12 @@ def write_netcdf(self):
             self.set_global_attr_netcdf(h_nc)
 
     def set_global_attr_netcdf(self, h_nc):
-        h_nc.title = 'Cyclonic' if self.sign_type == -1 else 'Anticyclonic' + ' eddy tracks'
+        """Set global attr
+        """
+        if self.sign_type == -1:
+            h_nc.title = 'Cyclonic'
+        else:
+            h_nc.title = 'Anticyclonic' + ' eddy tracks'
         #~ h_nc.grid_filename = self.grd.grid_filename
         #~ h_nc.grid_date = str(self.grd.grid_date)
         #~ h_nc.product = self.product
@@ -324,7 +413,7 @@ def set_global_attr_netcdf(self, h_nc):
         #~ h_nc.evolve_amp_max = self.evolve_amp_max
         #~ h_nc.evolve_area_min = self.evolve_area_min
         #~ h_nc.evolve_area_max = self.evolve_area_max
-#~ 
+
         #~ h_nc.llcrnrlon = self.grd.lonmin
         #~ h_nc.urcrnrlon = self.grd.lonmax
         #~ h_nc.llcrnrlat = self.grd.latmin