Extraction of association cost

Author: AntSimi

src/py_eddy_tracker/__init__.py

@@ -167,6 +167,16 @@ def parse_args(self, *args, **kwargs):
             description='Virtual observation: 0 for real',
             )
         ),
+    cost_association=dict(
+        attr_name=None,
+        nc_name='cost_association',
+        nc_type='float32',
+        nc_dims=('Nobs',),
+        nc_attr=dict(
+            long_name='cost_value_to_associate_with_next_observation',
+            description='Cost value to associate with the next observation',
+            )
+        ),
     lon=dict(
         attr_name='lon',
         compute_type='float64',

src/py_eddy_tracker/observations.py

@@ -164,7 +164,7 @@ class EddiesObservations(object):
 
     ELEMENTS = ['lon', 'lat', 'radius_s', 'radius_e', 'amplitude', 'speed_radius', 'time', 'eke',
                 'shape_error_e', 'shape_error_s', 'nb_contour_selected',
-                'height_max_speed_contour', 'height_external_contour', 'height_inner_contour']
+                'height_max_speed_contour', 'height_external_contour', 'height_inner_contour', 'cost_association']
 
     def __init__(self, size=0, track_extra_variables=None,
                  track_array_variables=0, array_variables=None):
@@ -452,22 +452,34 @@ def propagate(previous_obs, current_obs, obs_to_extend, dead_track, nb_next, mod
         return next_obs
 
     @staticmethod
-    def cost_function_common_area(records_in, records_out, distance):
+    def cost_function_common_area(records_in, records_out, distance, intern=False):
+        x_name, y_name = ('contour_lon_s', 'contour_lat_s') if intern else ('contour_lon_e', 'contour_lat_e')
         nb_records = records_in.shape[0]
         costs = ma.empty(nb_records,dtype='f4')
+        tolerance = 0.025
         for i_record in range(nb_records):
             poly_in = Polygon(
                 concatenate((
-                    (records_in[i_record]['contour_lon_e'],),
-                    (records_in[i_record]['contour_lat_e'],))
+                    (records_in[i_record][x_name],),
+                    (records_in[i_record][y_name],))
                     ).T
                 )
             poly_out = Polygon(
                 concatenate((
-                    (records_out[i_record]['contour_lon_e'],),
-                    (records_out[i_record]['contour_lat_e'],))
+                    (records_out[i_record][x_name],),
+                    (records_out[i_record][y_name],))
                     ).T
                 )
+            if not poly_in.is_valid:
+                poly_in = poly_in.simplify(tolerance=tolerance)
+                if not poly_in.is_valid:
+                    logging.warning('Need to simplify polygon in for a second time')
+                    poly_in = poly_in.buffer(0)
+            if not poly_out.is_valid:
+                poly_out = poly_out.simplify(tolerance=tolerance)
+                if not poly_out.is_valid:
+                    logging.warning('Need to simplify polygon out for a second time')
+                    poly_out = poly_out.buffer(0)
             try:
                 costs[i_record] = 1 - poly_in.intersection(poly_out).area / poly_in.area
             except TopologicalError:
@@ -725,7 +737,7 @@ def tracking(self, other):
 
         logging.debug('%d matched with previous', i_self.shape[0])
 
-        return i_self, i_other
+        return i_self, i_other, cost_mat[i_self, i_other]
 
     def to_netcdf(self, handler):
         eddy_size = len(self)
@@ -858,7 +870,7 @@ def filled_by_interpolation(self, mask):
 
         for field in self.obs.dtype.descr:
             var = field[0]
-            if var in ['n', 'virtual', 'track'] or var in self.array_variables:
+            if var in ['n', 'virtual', 'track', 'cost_association'] or var in self.array_variables:
                 continue
             # to normalize longitude before interpolation
             if var== 'lon':

src/py_eddy_tracker/tracking.py

@@ -30,8 +30,9 @@
 from matplotlib.dates import julian2num, num2date
 
 from py_eddy_tracker.observations import EddiesObservations, VirtualEddiesObservations, TrackEddiesObservations
-from numpy import bool_, array, arange, ones, setdiff1d, zeros, uint16, where, empty, isin, unique, concatenate, ma
-from netCDF4 import Dataset
+from numpy import bool_, array, arange, ones, setdiff1d, zeros, uint16, where, empty, isin, unique, concatenate, \
+    ma
+from netCDF4 import Dataset, default_fillvals
 import logging
 import platform
 
