- Cleanup

- Save always, short untracked
- Add propagation function for virtual

Author: AntSimi

src/py_eddy_tracker/featured_tracking/old_tracker_reference.py

@@ -18,19 +18,20 @@ def mask_function(self, other):
         """We mask link with ellips and ratio
         """
         # Compute Parameter of ellips
+        minor, major = 1.05, 1.5
         y = self.basic_formula_ellips_major_axis(
             self.obs['lat'],
             degrees=True,
-            c0=1.05,
-            cmin=1.05,
-            cmax=1.5,
+            c0=minor,
+            cmin=minor,
+            cmax=major,
             lat1=23,
             lat2=5,
         )
         # mask from ellips
         mask = self.shifted_ellipsoid_degrees_mask(
             other,
-            minor=1.5,
+            minor=minor, # Minor can be bigger than major??
             major=y)
 
         # We check ratio

src/py_eddy_tracker/grid/__init__.py

@@ -50,10 +50,10 @@ def browse_dataset_in(data_dir, files_model, date_regexp, date_model,
 
     filenames = bytes_(glob(full_path))
     dataset_list = empty(len(filenames),
-                         dtype=[('filename', 'S256'),
+                         dtype=[('filename', 'S500'),
                                 ('date', 'datetime64[D]'),
                                 ])
-    dataset_list['filename'] = bytes_(glob(full_path))
+    dataset_list['filename'] = filenames
 
     logging.info('%s grids available', dataset_list.shape[0])
     mode_attrs = False

src/py_eddy_tracker/observations.py

@@ -27,9 +27,9 @@
 ===========================================================================
 
 """
-from numpy import zeros, empty, nan, arange, interp, where, unique, \
+from numpy import zeros, empty, nan, arange, where, unique, \
     ma, concatenate, cos, radians, isnan, ones, ndarray, meshgrid, \
-    bincount, bool_, array, interp, int_, int32, round, maximum, floor
+    array, interp, int_, int32, round, maximum, floor
 from scipy.interpolate  import interp1d
 from netCDF4 import Dataset
 from py_eddy_tracker.tools import distance_matrix, distance_vector
@@ -38,7 +38,6 @@
 from . import VAR_DESCR, VAR_DESCR_inv
 import logging
 from datetime import datetime
-from scipy.interpolate  import RectBivariateSpline
 
 
 class GridDataset(object):
@@ -185,6 +184,10 @@ def __init__(self, size=0, track_extra_variables=None,
         self.active = True
         self.sign_type = None
 
+    @property
+    def sign_legend(self):
+        return 'Cyclonic' if self.sign_type == -1 else 'Anticyclonic'
+
     @property
     def shape(self):
         return self.observations.shape
@@ -222,7 +225,7 @@ def elements(self):
 
         if len(self.track_extra_variables):
             elements += self.track_extra_variables
-        return elements
+        return list(set(elements))
 
     def coherence(self, other):
         """Check coherence between two dataset
@@ -341,6 +344,8 @@ def index(self, index):
     @classmethod
     def load_from_netcdf(cls, filename):
         array_dim = 'NbSample'
+        if not isinstance(filename, str):
+            filename = filename.astype(str)
         with Dataset(filename) as h_nc:
             nb_obs = len(h_nc.dimensions['Nobs'])
             kwargs = dict()
@@ -386,6 +391,53 @@ def from_netcdf(cls, handler):
             eddies.obs[VAR_DESCR_inv[variable]] = handler.variables[variable][:]
         return eddies
 
+    @staticmethod
+    def propagate(previous_obs, current_obs, obs_to_extend, dead_track, nb_next, model):
+        """
+        Filled virtual obs (C)
+        Args:
+            previous_obs: previous obs from current (A)
+            current_obs: previous obs from virtual (B)
+            obs_to_extend:
+            dead_track:
+            nb_next:
+            model: 
+
+        Returns:
+            New position C = B + AB
+        """
+        next_obs = VirtualEddiesObservations(
+            size=nb_next,
+            track_extra_variables=model.track_extra_variables,
+            track_array_variables=model.track_array_variables,
+            array_variables=model.array_variables)
+        nb_dead = len(previous_obs)
+        nb_virtual_extend = nb_next - nb_dead
+
+        for key in model.elements:
+            if key in ['lon', 'lat', 'time'] or 'contour_' in key:
+                continue
+            next_obs[key][:nb_dead] = current_obs[key]
+        next_obs['dlon'][:nb_dead] = current_obs['lon'] - previous_obs['lon']
+        next_obs['dlat'][:nb_dead] = current_obs['lat'] - previous_obs['lat']
+        next_obs['lon'][:nb_dead] = current_obs['lon'] + next_obs['dlon'][:nb_dead]
+        next_obs['lat'][:nb_dead] = current_obs['lat'] + next_obs['dlat'][:nb_dead]
+        # Id which are extended
+        next_obs['track'][:nb_dead] = dead_track
+        # Add previous virtual
+        if nb_virtual_extend > 0:
+            for key in next_obs.elements:
+                if key in ['lon', 'lat', 'time', 'track', 'segment_size'] or 'contour_' in key:
+                    continue
+                next_obs[key][nb_dead:] = obs_to_extend[key]
+            next_obs['lon'][nb_dead:] = obs_to_extend['lon'] + obs_to_extend['dlon']
+            next_obs['lat'][nb_dead:] = obs_to_extend['lat'] + obs_to_extend['dlat']
+            next_obs['track'][nb_dead:] = obs_to_extend['track']
+            next_obs['segment_size'][nb_dead:] = obs_to_extend['segment_size']
+        # Count
+        next_obs['segment_size'][:] += 1
+        return next_obs
+
     @staticmethod
     def cost_function2(records_in, records_out, distance):
         nb_records = records_in.shape[0]
@@ -715,8 +767,7 @@ def write_netcdf(self, path='./', filename='%(path)s/%(sign_type)s.nc'):
         """Write a netcdf with eddy obs
         """
         eddy_size = len(self.observations)
-        sign_type = 'Cyclonic' if self.sign_type == -1 else 'Anticyclonic'
-        filename = filename % dict(path=path, sign_type=sign_type, prod_time=datetime.now().strftime('%Y%m%d'))
+        filename = filename % dict(path=path, sign_type=self.sign_legend, prod_time=datetime.now().strftime('%Y%m%d'))
         logging.info('Store in %s', filename)
         with Dataset(filename, 'w', format='NETCDF4') as h_nc:
             logging.info('Create file %s', filename)
@@ -763,7 +814,7 @@ class VirtualEddiesObservations(EddiesObservations):
     def elements(self):
         elements = super(VirtualEddiesObservations, self).elements
         elements.extend(['track', 'segment_size', 'dlon', 'dlat'])
-        return elements
+        return list(set(elements))
 
 
 class TrackEddiesObservations(EddiesObservations):
@@ -783,8 +834,16 @@ def filled_by_interpolation(self, mask):
             var = field[0]
             if var in ['n', 'virtual', 'track'] or var in self.array_variables:
                 continue
-            self.obs[var][mask] = interp(index[mask], index[~mask],
-                                         self.obs[var][~mask])
+            # to normalize longitude before interpolation
+            if var== 'lon':
+                lon = self.obs[var]
+                first = where(self.obs['n'] == 0)[0]
+                nb_obs = empty(first.shape, dtype='u4')
+                nb_obs[:-1] = first[1:] - first[:-1]
+                nb_obs[-1] = lon.shape[0] - first[-1]
+                lon0 = (lon[first] - 180).repeat(nb_obs)
+                self.obs[var] = (lon - lon0) % 360 + lon0
+            self.obs[var][mask] = interp(index[mask], index[~mask], self.obs[var][~mask])
 
     def extract_longer_eddies(self, nb_min, nb_obs, compress_id=True):
         """Select eddies which are longer than nb_min

src/py_eddy_tracker/tracking.py

@@ -260,12 +260,6 @@ def recense_dead_id_to_extend(self):
 
         # Save previous state to count virtual obs
         self.previous_virtual_obs = self.virtual_obs
-        # Creation of an virtual step for dead one
-        self.virtual_obs = VirtualEddiesObservations(
-            size=nb_dead + nb_virtual_extend,
-            track_extra_variables=self.previous_obs.track_extra_variables,
-            track_array_variables=self.previous_obs.track_array_variables,
-            array_variables=self.previous_obs.array_variables)
 
         # Find mask/index on previous correspondance to extrapolate
         # position
@@ -275,33 +269,13 @@ def recense_dead_id_to_extend(self):
         # Selection of observations on N-2 and N-1
         obs_a = self.previous2_obs.obs[self[-2][i_dead_id]['in']]
         obs_b = self.previous_obs.obs[self[-2][i_dead_id]['out']]
-        # Position N-2 : A
-        # Position N-1 : B
-        # Virtual Position : C
-        # New position C = B + AB
-        for key in self.previous_obs.elements:
-            if key in ['lon', 'lat', 'time'] or 'contour_' in key:
-                continue
-            self.virtual_obs[key][:nb_dead] = obs_b[key]
-        self.virtual_obs['dlon'][:nb_dead] = obs_b['lon'] - obs_a['lon']
-        self.virtual_obs['dlat'][:nb_dead] = obs_b['lat'] - obs_a['lat']
-        self.virtual_obs['lon'][:nb_dead] = obs_b['lon'] + self.virtual_obs['dlon'][:nb_dead]
-        self.virtual_obs['lat'][:nb_dead] = obs_b['lat'] + self.virtual_obs['dlat'][:nb_dead]
-        # Id which are extended
-        self.virtual_obs['track'][:nb_dead] = dead_id
-        # Add previous virtual
-        if nb_virtual_extend > 0:
-            obs_to_extend = self.previous_virtual_obs.obs[i_virtual_dead_id][alive_virtual_obs]
-            for key in self.virtual_obs.elements:
-                if key in ['lon', 'lat', 'time', 'track', 'segment_size'] or 'contour_' in key:
-                    continue
-                self.virtual_obs[key][nb_dead:] = obs_to_extend[key]
-            self.virtual_obs['lon'][nb_dead:] = obs_to_extend['lon'] + obs_to_extend['dlon']
-            self.virtual_obs['lat'][nb_dead:] = obs_to_extend['lat'] + obs_to_extend['dlat']
-            self.virtual_obs['track'][nb_dead:] = obs_to_extend['track']
-            self.virtual_obs['segment_size'][nb_dead:] = obs_to_extend['segment_size']
-        # Count
-        self.virtual_obs['segment_size'][:] += 1
+
+        self.virtual_obs = self.previous_obs.propagate(
+            obs_a, obs_b,
+            self.previous_virtual_obs.obs[i_virtual_dead_id][alive_virtual_obs] if nb_virtual_extend > 0 else None,
+            dead_track=dead_id,
+            nb_next=nb_dead + nb_virtual_extend,
+            model=self.previous_obs)
 
     def load_state(self):
         # If we have a previous file of correspondance, we will replay only recent part
@@ -616,20 +590,18 @@ def get_unused_data(self):
         has_virtual = 'virtual' in self[0].dtype.names
         for i, filename in enumerate(self.datasets):
             last_dataset = i == (nb_dataset - 1)
-            first_dataset = i == 0
-            if has_virtual:
-                if not last_dataset:
-                    m_in = ~self[i]['virtual']
-                if not first_dataset:
-                    m_out = ~self[i - 1]['virtual']
+            if has_virtual and not last_dataset:
+                m_in = ~self[i]['virtual']
             else:
-                m_in, m_out = slice(None), slice(None)
+                m_in = slice(None)
             if i == 0:
-                eddies_used = self[i]['in'][m_in]
+                eddies_used = self[i]['in']
             elif last_dataset:
-                eddies_used = self[i - 1]['out'][m_out]
+                eddies_used = self[i - 1]['out']
             else:
-                eddies_used = unique(concatenate((self[i - 1]['out'][m_out], self[i]['in'][m_in])))
+                eddies_used = unique(concatenate((self[i - 1]['out'], self[i]['in'][m_in])))
+            if not isinstance(filename, str):
+                filename = filename.astype(str)
             with Dataset(filename) as h:
                 nb_obs_day = len(h.dimensions['Nobs'])
             m = ones(nb_obs_day, dtype='bool')
@@ -645,7 +617,7 @@ def get_unused_data(self):
         j = 0
         for i, dataset in enumerate(self.datasets):
             current_obs = self.class_method.load_from_netcdf(dataset)
-            if i ==0:
+            if i == 0:
                 eddies.sign_type = current_obs.sign_type
             unused_obs = current_obs.observations[list_mask[i]]
             nb = unused_obs.shape[0]

src/scripts/EddyFinalTracking

@@ -8,7 +8,6 @@ from py_eddy_tracker.tracking import Correspondances
 from os.path import exists
 from os import mkdir
 import logging
-from numpy import unique
 import datetime as dt
 
 
@@ -24,9 +23,6 @@ def usage():
     parser.add_argument('--path_out',
                         default='./',
                         help='Path, where to write file')
-    parser.add_argument('--filename_model',
-                        default='%(path)s/%(sign_type)s.nc',
-                        help='Model to create file')
     parser.add_argument('--nb_obs_min',
                         type=int,
                         default=28,
@@ -40,7 +36,8 @@ if __name__ == '__main__':
     # Create output directory
     if not exists(CONFIG.path_out):
         mkdir(CONFIG.path_out)
-
+    SAVE_DIR = CONFIG.path_out
+    NB_OBS_MIN = CONFIG.nb_obs_min
     START_TIME = dt.datetime.now()
 
     CORRESPONDANCES = Correspondances.load(CONFIG.nc_file)
@@ -51,18 +48,30 @@ if __name__ == '__main__':
 
     logging.info('The longest tracks have %d observations', CORRESPONDANCES.nb_obs_by_tracks.max())
     logging.info('The mean length is %d observations before filtering', CORRESPONDANCES.nb_obs_by_tracks.mean())
-    CORRESPONDANCES.longer_than(size_min=CONFIG.nb_obs_min)
+
+    CORRESPONDANCES.get_unused_data().write_netcdf(path=SAVE_DIR, filename='%(path)s/%(sign_type)s_untracked.nc')
+
+    SHORT_CORRESPONDANCES = CORRESPONDANCES._copy()
+    SHORT_CORRESPONDANCES.shorter_than(size_max=NB_OBS_MIN)
+
+    CORRESPONDANCES.longer_than(size_min=NB_OBS_MIN)
+
     FINAL_EDDIES = CORRESPONDANCES.merge()
+    SHORT_TRACK = SHORT_CORRESPONDANCES.merge()
 
     # We flag obs
     if CORRESPONDANCES.virtual:
         FINAL_EDDIES['virtual'][:] = FINAL_EDDIES['time'] == 0
         FINAL_EDDIES.filled_by_interpolation(FINAL_EDDIES['virtual'] == 1)
+        SHORT_TRACK['virtual'][:] = SHORT_TRACK['time'] == 0
+        SHORT_TRACK.filled_by_interpolation(SHORT_TRACK['virtual'] == 1)
 
+    # Total running time
     FULL_TIME = dt.datetime.now() - START_TIME
     logging.info('Duration : %s', FULL_TIME)
 
-    logging.info('%d tracks will be saved', len(unique(FINAL_EDDIES['track'])))
+    logging.info('Longer track saved have %d obs', CORRESPONDANCES.nb_obs_by_tracks.max())
     logging.info('The mean length is %d observations after filtering', CORRESPONDANCES.nb_obs_by_tracks.mean())
 
-    FINAL_EDDIES.write_netcdf(path=CONFIG.path_out, filename=CONFIG.filename_model)
+    FINAL_EDDIES.write_netcdf(path=SAVE_DIR)
+    SHORT_TRACK.write_netcdf(filename='%(path)s/%(sign_type)s_track_too_short.nc', path=SAVE_DIR)

src/scripts/EddyIdentification

@@ -22,7 +22,6 @@ Email: [email protected]
 """
 
 import logging
-from glob import glob
 from yaml import load as yaml_load
 from datetime import datetime
 import datetime as dt
@@ -31,7 +30,7 @@ from matplotlib.dates import date2num, num2julian
 from matplotlib.figure import Figure
 from os.path import exists
 from os import mkdir
-from numpy import array, arange, atleast_1d, unique, round, int, interp, pi, \
+from numpy import arange, atleast_1d, round, int, interp, pi, \
     float_, errstate, ma
 
 from py_eddy_tracker import EddyParser