Merge branch 'master' into validation_particle

Author: ludwigVonKoopa

examples/07_cube_manipulation/pet_fsle_med.py

@@ -14,7 +14,7 @@
 
 from matplotlib import pyplot as plt
 from numba import njit
-from numpy import arange, arctan2, empty, isnan, log2, ma, meshgrid, ones, pi, zeros
+from numpy import arange, arctan2, empty, isnan, log, ma, meshgrid, ones, pi, zeros
 
 from py_eddy_tracker import start_logger
 from py_eddy_tracker.data import get_demo_path
@@ -71,7 +71,7 @@ def check_p(x, y, flse, theta, m_set, m, dt, dist_init=0.02, dist_max=0.6):
             s2 = ((dxn + dye) ** 2 + (dxe - dyn) ** 2) * (
                 (dxn - dye) ** 2 + (dxe + dyn) ** 2
             )
-            flse[i] = 1 / (2 * dt) * log2(1 / (2 * dist_init ** 2) * (s1 + s2 ** 0.5))
+            flse[i] = 1 / (2 * dt) * log(1 / (2 * dist_init ** 2) * (s1 + s2 ** 0.5))
             theta[i] = arctan2(at1, at2 + s2) * 180 / pi
             # To know where value are set
             m_set[i] = False
@@ -180,7 +180,7 @@ def build_triplet(x, y, step=0.02):
 ax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)
 ax.set_aspect("equal")
 ax.set_title("Finite size lyapunov exponent", weight="bold")
-kw = dict(cmap="viridis_r", vmin=-15, vmax=0)
+kw = dict(cmap="viridis_r", vmin=-20, vmax=0)
 m = fsle_custom.display(ax, 1 / fsle_custom.grid("fsle"), **kw)
 ax.grid()
 _ = plt.colorbar(m, cax=fig.add_axes([0.94, 0.05, 0.01, 0.9]))

notebooks/python_module/07_cube_manipulation/pet_fsle_med.ipynb

@@ -15,7 +15,7 @@
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "\n# FSLE experiment in med\n\nExample to build Finite Size Lyapunov Exponents, parameter values must be adapted for your case.\n\nExample use a method similar to `AVISO flse`_\n\n    https://www.aviso.altimetry.fr/en/data/products/value-added-products/\n    fsle-finite-size-lyapunov-exponents/fsle-description.html\n"
+        "\nFSLE experiment in med\n======================\n\nExample to build Finite Size Lyapunov Exponents, parameter values must be adapted for your case.\n\nExample use a method similar to `AVISO flse`_\n\n    https://www.aviso.altimetry.fr/en/data/products/value-added-products/\n    fsle-finite-size-lyapunov-exponents/fsle-description.html\n"
       ]
     },
     {
@@ -26,14 +26,14 @@
       },
       "outputs": [],
       "source": [
-        "from matplotlib import pyplot as plt\nfrom numba import njit\nfrom numpy import arange, arctan2, empty, isnan, log2, ma, meshgrid, ones, pi, zeros\n\nfrom py_eddy_tracker import start_logger\nfrom py_eddy_tracker.data import get_demo_path\nfrom py_eddy_tracker.dataset.grid import GridCollection, RegularGridDataset\n\nstart_logger().setLevel(\"ERROR\")"
+        "from matplotlib import pyplot as plt\nfrom numba import njit\nfrom numpy import arange, arctan2, empty, isnan, log, ma, meshgrid, ones, pi, zeros\n\nfrom py_eddy_tracker import start_logger\nfrom py_eddy_tracker.data import get_demo_path\nfrom py_eddy_tracker.dataset.grid import GridCollection, RegularGridDataset\n\nstart_logger().setLevel(\"ERROR\")"
       ]
     },
     {
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "## ADT in med\n:py:meth:`~py_eddy_tracker.dataset.grid.GridCollection.from_netcdf_cube` method is\nmade for data stores in time cube, you could use also\n:py:meth:`~py_eddy_tracker.dataset.grid.GridCollection.from_netcdf_list` method to\nload data-cube from multiple file.\n\n"
+        "ADT in med\n----------\n:py:meth:`~py_eddy_tracker.dataset.grid.GridCollection.from_netcdf_cube` method is\nmade for data stores in time cube, you could use also \n:py:meth:`~py_eddy_tracker.dataset.grid.GridCollection.from_netcdf_list` method to\nload data-cube from multiple file.\n\n"
       ]
     },
     {
@@ -51,7 +51,7 @@
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "## Methods to compute FSLE\n\n"
+        "Methods to compute FSLE\n-----------------------\n\n"
       ]
     },
     {
@@ -62,14 +62,14 @@
       },
       "outputs": [],
       "source": [
-        "@njit(cache=True, fastmath=True)\ndef check_p(x, y, flse, theta, m_set, m, dt, dist_init=0.02, dist_max=0.6):\n    \"\"\"\n    Check if distance between eastern or northern particle to center particle is bigger than `dist_max`\n    \"\"\"\n    nb_p = x.shape[0] // 3\n    delta = dist_max ** 2\n    for i in range(nb_p):\n        i0 = i * 3\n        i_n = i0 + 1\n        i_e = i0 + 2\n        # If particle already set, we skip\n        if m[i0] or m[i_n] or m[i_e]:\n            continue\n        # Distance with north\n        dxn, dyn = x[i0] - x[i_n], y[i0] - y[i_n]\n        dn = dxn ** 2 + dyn ** 2\n        # Distance with east\n        dxe, dye = x[i0] - x[i_e], y[i0] - y[i_e]\n        de = dxe ** 2 + dye ** 2\n\n        if dn >= delta or de >= delta:\n            s1 = dn + de\n            at1 = 2 * (dxe * dxn + dye * dyn)\n            at2 = de - dn\n            s2 = ((dxn + dye) ** 2 + (dxe - dyn) ** 2) * (\n                (dxn - dye) ** 2 + (dxe + dyn) ** 2\n            )\n            flse[i] = 1 / (2 * dt) * log2(1 / (2 * dist_init ** 2) * (s1 + s2 ** 0.5))\n            theta[i] = arctan2(at1, at2 + s2) * 180 / pi\n            # To know where value are set\n            m_set[i] = False\n            # To stop particle advection\n            m[i0], m[i_n], m[i_e] = True, True, True\n\n\n@njit(cache=True)\ndef build_triplet(x, y, step=0.02):\n    \"\"\"\n    Triplet building for each position we add east and north point with defined step\n    \"\"\"\n    nb_x = x.shape[0]\n    x_ = empty(nb_x * 3, dtype=x.dtype)\n    y_ = empty(nb_x * 3, dtype=y.dtype)\n    for i in range(nb_x):\n        i0 = i * 3\n        i_n, i_e = i0 + 1, i0 + 2\n        x__, y__ = x[i], y[i]\n        x_[i0], y_[i0] = x__, y__\n        x_[i_n], y_[i_n] = x__, y__ + step\n        x_[i_e], y_[i_e] = x__ + step, y__\n    return x_, y_"
+        "@njit(cache=True, fastmath=True)\ndef check_p(x, y, flse, theta, m_set, m, dt, dist_init=0.02, dist_max=0.6):\n    \"\"\"\n    Check if distance between eastern or northern particle to center particle is bigger than `dist_max`\n    \"\"\"\n    nb_p = x.shape[0] // 3\n    delta = dist_max ** 2\n    for i in range(nb_p):\n        i0 = i * 3\n        i_n = i0 + 1\n        i_e = i0 + 2\n        # If particle already set, we skip\n        if m[i0] or m[i_n] or m[i_e]:\n            continue\n        # Distance with north\n        dxn, dyn = x[i0] - x[i_n], y[i0] - y[i_n]\n        dn = dxn ** 2 + dyn ** 2\n        # Distance with east\n        dxe, dye = x[i0] - x[i_e], y[i0] - y[i_e]\n        de = dxe ** 2 + dye ** 2\n\n        if dn >= delta or de >= delta:\n            s1 = dn + de\n            at1 = 2 * (dxe * dxn + dye * dyn)\n            at2 = de - dn\n            s2 = ((dxn + dye) ** 2 + (dxe - dyn) ** 2) * (\n                (dxn - dye) ** 2 + (dxe + dyn) ** 2\n            )\n            flse[i] = 1 / (2 * dt) * log(1 / (2 * dist_init ** 2) * (s1 + s2 ** 0.5))\n            theta[i] = arctan2(at1, at2 + s2) * 180 / pi\n            # To know where value are set\n            m_set[i] = False\n            # To stop particle advection\n            m[i0], m[i_n], m[i_e] = True, True, True\n\n\n@njit(cache=True)\ndef build_triplet(x, y, step=0.02):\n    \"\"\"\n    Triplet building for each position we add east and north point with defined step\n    \"\"\"\n    nb_x = x.shape[0]\n    x_ = empty(nb_x * 3, dtype=x.dtype)\n    y_ = empty(nb_x * 3, dtype=y.dtype)\n    for i in range(nb_x):\n        i0 = i * 3\n        i_n, i_e = i0 + 1, i0 + 2\n        x__, y__ = x[i], y[i]\n        x_[i0], y_[i0] = x__, y__\n        x_[i_n], y_[i_n] = x__, y__ + step\n        x_[i_e], y_[i_e] = x__ + step, y__\n    return x_, y_"
       ]
     },
     {
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "## Settings\n\n"
+        "Settings\n--------\n\n"
       ]
     },
     {
@@ -87,7 +87,7 @@
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "## Particles\n\n"
+        "Particles\n---------\n\n"
       ]
     },
     {
@@ -105,7 +105,7 @@
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "## FSLE\n\n"
+        "FSLE\n----\n\n"
       ]
     },
     {
@@ -123,7 +123,7 @@
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "## Display FSLE\n\n"
+        "Display FSLE\n------------\n\n"
       ]
     },
     {
@@ -134,14 +134,14 @@
       },
       "outputs": [],
       "source": [
-        "fig = plt.figure(figsize=(13, 5), dpi=150)\nax = fig.add_axes([0.03, 0.03, 0.90, 0.94])\nax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)\nax.set_aspect(\"equal\")\nax.set_title(\"Finite size lyapunov exponent\", weight=\"bold\")\nkw = dict(cmap=\"viridis_r\", vmin=-15, vmax=0)\nm = fsle_custom.display(ax, 1 / fsle_custom.grid(\"fsle\"), **kw)\nax.grid()\n_ = plt.colorbar(m, cax=fig.add_axes([0.94, 0.05, 0.01, 0.9]))"
+        "fig = plt.figure(figsize=(13, 5), dpi=150)\nax = fig.add_axes([0.03, 0.03, 0.90, 0.94])\nax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)\nax.set_aspect(\"equal\")\nax.set_title(\"Finite size lyapunov exponent\", weight=\"bold\")\nkw = dict(cmap=\"viridis_r\", vmin=-20, vmax=0)\nm = fsle_custom.display(ax, 1 / fsle_custom.grid(\"fsle\"), **kw)\nax.grid()\n_ = plt.colorbar(m, cax=fig.add_axes([0.94, 0.05, 0.01, 0.9]))"
       ]
     },
     {
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "## Display Theta\n\n"
+        "Display Theta\n-------------\n\n"
       ]
     },
     {
@@ -172,7 +172,7 @@
       "name": "python",
       "nbconvert_exporter": "python",
       "pygments_lexer": "ipython3",
-      "version": "3.7.9"
+      "version": "3.7.7"
     }
   },
   "nbformat": 4,

