New way to set tracking

Author: adelepoulle

share/tracking.yaml

@@ -4,21 +4,18 @@ DIAGNOSTIC_TYPE: 'SLA'
 
 PATHS:
   # Files produces with EddyIdentification
-  FILES_PATTERN: /homelocal/adelepoulle/20151110_eddy_track/Anticyclonic_2014*.nc
+  FILES_PATTERN: id_track/acyc_eddy_stat_19930*.nc
   # Path and filename of Chelton et al (1998) Rossby radius data
   # Obtain file from:
   # http://www-po.coas.oregonstate.edu/research/po/research/rossby_radius/
   RW_PATH: 'rossrad.dat'
   # Path for saving of outputs
-  SAVE_DIR: 'test/'
+  SAVE_DIR: 'tracks_acyc_test/'
 
 # Minimum number of observations to store eddy
-TRACK_DURATION_MIN: 15
-VIRTUAL_LEGNTH_MAX: 2
-  
-SEPARATION_METHOD: 'ellipse'
-#SEPARATION_METHOD: 'sum_radii'
+TRACK_DURATION_MIN: 4
+VIRTUAL_LEGNTH_MAX: 0
 
-
-# Typical parameters
-DIST0: 25000. # m separation distance after ~7 days (see CSS11 fig 22)
+CLASS:
+    MODULE: py_eddy_tracker.featured_tracking.old_tracker_reference
+    CLASS: CheltonTracker

src/py_eddy_tracker/featured_tracking/__init__.py

@@ -0,0 +1,3 @@
+# -*- coding: utf-8 -*-
+"""
+"""

src/py_eddy_tracker/featured_tracking/mask_old.nc

@@ -0,0 +1,82 @@
+from ..observations import EddiesObservations as Model
+from ..observations import GridDataset
+from numpy import where, meshgrid, bincount, ones, unique, bool_, arange
+import logging
+from os import path
+
+
+class CheltonTracker(Model):
+    GROUND = GridDataset(path.join(path.dirname(__file__), 'mask_old.nc'), 'lon', 'lat')
+
+    @staticmethod
+    def cost_function(records_in, records_out, distance):
+        """We minimize on distance between two obs
+        """
+        return distance
+
+    def mask_function(self, other):
+        """We mask link with ellips and ratio
+        """
+        # Compute Parameter of ellips
+        y = self.basic_formula_ellips_major_axis(
+            self.obs['lat'],
+            degrees=True)
+        # mask from ellips
+        mask = self.shifted_ellipsoid_degrees_mask(
+            other,
+            minor=1.5,
+            major=y)
+
+        # We check ratio
+        other_amplitude, self_amplitude = meshgrid(
+            other.obs['amplitude'], self.obs['amplitude'])
+        other_radius, self_radius = meshgrid(
+            other.obs['radius_e'], self.obs['radius_e'])
+        ratio_amplitude = other_amplitude / self_amplitude
+        ratio_radius = other_radius / self_radius
+        mask *= \
+            (ratio_amplitude < 2.5) * \
+            (ratio_amplitude > 1 / 2.5) * \
+            (ratio_radius < 2.5) * \
+            (ratio_radius > 1 / 2.5)
+        indexs_closest = where(mask)
+        mask[indexs_closest] = self.across_ground(self.obs[indexs_closest[0]], other.obs[indexs_closest[1]])
+        return mask
+
+    @classmethod
+    def across_ground(cls, record0, record1):
+        i, j, d_pix = cls.GROUND.compute_pixel_path(
+            x0=record0['lon'],
+            y0=record0['lat'],
+            x1=record1['lon'],
+            y1=record1['lat'],
+            )
+
+        data = cls.GROUND.grid('mask')[i, j]
+        i_ground = unique(arange(len(record0)).repeat(d_pix + 1)[data==1])
+        mask = ones(record1.shape, dtype='bool')
+        mask[i_ground] = False
+        return mask
+
+    def solve_function(self, cost_matrix):
+        """Give the best link for each self obs
+        """
+        return where(self.solve_first(cost_matrix, multiple_link=True))
+
+    def post_process_link(self, other, i_self, i_other):
+        """When two self obs use the same other obs, we keep the self obs
+        with amplitude max
+        """
+        if unique(i_other).shape[0] != i_other.shape[0]:
+            nb_link = bincount(i_other)
+            mask = ones(i_self.shape, dtype=bool_)
+            for i in where(nb_link > 1)[0]:
+                m = i == i_other
+                multiple_in = i_self[m]
+                i_keep = self.obs['amplitude'][multiple_in].argmax()
+                m[where(m)[0][i_keep]] = False
+                mask[m] = False
+
+            i_self = i_self[mask]
+            i_other = i_other[mask]
+        return i_self, i_other

src/py_eddy_tracker/featured_tracking/old_tracker_reference.py

@@ -99,10 +99,12 @@ def __init__(self, aviso_file, the_domain,
     def set_filename(self, file_name):
         self.grid_filename = file_name
 
-    def get_aviso_data(self, aviso_file):
+    def get_aviso_data(self, aviso_file, dimensions=None):
         """
         Read nc data from AVISO file
         """
+        if dimensions is None:
+            dimensions = dict()
         self.grid_filename = aviso_file
         units = self.read_nc_att(self.grid_name, 'units')
         if units not in self.KNOWN_UNITS:
@@ -123,7 +125,8 @@ def get_aviso_data(self, aviso_file):
             elif grid_dim == lon_dim:
                 i_list.append(self.view_pad[1])
             else:
-                i_list.append(0)
+                i_list.append(
+                    dimensions[grid_dim] if grid_dim in dimensions.keys() else 0)
 
         zeta = self.read_nc(self.grid_name, indices=i_list)
         if transpose: