update example of propagation and lifetime

Author: AntSimi

examples/06_grid_manipulation/README.rst

@@ -1,2 +1,2 @@
 Grid Manipulation
-=======================
+=================

examples/08_tracking_manipulation/README.rst

@@ -1,2 +1,2 @@
 Tracking Manipulation
-=======================
+=====================

examples/10_tracking_diagnostics/README.rst

@@ -1,2 +1,2 @@
 Tracking diagnostics
-=======================
+====================

examples/10_tracking_diagnostics/pet_lifetime.py

@@ -6,7 +6,7 @@
 from matplotlib import pyplot as plt
 from py_eddy_tracker.observations.tracking import TrackEddiesObservations
 import py_eddy_tracker_sample
-from numpy import arange
+from numpy import arange, ones
 
 # %%
 # Load an experimental med atlas over a period of 26 years (1993-2019)
@@ -16,39 +16,41 @@
 c = TrackEddiesObservations.load_file(
     py_eddy_tracker_sample.get_path("eddies_med_adt_allsat_dt2018/Cyclonic.zarr")
 )
+nb_year = (a.period[1] - a.period[0] + 1) / 365.25
 
 # %%
-# Plot
-fig = plt.figure()
-ax_lifetime = fig.add_axes([0.05, 0.55, 0.4, 0.4])
-ax_cum_lifetime = fig.add_axes([0.55, 0.55, 0.4, 0.4])
-ax_ratio_lifetime = fig.add_axes([0.05, 0.05, 0.4, 0.4])
-ax_ratio_cum_lifetime = fig.add_axes([0.55, 0.05, 0.4, 0.4])
-
-cum_a, bins, _ = ax_cum_lifetime.hist(
-    a["n"], histtype="step", bins=arange(0, 800, 1), label="Anticyclonic", color="r"
-)
-cum_c, bins, _ = ax_cum_lifetime.hist(
-    c["n"], histtype="step", bins=arange(0, 800, 1), label="Cyclonic", color="b"
-)
-
-x = (bins[1:] + bins[:-1]) / 2.0
-ax_ratio_cum_lifetime.plot(x, cum_c / cum_a)
-
-nb_a, nb_c = cum_a[:-1] - cum_a[1:], cum_c[:-1] - cum_c[1:]
-ax_lifetime.plot(x[1:], nb_a, label="Anticyclonic", color="r")
-ax_lifetime.plot(x[1:], nb_c, label="Cyclonic", color="b")
-
-ax_ratio_lifetime.plot(x[1:], nb_c / nb_a)
-
-for ax in (ax_lifetime, ax_cum_lifetime, ax_ratio_cum_lifetime, ax_ratio_lifetime):
-    ax.set_xlim(0, 365)
-    if ax in (ax_lifetime, ax_cum_lifetime):
-        ax.set_ylim(1, None)
-        ax.set_yscale("log")
-        ax.legend()
+# Setup axes
+figure = plt.figure(figsize=(12, 8))
+ax_ratio_cum = figure.add_axes([0.55, 0.06, 0.42, 0.34])
+ax_ratio = figure.add_axes([0.07, 0.06, 0.46, 0.34])
+ax_cum = figure.add_axes([0.55, 0.43, 0.42, 0.54])
+ax = figure.add_axes([0.07, 0.43, 0.46, 0.54])
+ax.set_ylabel("Eddies by year")
+ax_ratio.set_ylabel("Ratio Cyclonic/Anticyclonic")
+for ax_ in (ax, ax_cum, ax_ratio_cum, ax_ratio):
+    ax_.set_xlim(0, 400)
+    if ax_ in (ax, ax_cum):
+        ax_.set_ylim(1e-1, 1e4), ax_.set_yscale("log")
     else:
-        ax.set_ylim(0, 2)
-        ax.set_ylabel("Ratio Cyclonic/Anticyclonic")
-    ax.set_xlabel("Lifetime (days)")
-    ax.grid()
+        ax_.set_xlabel("Lifetime in days (by week bins)")
+        ax_.set_ylim(0, 2)
+        ax_.axhline(1, color="g", lw=2)
+    ax_.grid()
+ax_cum.xaxis.set_ticklabels([]), ax_cum.yaxis.set_ticklabels([])
+ax.xaxis.set_ticklabels([]), ax_ratio_cum.yaxis.set_ticklabels([])
+
+# plot data
+bin_hist = arange(7, 2000, 7)
+x = (bin_hist[1:] + bin_hist[:-1]) / 2.0
+a_nb, c_nb = a.nb_obs_by_track, c.nb_obs_by_track
+a_nb, c_nb = a_nb[a_nb != 0], c_nb[c_nb != 0]
+w_a, w_c = ones(a_nb.shape) / nb_year, ones(c_nb.shape) / nb_year
+kwargs_a = dict(histtype="step", bins=bin_hist, x=a_nb, color="r", weights=w_a)
+kwargs_c = dict(histtype="step", bins=bin_hist, x=c_nb, color="b", weights=w_c)
+cum_a, _, _ = ax_cum.hist(cumulative=-1, **kwargs_a)
+cum_c, _, _ = ax_cum.hist(cumulative=-1, **kwargs_c)
+nb_a, _, _ = ax.hist(label="Anticyclonic", **kwargs_a)
+nb_c, _, _ = ax.hist(label="Cyclonic", **kwargs_c)
+ax_ratio_cum.plot(x, cum_c / cum_a)
+ax_ratio.plot(x, nb_c / nb_a)
+ax.legend()

examples/10_tracking_diagnostics/pet_propagation.py

@@ -7,7 +7,7 @@
 from py_eddy_tracker.observations.tracking import TrackEddiesObservations
 from py_eddy_tracker.generic import cumsum_by_track
 import py_eddy_tracker_sample
-from numpy import arange
+from numpy import arange, ones
 
 # %%
 # Load an experimental med atlas over a period of 26 years (1993-2019)
@@ -17,6 +17,7 @@
 c = TrackEddiesObservations.load_file(
     py_eddy_tracker_sample.get_path("eddies_med_adt_allsat_dt2018/Cyclonic.zarr")
 )
+nb_year = (a.period[1] - a.period[0] + 1) / 365.25
 
 # %%
 # Filtering position to remove noisy position
@@ -25,47 +26,42 @@
 
 # %%
 # Compute curvilign distance
-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
+i0, nb = a.index_from_track, a.nb_obs_by_track
+d_a = cumsum_by_track(a.distance_to_next(), a.tracks)[(i0 - 1 + nb)[nb != 0]] / 1000.0
+i0, nb = c.index_from_track, c.nb_obs_by_track
+d_c = cumsum_by_track(c.distance_to_next(), c.tracks)[(i0 - 1 + nb)[nb != 0]] / 1000.0
 
 # %%