share/tracking.yaml

@@ -8,6 +8,6 @@ PATHS:
 TRACK_DURATION_MIN: 4
 VIRTUAL_LENGTH_MAX: 0
 
-#CLASS:
-#    MODULE: py_eddy_tracker.featured_tracking.old_tracker_reference
-#    CLASS: CheltonTracker
+CLASS:
+   MODULE: py_eddy_tracker.featured_tracking.area_tracker
+   CLASS: AreaTracker

src/py_eddy_tracker/dataset/grid.py

@@ -617,7 +617,10 @@ def eddy_identification(
         :param float,None precision: Truncate values at the defined precision in m
         :param str force_height_unit: Unit used for height unit
         :param str force_speed_unit: Unit used for speed unit
-        :param dict kwargs: Argument given to amplitude
+        :param dict kwargs: Arguments given to amplitude (mle, nb_step_min, nb_step_to_be_mle).
+            Look at :py:meth:`py_eddy_tracker.eddy_feature.Amplitude`
+            The amplitude threshold is given by `step*nb_step_min`
+
 
         :return: Return a list of 2 elements: Anticyclones and Cyclones
         :rtype: py_eddy_tracker.observations.observation.EddiesObservations

src/py_eddy_tracker/eddy_feature.py

@@ -65,7 +65,7 @@ def __init__(
         :param float contour_height:
         :param array data:
         :param float interval:
-        :param int mle: maximum number of local maxima in contour
+        :param int mle: maximum number of local extrema in contour
         :param int nb_step_min: number of intervals to consider an eddy
         :param int nb_step_to_be_mle: number of intervals to be considered as an another maxima
         """

src/py_eddy_tracker/generic.py

@@ -89,8 +89,7 @@ def build_index(groups):
             first_index[group - i0 + 1 : next_group - i0 + 1] = i + 1
     last_index = zeros(amplitude, dtype=numba_types.int_)
     last_index[:-1] = first_index[1:]
-    # + 2 because we iterate only until -2 and we want upper bound ( 1 + 1)
-    last_index[-1] = i + 2
+    last_index[-1] = len(groups)
     return first_index, last_index, i0
 
 
@@ -310,7 +309,7 @@ def uniform_resample(x_val, y_val, num_fac=2, fixed_size=None):
     :param array_like x_val: input x contour coordinates
     :param array_like y_val: input y contour coordinates
     :param int num_fac: factor to increase lengths of output coordinates
-    :param int,None fixed_size: if define, it will used to set sampling
+    :param int,None fixed_size: if defined, will be used to set sampling
     """
     nb = x_val.shape[0]
     # Get distances

src/py_eddy_tracker/observations/observation.py

@@ -2057,7 +2057,9 @@ def contains(self, x, y, intern=False):
         xname, yname = self.intern(intern)
         m = ~(isnan(x) + isnan(y))
         i = -ones(x.shape, dtype="i4")
-        i[m] = poly_indexs(x[m], y[m], self[xname], self[yname])
+
+        if x.size != 0 and m.any():
+            i[m] = poly_indexs(x[m], y[m], self[xname], self[yname])
         return i
 
     def inside(self, x, y, intern=False):

src/py_eddy_tracker/poly.py

@@ -854,8 +854,8 @@ def poly_indexs(x_p, y_p, x_c, y_c):
     """
     Index of contour for each postion inside a contour, -1 in case of no contour
 
-    :param array x_p: longitude to test (must be define, no nan)
-    :param array y_p: latitude to test (must be define, no nan)
+    :param array x_p: longitude to test (must be defined, no nan)
+    :param array y_p: latitude to test (must be defined, no nan)
     :param array x_c: longitude of contours
     :param array y_c: latitude of contours
     """