Add short example about loopers and eddies

examples/12_external_data/pet_drifter_loopers.py

@@ -1,11 +1,9 @@
 """
-[draft] Loopers vs eddies from altimetry
+Colocate looper with eddy from altimetry
 ========================================
 
-All loopers datas used in this example are a subset from the dataset describe in this article
+All loopers data used in this example are a subset from the dataset describe in this article
 [Lumpkin, R. : Global characteristics of coherent vortices from surface drifter trajectories](https://doi.org/10.1002/2015JC011435)
-
-
 """
 
 import re
@@ -71,65 +69,85 @@ def update_axes(ax, mappable=None):
 )
 
 # %%
-#
+# Global view
+# ===========
 ax = start_axes("All drifter available in med from Lumpkin dataset")
 loopers_med.plot(ax, lw=0.5, color="r", ref=-10)
 update_axes(ax)
 
 # %%
-# Only long period drifter
-long_drifter = loopers_med.extract_with_length((400, -1))
-ax = start_axes("Long period drifter")
-long_drifter.plot(ax, lw=0.5, color="r", ref=-10)
-update_axes(ax)
-print(np.unique(long_drifter.tracks))
-
-# %%
-drifter_1 = long_drifter.extract_ids((3588,))
-fig = plt.figure(figsize=(8, 4))
+# One segment of drifter
+# ======================
+#
+# Get a drifter segment (index used have no correspondance with original dataset).
+looper = loopers_med.extract_ids((3588,))
+fig = plt.figure(figsize=(16, 6))
 ax = fig.add_subplot(111, aspect="equal")
+looper.plot(ax, lw=0.5, label="Original position of drifter")
+looper_filtered = looper.copy()
+looper_filtered.position_filter(1, 13)
+s = looper_filtered.scatter(
+    ax,
+    "time",
+    cmap=plt.get_cmap("Spectral_r", 20),
+    label="Filtered position of drifter",
+)
+plt.colorbar(s).set_label("time (days from 1/1/1950)")
+ax.legend()
 ax.grid()
-drifter_1.plot(ax, lw=0.5)
-# %%
-drifter_1.close_tracks(
-    anticyclonic_eddies, method="close_center", delta=0.5, nb_obs_min=25
-).tracks
 
 # %%
-# Animation which show drifter with colocated eddies
-ed = anticyclonic_eddies.extract_ids((4738, 4836, 4910,))
-x, y, t = drifter_1.lon, drifter_1.lat, drifter_1.time
-
-
-def update(frame):
-    # We display last 5 days of loopers trajectory
-    m = (t < frame) * (t > (frame - 5))
-    a.func_animation(frame)
-    line.set_data(x[m], y[m])
-
-
-a = Anim(ed, intern=True, figsize=(8, 8), cmap="magma_r", nb_step=10, dpi=60)
-line = a.ax.plot([], [], "r", lw=4, zorder=100)[0]
-kwargs = dict(frames=np.arange(*a.period, 1), interval=100)
-a.fig.suptitle("")
-_ = VideoAnimation(a.fig, update, **kwargs)
+# Try to find a detected eddies with adt at same place. We used filtered track to simulate an eddy center
+match = looper_filtered.close_tracks(
+    anticyclonic_eddies, method="close_center", delta=0.1, nb_obs_min=50
+)
+fig = plt.figure(figsize=(16, 6))
+ax = fig.add_subplot(111, aspect="equal")
+looper.plot(ax, lw=0.5, label="Original position of drifter")
+looper_filtered.plot(ax, lw=1.5, label="Filtered position of drifter")
+match.plot(ax, lw=1.5, label="Matched eddy")
+ax.legend()
+ax.grid()
 
 # %%
-fig = plt.figure(figsize=(12, 6))
+# Display radius of this 2 datasets.
+fig = plt.figure(figsize=(20, 8))
 ax = fig.add_subplot(111)
-ax.plot(drifter_1.time, drifter_1.radius_s / 1e3, label="loopers")
-drifter_1.median_filter(7, "time", "radius_s", inplace=True)
-drifter_1.loess_filter(15, "time", "radius_s", inplace=True)
+ax.plot(looper.time, looper.radius_s / 1e3, label="loopers")
+looper_radius = looper.copy()
+looper_radius.median_filter(1, "time", "radius_s", inplace=True)
+looper_radius.loess_filter(13, "time", "radius_s", inplace=True)
+ax.plot(
+    looper_radius.time,
+    looper_radius.radius_s / 1e3,
+    label="loopers (filtered half window 13 days)",
+)
+ax.plot(match.time, match.radius_s / 1e3, label="altimetry")
+match_radius = match.copy()
+match_radius.median_filter(1, "time", "radius_s", inplace=True)
+match_radius.loess_filter(13, "time", "radius_s", inplace=True)
 ax.plot(
-    drifter_1.time,
-    drifter_1.radius_s / 1e3,
-    label="loopers (filtered half window 15 days)",
+    match_radius.time,
+    match_radius.radius_s / 1e3,
+    label="altimetry (filtered half window 13 days)",
 )
-ax.plot(ed.time, ed.radius_s / 1e3, label="altimetry")
-ed.median_filter(7, "time", "radius_s", inplace=True)
-ed.loess_filter(15, "time", "radius_s", inplace=True)
-ax.plot(ed.time, ed.radius_s / 1e3, label="altimetry (filtered half window 15 days)")
 ax.set_ylabel("radius(km)"), ax.set_ylim(0, 100)
 ax.legend()
+ax.set_title("Radius from loopers and altimeter")
 ax.grid()
-_ = ax.set_title("Radius from loopers and altimeter")
+
+
+# %%
+# Animation which show drifter with colocated eddy
+def update(frame):
+    # We display last 5 days of loopers trajectory
+    m = (looper.time < frame) * (looper.time > (frame - 5))
+    anim.func_animation(frame)
+    line.set_data(looper.lon[m], looper.lat[m])
+
+
+anim = Anim(match, intern=True, figsize=(8, 8), cmap="magma_r", nb_step=10, dpi=75)
+# mappable to show drifter in red
+line = anim.ax.plot([], [], "r", lw=4, zorder=100)[0]
+anim.fig.suptitle("")
+_ = VideoAnimation(anim.fig, update, frames=np.arange(*anim.period, 1), interval=125)

notebooks/python_module/12_external_data/pet_drifter_loopers.ipynb

@@ -0,0 +1,191 @@
+{
+  "cells": [
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "%matplotlib inline"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "\nColocate looper with eddy from altimetry\n========================================\n\nAll loopers data used in this example are a subset from the dataset describe in this article\n[Lumpkin, R. : Global characteristics of coherent vortices from surface drifter trajectories](https://doi.org/10.1002/2015JC011435)\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "import re\n\nimport numpy as np\nimport py_eddy_tracker_sample\nfrom matplotlib import pyplot as plt\nfrom matplotlib.animation import FuncAnimation\n\nfrom py_eddy_tracker import data\nfrom py_eddy_tracker.appli.gui import Anim\nfrom py_eddy_tracker.observations.tracking import TrackEddiesObservations"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "class VideoAnimation(FuncAnimation):\n    def _repr_html_(self, *args, **kwargs):\n        \"\"\"To get video in html and have a player\"\"\"\n        content = self.to_html5_video()\n        return re.sub(\n            r'width=\"[0-9]*\"\\sheight=\"[0-9]*\"', 'width=\"100%\" height=\"100%\"', content\n        )\n\n    def save(self, *args, **kwargs):\n        if args[0].endswith(\"gif\"):\n            # In this case gif is use to create thumbnail which are not use but consume same time than video\n            # So we create an empty file, to save time\n            with open(args[0], \"w\") as _:\n                pass\n            return\n        return super().save(*args, **kwargs)\n\n\ndef start_axes(title):\n    fig = plt.figure(figsize=(13, 5))\n    ax = fig.add_axes([0.03, 0.03, 0.90, 0.94], aspect=\"equal\")\n    ax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)\n    ax.set_title(title, weight=\"bold\")\n    return ax\n\n\ndef update_axes(ax, mappable=None):\n    ax.grid()\n    if mappable:\n        plt.colorbar(mappable, cax=ax.figure.add_axes([0.94, 0.05, 0.01, 0.9]))"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Load eddies dataset\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "cyclonic_eddies = TrackEddiesObservations.load_file(\n    py_eddy_tracker_sample.get_demo_path(\"eddies_med_adt_allsat_dt2018/Cyclonic.zarr\")\n)\nanticyclonic_eddies = TrackEddiesObservations.load_file(\n    py_eddy_tracker_sample.get_demo_path(\n        \"eddies_med_adt_allsat_dt2018/Anticyclonic.zarr\"\n    )\n)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Load loopers dataset\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "loopers_med = TrackEddiesObservations.load_file(\n    data.get_demo_path(\"loopers_lumpkin_med.nc\")\n)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Global view\n===========\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "ax = start_axes(\"All drifter available in med from Lumpkin dataset\")\nloopers_med.plot(ax, lw=0.5, color=\"r\", ref=-10)\nupdate_axes(ax)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "One segment of drifter\n======================\n\nGet a drifter segment (index used have no correspondance with original dataset).\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "looper = loopers_med.extract_ids((3588,))\nfig = plt.figure(figsize=(16, 6))\nax = fig.add_subplot(111, aspect=\"equal\")\nlooper.plot(ax, lw=0.5, label=\"Original position of drifter\")\nlooper_filtered = looper.copy()\nlooper_filtered.position_filter(1, 13)\ns = looper_filtered.scatter(\n    ax,\n    \"time\",\n    cmap=plt.get_cmap(\"Spectral_r\", 20),\n    label=\"Filtered position of drifter\",\n)\nplt.colorbar(s).set_label(\"time (days from 1/1/1950)\")\nax.legend()\nax.grid()"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Try to find a detected eddies with adt at same place. We used filtered track to simulate an eddy center\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "match = looper_filtered.close_tracks(\n    anticyclonic_eddies, method=\"close_center\", delta=0.1, nb_obs_min=50\n)\nfig = plt.figure(figsize=(16, 6))\nax = fig.add_subplot(111, aspect=\"equal\")\nlooper.plot(ax, lw=0.5, label=\"Original position of drifter\")\nlooper_filtered.plot(ax, lw=1.5, label=\"Filtered position of drifter\")\nmatch.plot(ax, lw=1.5, label=\"Matched eddy\")\nax.legend()\nax.grid()"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Display radius of this 2 datasets.\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "fig = plt.figure(figsize=(20, 8))\nax = fig.add_subplot(111)\nax.plot(looper.time, looper.radius_s / 1e3, label=\"loopers\")\nlooper_radius = looper.copy()\nlooper_radius.median_filter(1, \"time\", \"radius_s\", inplace=True)\nlooper_radius.loess_filter(13, \"time\", \"radius_s\", inplace=True)\nax.plot(\n    looper_radius.time,\n    looper_radius.radius_s / 1e3,\n    label=\"loopers (filtered half window 13 days)\",\n)\nax.plot(match.time, match.radius_s / 1e3, label=\"altimetry\")\nmatch_radius = match.copy()\nmatch_radius.median_filter(1, \"time\", \"radius_s\", inplace=True)\nmatch_radius.loess_filter(13, \"time\", \"radius_s\", inplace=True)\nax.plot(\n    match_radius.time,\n    match_radius.radius_s / 1e3,\n    label=\"altimetry (filtered half window 13 days)\",\n)\nax.set_ylabel(\"radius(km)\"), ax.set_ylim(0, 100)\nax.legend()\nax.set_title(\"Radius from loopers and altimeter\")\nax.grid()"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Animation which show drifter with colocated eddy\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "def update(frame):\n    # We display last 5 days of loopers trajectory\n    m = (looper.time < frame) * (looper.time > (frame - 5))\n    anim.func_animation(frame)\n    line.set_data(looper.lon[m], looper.lat[m])\n\n\nanim = Anim(match, intern=True, figsize=(8, 8), cmap=\"magma_r\", nb_step=10, dpi=75)\n# mappable to show drifter in red\nline = anim.ax.plot([], [], \"r\", lw=4, zorder=100)[0]\nanim.fig.suptitle(\"\")\n_ = VideoAnimation(anim.fig, update, frames=np.arange(*anim.period, 1), interval=125)"
+      ]
+    }
+  ],
+  "metadata": {
+    "kernelspec": {
+      "display_name": "Python 3",
+      "language": "python",
+      "name": "python3"
+    },
+    "language_info": {
+      "codemirror_mode": {
+        "name": "ipython",
+        "version": 3
+      },
+      "file_extension": ".py",
+      "mimetype": "text/x-python",
+      "name": "python",
+      "nbconvert_exporter": "python",
+      "pygments_lexer": "ipython3",
+      "version": "3.7.7"
+    }
+  },
+  "nbformat": 4,
+  "nbformat_minor": 0
+}

src/py_eddy_tracker/observations/network.py

@@ -1327,7 +1327,7 @@ def extract_light_with_mask(self, mask):
         )
         new.sign_type = self.sign_type
         if nb_obs == 0:
-            logger.warning("Empty dataset will be created")
+            logger.info("Empty dataset will be created")
         else:
             logger.info(
                 f"{nb_obs} observations will be extracted ({nb_obs / self.shape[0]:.3%})"
@@ -1353,7 +1353,7 @@ def extract_with_mask(self, mask):
         new = self.__class__.new_like(self, nb_obs)
         new.sign_type = self.sign_type
         if nb_obs == 0:
-            logger.warning("Empty dataset will be created")
+            logger.info("Empty dataset will be created")
         else:
             logger.debug(
                 f"{nb_obs} observations will be extracted ({nb_obs / self.shape[0]:.3%})"

src/py_eddy_tracker/observations/tracking.py

@@ -502,7 +502,7 @@ def extract_with_mask(
         new = self.__class__.new_like(self, nb_obs)
         new.sign_type = self.sign_type
         if nb_obs == 0:
-            logger.warning("Empty dataset will be created")
+            logger.info("Empty dataset will be created")
         else:
             for field in self.obs.dtype.descr:
                 logger.debug("Copy of field %s ...", field)
@@ -568,6 +568,9 @@ def close_tracks(self, other, nb_obs_min=10, **kwargs):
             It could be a costly operation for huge dataset
         """
         p0, p1 = self.period
+        p0_other, p1_other = other.period
+        if p1_other < p0 or p1 < p0_other:
+            return other.__class__.new_like(other, 0)
         indexs = list()
         for i_self, i_other, t0, t1 in self.align_on(other, bins=arange(p0, p1 + 2)):
             i, j, s = self.match(other, i_self=i_self, i_other=i_other, **kwargs)