-# Plot
-fig = plt.figure()
-ax_propagation = fig.add_axes([0.05, 0.55, 0.4, 0.4])
-ax_cum_propagation = fig.add_axes([0.55, 0.55, 0.4, 0.4])
-ax_ratio_propagation = fig.add_axes([0.05, 0.05, 0.4, 0.4])
-ax_ratio_cum_propagation = fig.add_axes([0.55, 0.05, 0.4, 0.4])
-
-bins = arange(0, 1500, 10)
-cum_a, bins, _ = ax_cum_propagation.hist(
-    d_a, histtype="step", bins=bins, label="Anticyclonic", color="r"
-)
-cum_c, bins, _ = ax_cum_propagation.hist(
-    d_c, histtype="step", bins=bins, label="Cyclonic", color="b"
-)
-
-x = (bins[1:] + bins[:-1]) / 2.0
-ax_ratio_cum_propagation.plot(x, cum_c / cum_a)
-
-nb_a, nb_c = cum_a[:-1] - cum_a[1:], cum_c[:-1] - cum_c[1:]
-ax_propagation.plot(x[1:], nb_a, label="Anticyclonic", color="r")
-ax_propagation.plot(x[1:], nb_c, label="Cyclonic", color="b")
-
-ax_ratio_propagation.plot(x[1:], nb_c / nb_a)
-
-for ax in (
-    ax_propagation,
-    ax_cum_propagation,
-    ax_ratio_cum_propagation,
-    ax_ratio_propagation,
-):
-    ax.set_xlim(0, 1000)
-    if ax in (ax_propagation, ax_cum_propagation):
-        ax.set_ylim(1, None)
-        ax.set_yscale("log")
-        ax.legend()
+# Setup axes
+figure = plt.figure(figsize=(12, 8))
+ax_ratio_cum = figure.add_axes([0.55, 0.06, 0.42, 0.34])
+ax_ratio = figure.add_axes([0.07, 0.06, 0.46, 0.34])
+ax_cum = figure.add_axes([0.55, 0.43, 0.42, 0.54])
+ax = figure.add_axes([0.07, 0.43, 0.46, 0.54])
+ax.set_ylabel("Eddies by year")
+ax_ratio.set_ylabel("Ratio Cyclonic/Anticyclonic")
+for ax_ in (ax, ax_cum, ax_ratio_cum, ax_ratio):
+    ax_.set_xlim(0, 1000)
+    if ax_ in (ax, ax_cum):
+        ax_.set_ylim(1e-1, 1e4), ax_.set_yscale("log")
     else:
-        ax.set_ylim(0, 2)
-        ax.set_ylabel("Ratio Cyclonic/Anticyclonic")
-    ax.set_xlabel("Propagation (km)")
-    ax.grid()
+        ax_.set_xlabel("Propagation in km (with bins of 20 km)")
+        ax_.set_ylim(0, 2)
+        ax_.axhline(1, color="g", lw=2)
+    ax_.grid()
+ax_cum.xaxis.set_ticklabels([]), ax_cum.yaxis.set_ticklabels([])
+ax.xaxis.set_ticklabels([]), ax_ratio_cum.yaxis.set_ticklabels([])
+
+# plot data
+bin_hist = arange(0, 2000, 20)
+x = (bin_hist[1:] + bin_hist[:-1]) / 2.0
+w_a, w_c = ones(d_a.shape) / nb_year, ones(d_c.shape) / nb_year
+kwargs_a = dict(histtype="step", bins=bin_hist, x=d_a, color="r", weights=w_a)
+kwargs_c = dict(histtype="step", bins=bin_hist, x=d_c, color="b", weights=w_c)
+cum_a, _, _ = ax_cum.hist(cumulative=-1, **kwargs_a)
+cum_c, _, _ = ax_cum.hist(cumulative=-1, **kwargs_c)
+nb_a, _, _ = ax.hist(label="Anticyclonic", **kwargs_a)
+nb_c, _, _ = ax.hist(label="Cyclonic", **kwargs_c)
+ax_ratio_cum.plot(x, cum_c / cum_a)
+ax_ratio.plot(x, nb_c / nb_a)
+ax.legend()

notebooks/python_module/10_tracking_diagnostics/pet_lifetime.ipynb

