add a fonction to get distance to next obs

Author: AntSimi

doc/autodoc/tracking.rst

@@ -1,5 +1,5 @@
 Correspondances
-============
+===============
 
 .. automodule:: py_eddy_tracker.tracking
     :members:

doc/index.rst

@@ -42,6 +42,7 @@ Welcome to py-eddy-tracker's documentation!
    autodoc/observations
    autodoc/eddy_feature
    autodoc/featured_tracking
+   autodoc/tracking
    autodoc/poly
    autodoc/generic
 

examples/10_tracking_diagnostics/pet_propagation.py

@@ -5,29 +5,9 @@
 """
 from matplotlib import pyplot as plt
 from py_eddy_tracker.observations.tracking import TrackEddiesObservations
-from py_eddy_tracker.generic import distance
+from py_eddy_tracker.generic import cumsum_by_track
 import py_eddy_tracker_sample
-from numpy import arange, empty
-from numba import njit
-
-
-# %%
-# We will create a function compile with numba, to compute a field which contains curvilign distance
-@njit(cache=True)
-def cum_distance_by_track(distance, track):
-    tr_previous = 0
-    d_cum = 0
-    new_distance = empty(track.shape, dtype=distance.dtype)
-    for i in range(distance.shape[0]):
-        tr = track[i]
-        if i != 0 and tr != tr_previous:
-            d_cum = 0
-        new_distance[i] = d_cum
-        d_cum += distance[i]
-        tr_previous = tr
-    new_distance[i + 1] = d_cum
-    return new_distance
-
+from numpy import arange
 
 # %%
 # Load an experimental med atlas over a period of 26 years (1993-2019)
@@ -45,10 +25,8 @@ def cum_distance_by_track(distance, track):
 
 # %%
 # Compute curvilign distance
-d_a = distance(a.longitude[:-1], a.latitude[:-1], a.longitude[1:], a.latitude[1:])
-d_c = distance(c.longitude[:-1], c.latitude[:-1], c.longitude[1:], c.latitude[1:])
-d_a = cum_distance_by_track(d_a, a["track"]) / 1000.0
-d_c = cum_distance_by_track(d_c, c["track"]) / 1000.0
+d_a = cumsum_by_track(a.distance_to_next(), a.tracks) / 1000.0
+d_c = cumsum_by_track(c.distance_to_next(), c.tracks) / 1000.0
 
 # %%
 # Plot

notebooks/python_module/10_tracking_diagnostics/pet_propagation.ipynb

@@ -15,7 +15,7 @@
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "\nPropagation Histogram\n===================\n"
+        "\n# Propagation Histogram\n"
       ]
     },
     {
@@ -26,25 +26,7 @@
       },
       "outputs": [],
       "source": [
-        "from matplotlib import pyplot as plt\nfrom py_eddy_tracker.observations.tracking import TrackEddiesObservations\nfrom py_eddy_tracker.generic import distance\nimport py_eddy_tracker_sample\nfrom numpy import arange, empty\nfrom numba import njit"
-      ]
-    },
-    {
-      "cell_type": "markdown",
-      "metadata": {},
-      "source": [
-        "We will create a function compile with numba, to compute a field which contains curvilign distance\n\n"
-      ]
-    },
-    {
-      "cell_type": "code",
-      "execution_count": null,
-      "metadata": {
-        "collapsed": false
-      },
-      "outputs": [],
-      "source": [
-        "@njit(cache=True)\ndef cum_distance_by_track(distance, track):\n    tr_previous = 0\n    d_cum = 0\n    new_distance = empty(track.shape, dtype=distance.dtype)\n    for i in range(distance.shape[0]):\n        tr = track[i]\n        if i != 0 and tr != tr_previous:\n            d_cum = 0\n        new_distance[i] = d_cum\n        d_cum += distance[i]\n        tr_previous = tr\n    new_distance[i + 1] = d_cum\n    return new_distance"
+        "from matplotlib import pyplot as plt\nfrom py_eddy_tracker.observations.tracking import TrackEddiesObservations\nfrom py_eddy_tracker.generic import cumsum_by_track\nimport py_eddy_tracker_sample\nfrom numpy import arange"
       ]
     },
     {
@@ -98,7 +80,7 @@
       },
       "outputs": [],
       "source": [
-        "d_a = distance(a.longitude[:-1], a.latitude[:-1], a.longitude[1:], a.latitude[1:])\nd_c = distance(c.longitude[:-1], c.latitude[:-1], c.longitude[1:], c.latitude[1:])\nd_a = cum_distance_by_track(d_a, a[\"track\"]) / 1000.0\nd_c = cum_distance_by_track(d_c, c[\"track\"]) / 1000.0"
+        "d_a = cumsum_by_track(a.distance_to_next(), a.tracks) / 1000.0\nd_c = cumsum_by_track(c.distance_to_next(), c.tracks) / 1000.0"
       ]
     },
     {

src/py_eddy_tracker/__init__.py

@@ -25,6 +25,10 @@
 import logging
 import numpy
 import zarr
+from ._version import get_versions
+
+__version__ = get_versions()["version"]
+del get_versions
 
 
 def start_logger():
@@ -140,6 +144,18 @@ def parse_args(self, *args, **kwargs):
         nc_dims=("obs",),
         nc_attr=dict(),
     ),
+    distance_next=dict(
+        attr_name=None,
+        nc_name="distance_next",
+        nc_type="float32",
+        output_type="uint16",
+        scale_factor=50.,
+        nc_dims=("obs",),
+        nc_attr=dict(
+            long_name="Distance to next position",
+            units='m',
+        ),
+    ),
     virtual=dict(
         attr_name=None,
         nc_name="observation_flag",
@@ -622,8 +638,3 @@ def parse_args(self, *args, **kwargs):
     VAR_DESCR_inv[VAR_DESCR[key]["nc_name"]] = key
     for key_old in VAR_DESCR[key].get("old_nc_name", list()):
         VAR_DESCR_inv[key_old] = key
-
-from ._version import get_versions
-
-__version__ = get_versions()["version"]
-del get_versions

src/py_eddy_tracker/appli/eddies.py

@@ -109,7 +109,7 @@ def display_infos():
         print(f"-- {filename} -- ")
         if track:
             vars_ = vars.copy()
-            vars_.extend(('track', 'observation_number', 'observation_flag'))
+            vars_.extend(('track', 'observation_number', 'observation_flag', 'longitude'))
             e = TrackEddiesObservations.load_file(filename, include_vars=vars_)
         else:
             e = EddiesObservations.load_file(filename, include_vars=vars)

src/py_eddy_tracker/appli/gui.py

@@ -58,7 +58,7 @@ def setup(self, cmap="jet", nb_step=25, figsize=(8, 6)):
         # plot
         self.fig = pyplot.figure(figsize=figsize)
         t0, t1 = self.period
-        self.fig.suptitle(f'{t0} -> {t1}')
+        self.fig.suptitle(f"{t0} -> {t1}")
         self.ax = self.fig.add_axes((0.05, 0.05, 0.9, 0.9))
         self.ax.set_xlim(x_min, x_max), self.ax.set_ylim(y_min, y_max)
         self.ax.set_aspect("equal")

src/py_eddy_tracker/generic.py

@@ -145,6 +145,30 @@ def distance(lon0, lat0, lon1, lat1):
     return 6370997.0 * 2 * arctan2(a_val ** 0.5, (1 - a_val) ** 0.5)
 
 
+@njit(cache=True)
+def cumsum_by_track(field, track):
+    """
+    Cumsum by track
+
+    :param array field: data to sum
+    :pram array(int) track: id of track to separate data
+    :return: cumsum with a reset at each start of track
+    :rtype: array
+    """
+    tr_previous = 0
+    d_cum = 0
+    cumsum_array = empty(track.shape, dtype=field.dtype)
+    for i in range(field.shape[0]):
+        tr = track[i]
+        if tr != tr_previous:
+            d_cum = 0
+        cumsum_array[i] = d_cum
+        d_cum += field[i]
+        tr_previous = tr
+    cumsum_array[i + 1] = d_cum
+    return cumsum_array
+
+
 @njit(cache=True, fastmath=True)
 def interp2d_geo(x_g, y_g, z_g, m_g, x, y):
     """

src/py_eddy_tracker/observations/observation.py

@@ -60,7 +60,7 @@
 from tokenize import TokenError
 from tarfile import ExFileObject
 from matplotlib.path import Path as BasePath
-from .. import VAR_DESCR, VAR_DESCR_inv
+from .. import VAR_DESCR, VAR_DESCR_inv, __version__
 from ..generic import (
     distance_grid,
     distance,
@@ -1318,6 +1318,7 @@ def global_attr(self):
             Metadata_Conventions="Unidata Dataset Discovery v1.0",
             comment="Surface product; mesoscale eddies",
             framework_used="https://github.com/AntSimi/py-eddy-tracker",
+            framework_version=__version__,
             standard_name_vocabulary="NetCDF Climate and Forecast (CF) Metadata Convention Standard Name Table",
             rotation_type=self.sign_type,
         )

src/py_eddy_tracker/observations/tracking.py

@@ -38,14 +38,14 @@
     zeros,
     array,
     median,
-    histogram
+    histogram,
 )
 from datetime import datetime, timedelta
 from numba import njit
 from Polygon import Polygon
 from .observation import EddiesObservations
 from .. import VAR_DESCR_inv
-from ..generic import split_line, wrap_longitude, build_index
+from ..generic import split_line, wrap_longitude, build_index, distance, cumsum_by_track
 from ..poly import polygon_overlap, create_vertice_from_2darray
 
 