@@ -55,7 +56,9 @@ def __init__(self, datasets, virtual=0, class_method=None, previous_correspondan
         # Correspondance dtype
         self.correspondance_dtype = [('in', 'u2'),
                                      ('out', 'u2'),
-                                     ('id', self.ID_DTYPE)]
+                                     ('id', self.ID_DTYPE),
+                                     ('cost_value', 'f4')
+                                     ]
         if class_method is None:
             self.class_method = EddiesObservations
         else:
@@ -168,7 +171,7 @@ def merge_correspondance(self, other):
         # We set new id available
         self.current_id = translate[-1] + 1
 
-    def store_correspondance(self, i_previous, i_current, nb_real_obs):
+    def store_correspondance(self, i_previous, i_current, nb_real_obs, association_cost):
         """Storing correspondance in an array
         """
         # Create array to store correspondance data
@@ -177,6 +180,7 @@ def store_correspondance(self, i_previous, i_current, nb_real_obs):
             correspondance['virtual_length'][:] = 255
         # index from current_obs
         correspondance['out'] = i_current
+        correspondance['cost_value'] = association_cost
 
         if self.virtual:
             # if index in previous dataset is bigger than real obs number
@@ -314,10 +318,10 @@ def track(self):
             if flg_virtual:
                 logging.debug('%d virtual obs will be add to previous', len(self.virtual_obs))
                 self.previous_obs = self.previous_obs.merge(self.virtual_obs)
-            i_previous, i_current = self.previous_obs.tracking(self.current_obs)
+            i_previous, i_current, association_cost = self.previous_obs.tracking(self.current_obs)
 
             # return true if the first time (previous2obs is none)
-            if self.store_correspondance(i_previous, i_current, nb_real_obs):
+            if self.store_correspondance(i_previous, i_current, nb_real_obs, association_cost):
                 continue
 
             self.recense_dead_id_to_extend()
@@ -376,7 +380,11 @@ def save(self, filename, dict_completion=None):
                 for name, _ in self.correspondance_dtype:
                     datas[name][i, :nb_elt] = correspondance[name]
             for name, data in datas.items():
-                h_nc.variables[name][:] = data
+                h_v = h_nc.variables[name]
+                h_v[:] = data
+                if 'File' not in name:
+                    h_v.min = h_v[:].min()
+                    h_v.max = h_v[:].max()
 
             h_nc.virtual_use = str(self.virtual)
             h_nc.virtual_max_segment = self.nb_virtual
@@ -512,6 +520,8 @@ def merge(self, until=-1):
             track_array_variables=self.current_obs.track_array_variables,
             array_variables=self.current_obs.array_variables)
 
+        # All the value put at nan, necessary only for all end of track
+        eddies['cost_association'][:] = default_fillvals['f4']
         # Calculate the index in each tracks, we compute in u4 and translate
         # in u2 (which are limited to 65535)
         logging.debug('Compute global index array (N)')
@@ -547,6 +557,8 @@ def merge(self, until=-1):
                 # Copy all variable
                 for field in fields:
                     var = field[0]
+                    if var == 'cost_association':
+                        continue
                     eddies[var][index_final[m_first_obs]] = self.previous_obs[var][index_in]
                 # Increment
                 self.i_current_by_tracks[i_id[m_first_obs]] += 1
@@ -570,7 +582,10 @@ def merge(self, until=-1):
             # Copy all variable
             for field in fields:
                 var = field[0]
-                eddies[var][index_final] = self.current_obs[var][index_current]
+                if var == 'cost_association':
+                    eddies[var][index_final-1] = self[i]['cost_value']
+                else:
+                    eddies[var][index_final] = self.current_obs[var][index_current]
 
             # Add increment for each index used
             self.i_current_by_tracks[i_id] += 1

src/scripts/EddyTracking

@@ -104,7 +104,6 @@ if __name__ == '__main__':
     DICT_COMPLETION = dict(date_start=DATE_START, date_stop=DATE_STOP, date_prod=START_TIME,
                            path=SAVE_DIR, sign_type=CORRESPONDANCES.current_obs.sign_legend)
 
-
     CORRESPONDANCES.save(CORRESPONDANCES_OUT, DICT_COMPLETION)
     if SAVE_STOP:
         exit()