@@ -15,7 +15,7 @@
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "\nLifetime Histogram\n===================\n"
+        "\n# Lifetime Histogram\n"
       ]
     },
     {
@@ -26,7 +26,7 @@
       },
       "outputs": [],
       "source": [
-        "from matplotlib import pyplot as plt\nfrom py_eddy_tracker.observations.tracking import TrackEddiesObservations\nimport py_eddy_tracker_sample\nfrom numpy import arange"
+        "from matplotlib import pyplot as plt\nfrom py_eddy_tracker.observations.tracking import TrackEddiesObservations\nimport py_eddy_tracker_sample\nfrom numpy import arange, ones"
       ]
     },
     {
@@ -44,14 +44,14 @@
       },
       "outputs": [],
       "source": [
-        "a = TrackEddiesObservations.load_file(\n    py_eddy_tracker_sample.get_path(\"eddies_med_adt_allsat_dt2018/Anticyclonic.zarr\")\n)\nc = TrackEddiesObservations.load_file(\n    py_eddy_tracker_sample.get_path(\"eddies_med_adt_allsat_dt2018/Cyclonic.zarr\")\n)"
+        "a = TrackEddiesObservations.load_file(\n    py_eddy_tracker_sample.get_path(\"eddies_med_adt_allsat_dt2018/Anticyclonic.zarr\")\n)\nc = TrackEddiesObservations.load_file(\n    py_eddy_tracker_sample.get_path(\"eddies_med_adt_allsat_dt2018/Cyclonic.zarr\")\n)\nnb_year = (a.period[1] - a.period[0] + 1) / 365.25"
       ]
     },
     {
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "Plot\n\n"
+        "Setup axes\n\n"
       ]
     },
     {
@@ -62,7 +62,7 @@
       },
       "outputs": [],
       "source": [
-        "fig = plt.figure()\nax_lifetime = fig.add_axes([0.05, 0.55, 0.4, 0.4])\nax_cum_lifetime = fig.add_axes([0.55, 0.55, 0.4, 0.4])\nax_ratio_lifetime = fig.add_axes([0.05, 0.05, 0.4, 0.4])\nax_ratio_cum_lifetime = fig.add_axes([0.55, 0.05, 0.4, 0.4])\n\ncum_a, bins, _ = ax_cum_lifetime.hist(\n    a[\"n\"], histtype=\"step\", bins=arange(0, 800, 1), label=\"Anticyclonic\", color=\"r\"\n)\ncum_c, bins, _ = ax_cum_lifetime.hist(\n    c[\"n\"], histtype=\"step\", bins=arange(0, 800, 1), label=\"Cyclonic\", color=\"b\"\n)\n\nx = (bins[1:] + bins[:-1]) / 2.0\nax_ratio_cum_lifetime.plot(x, cum_c / cum_a)\n\nnb_a, nb_c = cum_a[:-1] - cum_a[1:], cum_c[:-1] - cum_c[1:]\nax_lifetime.plot(x[1:], nb_a, label=\"Anticyclonic\", color=\"r\")\nax_lifetime.plot(x[1:], nb_c, label=\"Cyclonic\", color=\"b\")\n\nax_ratio_lifetime.plot(x[1:], nb_c / nb_a)\n\nfor ax in (ax_lifetime, ax_cum_lifetime, ax_ratio_cum_lifetime, ax_ratio_lifetime):\n    ax.set_xlim(0, 365)\n    if ax in (ax_lifetime, ax_cum_lifetime):\n        ax.set_ylim(1, None)\n        ax.set_yscale(\"log\")\n        ax.legend()\n    else:\n        ax.set_ylim(0, 2)\n        ax.set_ylabel(\"Ratio Cyclonic/Anticyclonic\")\n    ax.set_xlabel(\"Lifetime (days)\")\n    ax.grid()"
+        "figure = plt.figure(figsize=(12, 8))\nax_ratio_cum = figure.add_axes([0.55, 0.06, 0.42, 0.34])\nax_ratio = figure.add_axes([0.07, 0.06, 0.46, 0.34])\nax_cum = figure.add_axes([0.55, 0.43, 0.42, 0.54])\nax = figure.add_axes([0.07, 0.43, 0.46, 0.54])\nax.set_ylabel(\"Eddies by year\")\nax_ratio.set_ylabel(\"Ratio Cyclonic/Anticyclonic\")\nfor ax_ in (ax, ax_cum, ax_ratio_cum, ax_ratio):\n    ax_.set_xlim(0, 400)\n    if ax_ in (ax, ax_cum):\n        ax_.set_ylim(1e-1, 1e4), ax_.set_yscale(\"log\")\n    else:\n        ax_.set_xlabel(\"Lifetime in days (by week bins)\")\n        ax_.set_ylim(0, 2)\n        ax_.axhline(1, color=\"g\", lw=2)\n    ax_.grid()\nax_cum.xaxis.set_ticklabels([]), ax_cum.yaxis.set_ticklabels([])\nax.xaxis.set_ticklabels([]), ax_ratio_cum.yaxis.set_ticklabels([])\n\n# plot data\nbin_hist = arange(7, 2000, 7)\nx = (bin_hist[1:] + bin_hist[:-1]) / 2.0\na_nb, c_nb = a.nb_obs_by_track, c.nb_obs_by_track\na_nb, c_nb = a_nb[a_nb != 0], c_nb[c_nb != 0]\nw_a, w_c = ones(a_nb.shape) / nb_year, ones(c_nb.shape) / nb_year\nkwargs_a = dict(histtype=\"step\", bins=bin_hist, x=a_nb, color=\"r\", weights=w_a)\nkwargs_c = dict(histtype=\"step\", bins=bin_hist, x=c_nb, color=\"b\", weights=w_c)\ncum_a, _, _ = ax_cum.hist(cumulative=-1, **kwargs_a)\ncum_c, _, _ = ax_cum.hist(cumulative=-1, **kwargs_c)\nnb_a, _, _ = ax.hist(label=\"Anticyclonic\", **kwargs_a)\nnb_c, _, _ = ax.hist(label=\"Cyclonic\", **kwargs_c)\nax_ratio_cum.plot(x, cum_c / cum_a)\nax_ratio.plot(x, nb_c / nb_a)\nax.legend()"
       ]
     }
   ],