@@ -56,7 +56,7 @@ class TrackEddiesObservations(EddiesObservations):
     """Class to practice Tracking on observations
     """
 
-    __slots__ = ("__obs_by_track", "__first_index_of_track")
+    __slots__ = ("__obs_by_track", "__first_index_of_track", "__nb_track")
 
     ELEMENTS = [
         "lon",
@@ -85,6 +85,16 @@ def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
         self.__first_index_of_track = None
         self.__obs_by_track = None
+        self.__nb_track = None
+
+    @property
+    def nb_tracks(self):
+        """
+        Will count and send number of track
+        """
+        if self.__nb_track is None:
+            self.__nb_track = (self.nb_obs_by_track != 0).sum()
+        return self.__nb_track
 
     def __repr__(self):
         content = super().__repr__()
@@ -94,23 +104,56 @@ def __repr__(self):
         nb_obs = self.observations.shape[0]
         m = self["virtual"].astype("bool")
         nb_m = m.sum()
-        bins_t = (0, 20, 50, 100, 200, 1000, 10000)
+        bins_t = (1, 20, 50, 100, 200, 1000, 10000)
         nb_tracks_by_t = histogram(nb, bins=bins_t)[0]
         nb_obs_by_t = histogram(nb, bins=bins_t, weights=nb)[0]
-        pct_tracks_by_t = nb_tracks_by_t / nb_tracks_by_t.sum() * 100.
-        pct_obs_by_t = nb_obs_by_t / nb_obs_by_t.sum() * 100.
+        pct_tracks_by_t = nb_tracks_by_t / nb_tracks_by_t.sum() * 100.0
+        pct_obs_by_t = nb_obs_by_t / nb_obs_by_t.sum() * 100.0
+        d = self.distance_to_next() / 1000.0
+        cum_d = cumsum_by_track(d, self.tracks)
+        m_last = ones(d.shape, dtype="bool")
+        m_last[-1] = False
+        m_last[self.index_from_track[1:] - 1] = False
         content += f"""
-    | {nb.shape[0]} tracks ({nb.mean():.2f} obs/tracks, shorter {nb.min()} obs, longer {nb.max()} obs)
-    |   {nb_m} filled observations ({nb_m / nb.shape[0]:.2f} obs/tracks, {nb_m / nb_obs * 100:.2f} % of total)
-    |   Intepolated speed area : {self["speed_area"][m].sum() / period / 1e12:.2f} Mkm²/day
-    |   Intepolated effective area : {self["effective_area"][m].sum() / period / 1e12:.2f} Mkm²/day
+    | {self.nb_tracks} tracks ({
+        nb_obs / self.nb_tracks:.2f} obs/tracks, shorter {nb[nb!=0].min()} obs, longer {nb.max()} obs)
+    |   {nb_m} filled observations ({nb_m / self.nb_tracks:.2f} obs/tracks, {nb_m / nb_obs * 100:.2f} % of total)
+    |   Intepolated speed area      : {self["speed_area"][m].sum() / period / 1e12:.2f} Mkm²/day
+    |   Intepolated effective area  : {self["effective_area"][m].sum() / period / 1e12:.2f} Mkm²/day
+    |   Distance by day             : Mean {d[m_last].mean():.2f} , Median {median(d[m_last]):.2f} km/day
+    |   Distance by track           : Mean {cum_d[~m_last].mean():.2f} , Median {median(cum_d[~m_last]):.2f} km/track
     ----Distribution in lifetime:
     |   Lifetime (days  )      {self.box_display(bins_t)}
     |   Percent of tracks         : {self.box_display(pct_tracks_by_t)}
-    |   Percent of eddies         : {self.box_display(pct_obs_by_t)}
-"""
+    |   Percent of eddies         : {self.box_display(pct_obs_by_t)}"""
         return content
 
+    def add_distance(self):
+        """Add a field of distance (m) between to consecutive observation, 0 for the last observation of each track
+        """
+        if "distance_next" in self.observations.dtype.descr:
+            return self
+        new = self.add_fields(("distance_next",))
+        new["distance_next"][:1] = self.distance_to_next()
+        return new
+
+    def distance_to_next(self):
+        """
+        :return: array of distance in m, 0 when next obs if from another track
+        :rtype: array
+        """
+        d = distance(
+            self.longitude[:-1],
+            self.latitude[:-1],
+            self.longitude[1:],
+            self.latitude[1:],
+        )
+        d[self.index_from_track[1:] - 1] = 0
+        d_ = empty(d.shape[0] + 1, dtype=d.dtype)
+        d_[:-1] = d
+        d_[-1] = 0
+        return d_
+
     def filled_by_interpolation(self, mask):
         """Filled selected values by interpolation
         """
@@ -329,15 +372,14 @@ def __extract_with_mask(
     ):
         """
         Extract a subset of observations
-        Args:
-            mask: mask to select observations
-            full_path: extract full path if only one part is selected
-            remove_incomplete: delete path which are not fully selected
-            compress_id: resample track number to use a little range
-            reject_virtual: if track are only virtual in selection we remove track
-
-        Returns:
-            same object with selected observations
+
+        :param array(bool) mask: mask to select observations
+        :param full_path: extract full path if only one part is selected
+        :param remove_incomplete: delete path which are not fully selected
+        :param compress_id: resample track number to use a little range
+        :param reject_virtual: if track are only virtual in selection we remove track
+        :return: same object with selected observations
+        :rtype: self
         """
         if full_path and remove_incomplete:
             logger.warning(
@@ -373,6 +415,14 @@ def __extract_with_mask(
         return new
 
     def plot(self, ax, ref=None, **kwargs):
+        """
+        This function will draw path of each track
+
+        :param matplotlib.axes.Axes ax: ax where drawed
+        :param float,int ref: if defined all coordinates will be wrapped with ref like west boundary
+        :param dict kwargs: keyword arguments for Axes.plot
+        :return: matplotlib mappable
+        """
         if "label" in kwargs:
             kwargs["label"] += " (%s eddies)" % (self.nb_obs_by_track != 0).sum()
         x, y = split_line(self.longitude, self.latitude, self.tracks)
@@ -436,7 +486,15 @@ def split_network(self, intern=True, **kwargs):
         # )
 
     def set_tracks(self, x, y, ids, window):
-        # Will split one group in tracks
+        """
+        Will split one group in tracks
+
+        :param array x: coordinates of group
+        :param array y: coordinates of group
+        :param ndarray ids: several fields like time, group, ...
+        :param int windows: number of days where observations could missed
+        """
+
         time_index = build_index(ids["time"])
         nb = x.shape[0]
         used = zeros(nb, dtype="bool")