- Add period to cube

- Add some methods for display
- Speed up overlap

Author: AntSimi

README.md

@@ -1,11 +1,16 @@
 [![PyPI version](https://badge.fury.io/py/pyEddyTracker.svg)](https://badge.fury.io/py/pyEddyTracker)
+[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.6333988.svg)](https://doi.org/10.5281/zenodo.6333988)
 [![Documentation Status](https://readthedocs.org/projects/py-eddy-tracker/badge/?version=stable)](https://py-eddy-tracker.readthedocs.io/en/stable/?badge=stable)
 [![Gitter](https://badges.gitter.im/py-eddy-tracker/community.svg)](https://gitter.im/py-eddy-tracker/community?utm_source=badge&utm_medium=badge&utm_campaign=pr-badge) 
 [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/AntSimi/py-eddy-tracker/master?urlpath=lab/tree/notebooks/python_module/)
 [![pytest](https://github.com/AntSimi/py-eddy-tracker/actions/workflows/python-app.yml/badge.svg)](https://github.com/AntSimi/py-eddy-tracker/actions/workflows/python-app.yml)
 
 # README #
 
+### How to cite code? ###
+
+Zenodo provide DOI for each tagged version, [all DOI are available here](https://doi.org/10.5281/zenodo.6333988)
+
 ### Method ###
 
 Method was described in :

src/py_eddy_tracker/appli/network.py

@@ -128,6 +128,9 @@ def subset_network():
         action="store_true",
         help="Remove trash (network id == 0)",
     )
+    parser.add_argument(
+        "-i", "--ids", nargs="+", type=int, help="List of network which will be extract"
+    )
     parser.add_argument(
         "-p",
         "--period",
@@ -138,6 +141,8 @@ def subset_network():
     )
     args = parser.parse_args()
     n = NetworkObservations.load_file(args.input, raw_data=True)
+    if args.ids is not None:
+        n = n.networks(args.ids)
     if args.length is not None:
         n = n.longer_than(*args.length)
     if args.remove_dead_end is not None:

src/py_eddy_tracker/dataset/grid.py

@@ -2,45 +2,20 @@
 """
 Class to load and manipulate RegularGrid and UnRegularGrid
 """
-from datetime import datetime
 import logging
+from datetime import datetime
 
 from cv2 import filter2D
 from matplotlib.path import Path as BasePath
 from netCDF4 import Dataset
-from numba import njit, prange, types as numba_types
-from numpy import (
-    arange,
-    array,
-    ceil,
-    concatenate,
-    cos,
-    deg2rad,
-    empty,
-    errstate,
-    exp,
-    float_,
-    floor,
-    histogram2d,
-    int_,
-    interp,
-    isnan,
-    linspace,
-    ma,
-    mean as np_mean,
-    meshgrid,
-    nan,
-    nanmean,
-    ones,
-    percentile,
-    pi,
-    radians,
-    round_,
-    sin,
-    sinc,
-    where,
-    zeros,
-)
+from numba import njit, prange
+from numba import types as numba_types
+from numpy import (arange, array, ceil, concatenate, cos, deg2rad, empty,
+                   errstate, exp, float_, floor, histogram2d, int_, interp,
+                   isnan, linspace, ma)
+from numpy import mean as np_mean
+from numpy import (meshgrid, nan, nanmean, ones, percentile, pi, radians,
+                   round_, sin, sinc, where, zeros)
 from pint import UnitRegistry
 from scipy.interpolate import RectBivariateSpline, interp1d
 from scipy.ndimage import gaussian_filter
@@ -49,26 +24,15 @@
 from scipy.special import j1
 
 from .. import VAR_DESCR
+from ..data import get_demo_path
 from ..eddy_feature import Amplitude, Contours
-from ..generic import (
-    bbox_indice_regular,
-    coordinates_to_local,
-    distance,
-    interp2d_geo,
-    local_to_coordinates,
-    nearest_grd_indice,
-    uniform_resample,
-)
+from ..generic import (bbox_indice_regular, coordinates_to_local, distance,
+                       interp2d_geo, local_to_coordinates, nearest_grd_indice,
+                       uniform_resample)
 from ..observations.observation import EddiesObservations
-from ..poly import (
-    create_vertice,
-    fit_circle,
-    get_pixel_in_regular,
-    poly_area,
-    poly_contain_poly,
-    visvalingam,
-    winding_number_poly,
-)
+from ..poly import (create_vertice, fit_circle, get_pixel_in_regular,
+                    poly_area, poly_contain_poly, visvalingam,
+                    winding_number_poly)
 
 logger = logging.getLogger("pet")
 
@@ -1318,9 +1282,13 @@ def compute_pixel_path(self, x0, y0, x1, y1):
             self.x_size,
         )
 
-    def clean_land(self):
+    def clean_land(self, name):
         """Function to remove all land pixel"""
-        pass
+        mask_land = self.__class__(get_demo_path("mask_1_60.nc"), "lon", "lat")
+        x,y = meshgrid(self.x_c, self.y_c)
+        m = mask_land.interp('mask', x.reshape(-1), y.reshape(-1), 'nearest')
+        data = self.grid(name)
+        self.vars[name] = ma.array(data, mask=m.reshape(x.shape).T)
 
     def is_circular(self):
         """Check if the grid is circular"""
@@ -2392,6 +2360,15 @@ def __iter__(self):
         for _, d in self.datasets:
             yield d
 
+    @property
+    def time(self):
+        return array([t for t, _ in self.datasets])
+
+    @property
+    def period(self):
+        t = self.time
+        return t.min(), t.max()
+
     def __getitem__(self, item):
         for t, d in self.datasets:
             if t == item:

src/py_eddy_tracker/generic.py

@@ -3,27 +3,11 @@
 Tool method which use mostly numba
 """
 
-from numba import njit, prange, types as numba_types
-from numpy import (
-    absolute,
-    arcsin,
-    arctan2,
-    bool_,
-    cos,
-    empty,
-    floor,
-    histogram,
-    interp,
-    isnan,
-    linspace,
-    nan,
-    ones,
-    pi,
-    radians,
-    sin,
-    where,
-    zeros,
-)
+from numba import njit, prange
+from numba import types as numba_types
+from numpy import (absolute, arcsin, arctan2, bool_, cos, empty, floor,
+                   histogram, interp, isnan, linspace, nan, ones, pi, radians,
+                   sin, where, zeros)
 
 
 @njit(cache=True)
@@ -426,7 +410,7 @@ def split_line(x, y, i):
     """
     nb_jump = len(where(i[1:] - i[:-1] != 0)[0])
     nb_value = x.shape[0]
-    final_size = (nb_jump - 1) + nb_value
+    final_size = nb_jump + nb_value
     new_x = empty(final_size, dtype=x.dtype)
     new_y = empty(final_size, dtype=y.dtype)
     new_j = 0

src/py_eddy_tracker/observations/network.py

@@ -2,37 +2,26 @@
 """
 Class to create network of observations
 """
-from glob import glob
 import logging
 import time
+from glob import glob
 
 import netCDF4
-from numba import njit, types as nb_types
-from numba.typed import List
-from numpy import (
-    arange,
-    array,
-    bincount,
-    bool_,
-    concatenate,
-    empty,
-    nan,
-    ones,
-    percentile,
-    uint16,
-    uint32,
-    unique,
-    where,
-    zeros,
-)
 import zarr
+from numba import njit
+from numba import types as nb_types
+from numba.typed import List
+from numpy import (arange, array, bincount, bool_, concatenate, empty, nan,
+                   ones, percentile, uint16, uint32, unique, where, zeros)
 
 from ..dataset.grid import GridCollection
 from ..generic import build_index, wrap_longitude
 from ..poly import bbox_intersection, vertice_overlap
-from .groups import GroupEddiesObservations, get_missing_indices, particle_candidate
+from .groups import (GroupEddiesObservations, get_missing_indices,
+                     particle_candidate)
 from .observation import EddiesObservations
-from .tracking import TrackEddiesObservations, track_loess_filter, track_median_filter
+from .tracking import (TrackEddiesObservations, track_loess_filter,
+                       track_median_filter)
 
 logger = logging.getLogger("pet")
 
@@ -280,6 +269,15 @@ def longer_than(self, nb_day_min=-1, nb_day_max=-1):
         """
         Select network on time duration
 
+        :param int nb_day_min: Minimal number of days covered by one network, if negative -> not used
+        :param int nb_day_max: Maximal number of days covered by one network, if negative -> not used
+        """
+        return self.extract_with_mask(self.mask_longer_than(nb_day_min, nb_day_max))
+
+    def mask_longer_than(self, nb_day_min=-1, nb_day_max=-1):
+        """
+        Select network on time duration
+
         :param int nb_day_min: Minimal number of days covered by one network, if negative -> not used
         :param int nb_day_max: Maximal number of days covered by one network, if negative -> not used
         """
@@ -293,7 +291,7 @@ def longer_than(self, nb_day_min=-1, nb_day_max=-1):
                 continue
             if nb_day_min <= (ptp(t[i]) + 1) <= nb_day_max:
                 mask[i] = True
-        return self.extract_with_mask(mask)
+        return mask
 
     @classmethod
     def from_split_network(cls, group_dataset, indexs, **kwargs):
@@ -800,7 +798,7 @@ def display_timeline(
         if field is not None:
             field = self.parse_varname(field)
         for i, b0, b1 in self.iter_on("segment"):
-            x = self.time[i]
+            x = self.time_datetime64[i]
             if x.shape[0] == 0:
                 continue
             if field is None:
@@ -831,7 +829,7 @@ def event_timeline(self, ax, field=None, method=None, factor=1, colors_mode="rol
 
         # TODO : fill mappables dict
         y_seg = dict()
-        _time = self.time
+        _time = self.time_datetime64
 
         if field is not None and method != "all":
             for i, b0, _ in self.iter_on("segment"):
@@ -1011,7 +1009,7 @@ def scatter_timeline(
         if "c" not in kwargs:
             v = self.parse_varname(name)
             kwargs["c"] = v * factor
-        mappables["scatter"] = ax.scatter(self.time, y, **kwargs)
+        mappables["scatter"] = ax.scatter(self.time_datetime64, y, **kwargs)
         return mappables
 
     def event_map(self, ax, **kwargs):
@@ -1244,7 +1242,7 @@ def networks_mask(self, id_networks, segment=False):
 
     def networks(self, id_networks):
         return self.extract_with_mask(
-            generate_mask_from_ids(id_networks, self.track.size, *self.index_network)
+            generate_mask_from_ids(array(id_networks), self.track.size, *self.index_network)
         )
 
     @property
@@ -1423,10 +1421,10 @@ def plot(self, ax, ref=None, color_cycle=None, **kwargs):
         :param dict kwargs: keyword arguments for Axes.plot
         :return: a list of matplotlib mappables
         """
-        nb_colors = 0
-        if color_cycle is not None:
-            kwargs = kwargs.copy()
-            nb_colors = len(color_cycle)
+        kwargs = kwargs.copy()
+        if color_cycle is None:
+            color_cycle = self.COLORS
+        nb_colors = len(color_cycle)
         mappables = list()
         if "label" in kwargs:
             kwargs["label"] = self.format_label(kwargs["label"])