notebooks/python_module/10_tracking_diagnostics/pet_propagation.ipynb

@@ -26,7 +26,7 @@
       },
       "outputs": [],
       "source": [
-        "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"
+        "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, ones"
       ]
     },
     {
@@ -44,7 +44,7 @@
       },
       "outputs": [],
       "source": [
-        "a = TrackEddiesObservations.load_file(\n    py_eddy_tracker_sample.get_path(\"eddies_med_adt_allsat_dt2018/Anticyclonic.zarr\")\n)\nc = TrackEddiesObservations.load_file(\n    py_eddy_tracker_sample.get_path(\"eddies_med_adt_allsat_dt2018/Cyclonic.zarr\")\n)"
+        "a = TrackEddiesObservations.load_file(\n    py_eddy_tracker_sample.get_path(\"eddies_med_adt_allsat_dt2018/Anticyclonic.zarr\")\n)\nc = TrackEddiesObservations.load_file(\n    py_eddy_tracker_sample.get_path(\"eddies_med_adt_allsat_dt2018/Cyclonic.zarr\")\n)\nnb_year = (a.period[1] - a.period[0] + 1) / 365.25"
       ]
     },
     {
@@ -80,14 +80,14 @@
       },
       "outputs": [],
       "source": [
-        "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"
+        "i0, nb = a.index_from_track, a.nb_obs_by_track\nd_a = cumsum_by_track(a.distance_to_next(), a.tracks)[(i0 - 1 + nb)[nb != 0]] / 1000.0\ni0, nb = c.index_from_track, c.nb_obs_by_track\nd_c = cumsum_by_track(c.distance_to_next(), c.tracks)[(i0 - 1 + nb)[nb != 0]] / 1000.0"
       ]
     },
     {
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "Plot\n\n"
+        "Setup axes\n\n"
       ]
     },
     {
@@ -98,7 +98,7 @@
       },
       "outputs": [],
       "source": [
-        "fig = plt.figure()\nax_propagation = fig.add_axes([0.05, 0.55, 0.4, 0.4])\nax_cum_propagation = fig.add_axes([0.55, 0.55, 0.4, 0.4])\nax_ratio_propagation = fig.add_axes([0.05, 0.05, 0.4, 0.4])\nax_ratio_cum_propagation = fig.add_axes([0.55, 0.05, 0.4, 0.4])\n\nbins = arange(0, 1500, 10)\ncum_a, bins, _ = ax_cum_propagation.hist(\n    d_a, histtype=\"step\", bins=bins, label=\"Anticyclonic\", color=\"r\"\n)\ncum_c, bins, _ = ax_cum_propagation.hist(\n    d_c, histtype=\"step\", bins=bins, label=\"Cyclonic\", color=\"b\"\n)\n\nx = (bins[1:] + bins[:-1]) / 2.0\nax_ratio_cum_propagation.plot(x, cum_c / cum_a)\n\nnb_a, nb_c = cum_a[:-1] - cum_a[1:], cum_c[:-1] - cum_c[1:]\nax_propagation.plot(x[1:], nb_a, label=\"Anticyclonic\", color=\"r\")\nax_propagation.plot(x[1:], nb_c, label=\"Cyclonic\", color=\"b\")\n\nax_ratio_propagation.plot(x[1:], nb_c / nb_a)\n\nfor ax in (\n    ax_propagation,\n    ax_cum_propagation,\n    ax_ratio_cum_propagation,\n    ax_ratio_propagation,\n):\n    ax.set_xlim(0, 1000)\n    if ax in (ax_propagation, ax_cum_propagation):\n        ax.set_ylim(1, None)\n        ax.set_yscale(\"log\")\n        ax.legend()\n    else:\n        ax.set_ylim(0, 2)\n        ax.set_ylabel(\"Ratio Cyclonic/Anticyclonic\")\n    ax.set_xlabel(\"Propagation (km)\")\n    ax.grid()"
+        "figure = plt.figure(figsize=(12, 8))\nax_ratio_cum = figure.add_axes([0.55, 0.06, 0.42, 0.34])\nax_ratio = figure.add_axes([0.07, 0.06, 0.46, 0.34])\nax_cum = figure.add_axes([0.55, 0.43, 0.42, 0.54])\nax = figure.add_axes([0.07, 0.43, 0.46, 0.54])\nax.set_ylabel(\"Eddies by year\")\nax_ratio.set_ylabel(\"Ratio Cyclonic/Anticyclonic\")\nfor ax_ in (ax, ax_cum, ax_ratio_cum, ax_ratio):\n    ax_.set_xlim(0, 1000)\n    if ax_ in (ax, ax_cum):\n        ax_.set_ylim(1e-1, 1e4), ax_.set_yscale(\"log\")\n    else:\n        ax_.set_xlabel(\"Propagation in km (with bins of 20 km)\")\n        ax_.set_ylim(0, 2)\n        ax_.axhline(1, color=\"g\", lw=2)\n    ax_.grid()\nax_cum.xaxis.set_ticklabels([]), ax_cum.yaxis.set_ticklabels([])\nax.xaxis.set_ticklabels([]), ax_ratio_cum.yaxis.set_ticklabels([])\n\n# plot data\nbin_hist = arange(0, 2000, 20)\nx = (bin_hist[1:] + bin_hist[:-1]) / 2.0\nw_a, w_c = ones(d_a.shape) / nb_year, ones(d_c.shape) / nb_year\nkwargs_a = dict(histtype=\"step\", bins=bin_hist, x=d_a, color=\"r\", weights=w_a)\nkwargs_c = dict(histtype=\"step\", bins=bin_hist, x=d_c, color=\"b\", weights=w_c)\ncum_a, _, _ = ax_cum.hist(cumulative=-1, **kwargs_a)\ncum_c, _, _ = ax_cum.hist(cumulative=-1, **kwargs_c)\nnb_a, _, _ = ax.hist(label=\"Anticyclonic\", **kwargs_a)\nnb_c, _, _ = ax.hist(label=\"Cyclonic\", **kwargs_c)\nax_ratio_cum.plot(x, cum_c / cum_a)\nax_ratio.plot(x, nb_c / nb_a)\nax.legend()"
       ]
     }
   ],