update examples

Author: AntSimi

examples/06_grid_manipulation/pet_advect.py

@@ -4,9 +4,11 @@
 
 Dummy advection which use only static geostrophic current, which didn't resolve the complex circulation of the ocean.
 """
-import numpy as np
+import re
+
 from matplotlib import pyplot as plt
 from matplotlib.animation import FuncAnimation
+from numpy import arange, isnan, meshgrid, ones
 
 import py_eddy_tracker.gui
 from py_eddy_tracker import data
@@ -32,7 +34,7 @@
 fig = plt.figure(figsize=(10, 5))
 ax = fig.add_axes([0, 0, 1, 1], projection="full_axes")
 ax.set_xlim(19, 30), ax.set_ylim(31, 36.5), ax.grid()
-x, y = np.meshgrid(g.x_c, g.y_c)
+x, y = meshgrid(g.x_c, g.y_c)
 a.filled(ax, facecolors="r", alpha=0.1), c.filled(ax, facecolors="b", alpha=0.1)
 _ = ax.quiver(x.T, y.T, g.grid("u"), g.grid("v"), scale=20)
 
@@ -41,7 +43,10 @@
 class VideoAnimation(FuncAnimation):
     def _repr_html_(self, *args, **kwargs):
         """To get video in html and have a player"""
-        return self.to_html5_video()
+        content = self.to_html5_video()
+        return re.sub(
+            r'width="[0-9]*"\sheight="[0-9]*"', 'width="100%" height="100%"', content
+        )
 
     def save(self, *args, **kwargs):
         if args[0].endswith("gif"):
@@ -56,26 +61,27 @@ def save(self, *args, **kwargs):
 # %%
 # Anim
 # ----
-# Particules positions
-x, y = np.meshgrid(np.arange(13, 36, 0.125), np.arange(28, 40, 0.125))
+# Particles setup
+step_p = 1 / 8
+x, y = meshgrid(arange(13, 36, step_p), arange(28, 40, step_p))
 x, y = x.reshape(-1), y.reshape(-1)
 # Remove all original position that we can't advect at first place
-m = ~np.isnan(g.interp("u", x, y))
-x, y = x[m], y[m]
+m = ~isnan(g.interp("u", x, y))
+x0, y0 = x[m], y[m]
+x, y = x0.copy(), y0.copy()
 
+# %%
 # Movie properties
-kwargs = dict(frames=np.arange(51), interval=90)
+kwargs = dict(frames=arange(51), interval=100)
 kw_p = dict(nb_step=2, time_step=21600)
 frame_t = kw_p["nb_step"] * kw_p["time_step"] / 86400.0
 
 
-def anim_ax(generator, **kw):
+# %%
+# Method
+def anim_ax(**kw):
     t = 0
-    for _ in range(4):
-        generator.__next__()
-        t += frame_t
-
-    fig = plt.figure(figsize=(10, 5), dpi=64)
+    fig = plt.figure(figsize=(10, 5), dpi=55)
     ax = fig.add_axes([0, 0, 1, 1], projection="full_axes")
     ax.set_xlim(19, 30), ax.set_ylim(31, 36.5), ax.grid()
     a.filled(ax, facecolors="r", alpha=0.1), c.filled(ax, facecolors="b", alpha=0.1)
@@ -94,20 +100,74 @@ def update(i_frame, t_step):
 # %%
 # Filament forward
 # ^^^^^^^^^^^^^^^^
-p = g.filament(x, y, "u", "v", **kw_p, filament_size=4, rk4=True)
-fig, txt, l, t = anim_ax(p, lw=0.5)
+# Draw 3 last position in one path for each particles.,
+# it could be run backward with `backward=True` option in filament method
+p = g.filament(x, y, "u", "v", **kw_p, filament_size=3)
+fig, txt, l, t = anim_ax(lw=0.5)
 ani = VideoAnimation(fig, update, **kwargs, fargs=(frame_t,))
 
 # %%
-# Filament backward
+# Particle forward
 # ^^^^^^^^^^^^^^^^^
-p = g.filament(x, y, "u", "v", **kw_p, filament_size=4, backward=True, rk4=True)
-fig, txt, l, t = anim_ax(p, lw=0.5)
+# Forward advection of particles
+p = g.advect(x, y, "u", "v", **kw_p)
+fig, txt, l, t = anim_ax(ls="", marker=".", markersize=1)
+ani = VideoAnimation(fig, update, **kwargs, fargs=(frame_t,))
+
+# %%
+# We get last position and run backward until original position
+p = g.advect(x, y, "u", "v", **kw_p, backward=True)
+fig, txt, l, _ = anim_ax(ls="", marker=".", markersize=1)
 ani = VideoAnimation(fig, update, **kwargs, fargs=(-frame_t,))
 
 # %%
-# Particule forward
-# ^^^^^^^^^^^^^^^^^
-p = g.advect(x, y, "u", "v", **kw_p, rk4=True)
-fig, txt, l, t = anim_ax(p, ls="", marker=".", markersize=1)
-ani = VideoAnimation(fig, update, **kwargs, fargs=(frame_t,))
+# Particles stat
+# --------------
+
+# %%
+# Time_step settings
+# ^^^^^^^^^^^^^^^^^^
+# Dummy experiment to test advection precision, we run particles 50 days forward and backward ad different time step
+# and we measure distance between new positions and original positions.
+fig = plt.figure()
+ax = fig.add_subplot(111)
+ax.grid()
+kw = dict(
+    bins=arange(0, 50, 0.001),
+    cumulative=True,
+    weights=ones(x0.shape) / x0.shape[0] * 100.0,
+    histtype="step",
+)
+for time_step in (10800, 21600, 43200, 86400):
+    x, y = x0.copy(), y0.copy()
+    kw_advect = dict(nb_step=int(50 * 86400 / time_step), time_step=time_step)
+    g.advect(x, y, "u", "v", **kw_advect).__next__()
+    g.advect(x, y, "u", "v", **kw_advect, backward=True).__next__()
+    d = ((x - x0) ** 2 + (y - y0) ** 2) ** 0.5
+    ax.hist(d, **kw, label=f"{86400. / time_step:.0f} time step by day")
+ax.set_xlim(0, 0.25), ax.set_ylim(0, 100), ax.legend(loc="lower right")
+ax.set_title("Distance after 50 days forward and 50 days backward")
+ax.set_xlabel("Distance between original position and final position (in degrees)")
+_ = ax.set_ylabel("Percent of particles with distance lesser than")
+
+# %%
+# Time duration
+# ^^^^^^^^^^^^^
+# We keep same time_step but change time duration
+fig = plt.figure()
+ax = fig.add_subplot(111)
+ax.grid()
+time_step = 10800
+for duration in (5, 50, 100):
+    x, y = x0.copy(), y0.copy()
+    kw_advect = dict(nb_step=int(duration * 86400 / time_step), time_step=time_step)
+    g.advect(x, y, "u", "v", **kw_advect).__next__()
+    g.advect(x, y, "u", "v", **kw_advect, backward=True).__next__()
+    d = ((x - x0) ** 2 + (y - y0) ** 2) ** 0.5
+    ax.hist(d, **kw, label=f"Time duration {duration} days")
+ax.set_xlim(0, 0.25), ax.set_ylim(0, 100), ax.legend(loc="lower right")
+ax.set_title(
+    "Distance after N days forward and N days backward\nwith a time step of 1/8 days"
+)
+ax.set_xlabel("Distance between original position and final position (in degrees)")
+_ = ax.set_ylabel("Percent of particles with distance lesser than ")

examples/06_grid_manipulation/pet_lavd.py

@@ -7,8 +7,12 @@
 
 Method are described here:
 
-    - Abernathey, Ryan, and George Haller. " Transport by Lagrangian Vortices in the Eastern Pacific", Journal of Physical Oceanography 48, 3 (2018): 667-685, accessed Feb 16, 2021, https://doi.org/10.1175/JPO-D-17-0102.1
-    - `Transport by Coherent Lagrangian Vortices`_, R. Abernathey, Sinha A., Tarshish N., Liu T., Zhang C., Haller G., 2019, Talk a t the Sources and Sinks of Ocean Mesoscale Eddy Energy CLIVAR Workshop
+    - Abernathey, Ryan, and George Haller. "Transport by Lagrangian Vortices in the Eastern Pacific",
+    Journal of Physical Oceanography 48, 3 (2018): 667-685, accessed Feb 16, 2021,
+    https://doi.org/10.1175/JPO-D-17-0102.1
+    - `Transport by Coherent Lagrangian Vortices`_,
+    R. Abernathey, Sinha A., Tarshish N., Liu T., Zhang C., Haller G., 2019,
+    Talk a t the Sources and Sinks of Ocean Mesoscale Eddy Energy CLIVAR Workshop
 
 .. _Transport by Coherent Lagrangian Vortices:
     https://usclivar.org/sites/default/files/meetings/2019/presentations/Aberernathey_CLIVAR.pdf

examples/10_tracking_diagnostics/pet_center_count.py

@@ -62,7 +62,7 @@
 
 g_c.vars["count"] = ratio
 m = g_c.display(
-    ax_ratio, name="count", vmin=0.1, vmax=10, norm=LogNorm(), cmap="coolwarm_r"
+    ax_ratio, name="count", norm=LogNorm(vmin=0.1, vmax=10), cmap="coolwarm_r"
 )
 plt.colorbar(m, cax=fig.add_axes([0.94, 0.02, 0.01, 0.2]))
 

examples/10_tracking_diagnostics/pet_pixel_used.py

@@ -61,7 +61,7 @@
 
 g_c.vars["count"] = ratio
 m = g_c.display(
-    ax_ratio, name="count", vmin=0.1, vmax=10, norm=LogNorm(), cmap="coolwarm_r"
+    ax_ratio, name="count", norm=LogNorm(vmin=0.1, vmax=10), cmap="coolwarm_r"
 )
 plt.colorbar(m, cax=fig.add_axes([0.95, 0.02, 0.01, 0.2]))
 

examples/14_generic_tools/pet_fit_contour.py

@@ -21,6 +21,7 @@
 # Load example identification file
 a = EddiesObservations.load_file(data.get_path("Anticyclonic_20190223.nc"))
 
+
 # %%
 # Function to draw circle or ellips from parameter
 def build_circle(x0, y0, r):

examples/16_network/pet_atlas.py

@@ -6,18 +6,16 @@
 from numpy import ma
 
 import py_eddy_tracker.gui
+from py_eddy_tracker.data import get_path
 from py_eddy_tracker.observations.network import NetworkObservations
 
-n = NetworkObservations.load_file(
-    "/data/adelepoulle/work/Eddies/20201217_network_build/tracking/med/Anticyclonic_seg.nc"
-)
-
+n = NetworkObservations.load_file(get_path("Anticyclonic_seg.nc"))
 # %%
 # Parameters
 step = 1 / 10.0
 bins = ((-10, 37, step), (30, 46, step))
 kw_time = dict(cmap="terrain_r", factor=100.0 / n.nb_days, name="count")
-kw_ratio = dict(cmap=plt.get_cmap("magma_r", 10))
+kw_ratio = dict(cmap=plt.get_cmap("viridis", 10))
 
 
 # %%
@@ -149,16 +147,26 @@ def update_axes(ax, mappable=None):
     ax,
     **kw_ratio,
     vmin=0,
-    vmax=1,
+    vmax=5,
     name=g_all_merging.vars["count"] * 100.0 / g_all.vars["count"]
 )
 update_axes(ax, m).set_label("Pixel used in % all atlas")
 
+# %%
+# Merging in networks longer than 10 days, with dead end remove (shorter than 5 days)
+# -----------------------------------------------------------------------------------
+ax = start_axes("")
+merger = n10.remove_dead_end(nobs=10).merging_event()
+g_10_merging = merger.grid_count(bins)
+m = g_10_merging.display(ax, **kw_time, vmin=0, vmax=1)
+update_axes(ax, m).set_label("Pixel used in % of time")
+
 # %%
 # Merging in networks longer than 10 days
 # ---------------------------------------
 ax = start_axes("")
-g_10_merging = n10.merging_event().grid_count(bins)
+merger = n10.merging_event()
+g_10_merging = merger.grid_count(bins)
 m = g_10_merging.display(ax, **kw_time, vmin=0, vmax=1)
 update_axes(ax, m).set_label("Pixel used in % of time")
 # %%
@@ -167,10 +175,52 @@ def update_axes(ax, mappable=None):
     ax,
     **kw_ratio,
     vmin=0,
-    vmax=2,
+    vmax=5,
     name=ma.array(
         g_10_merging.vars["count"] * 100.0 / g_10.vars["count"],
         mask=g_10.vars["count"] < 365,
     )
 )
 update_axes(ax, m).set_label("Pixel used in % all atlas")
+
+# %%
+# All Spliting
+# ------------
+# Display the occurence of spliting events
+ax = start_axes("")
+g_all_spliting = n.spliting_event().grid_count(bins)
+m = g_all_spliting.display(ax, **kw_time, vmin=0, vmax=1)
+update_axes(ax, m).set_label("Pixel used in % of time")
+
+# %%
+# Ratio merging events / eddy presence
+ax = start_axes("")
+m = g_all_spliting.display(
+    ax,
+    **kw_ratio,
+    vmin=0,
+    vmax=5,
+    name=g_all_spliting.vars["count"] * 100.0 / g_all.vars["count"]
+)
+update_axes(ax, m).set_label("Pixel used in % all atlas")
+
+# %%
+# Spliting in networks longer than 10 days
+# ----------------------------------------
+ax = start_axes("")
+g_10_spliting = n10.spliting_event().grid_count(bins)
+m = g_10_spliting.display(ax, **kw_time, vmin=0, vmax=1)
+update_axes(ax, m).set_label("Pixel used in % of time")
+# %%
+ax = start_axes("")
+m = g_10_spliting.display(
+    ax,
+    **kw_ratio,
+    vmin=0,
+    vmax=5,
+    name=ma.array(
+        g_10_spliting.vars["count"] * 100.0 / g_10.vars["count"],
+        mask=g_10.vars["count"] < 365,
+    )
+)
+update_axes(ax, m).set_label("Pixel used in % all atlas")

examples/16_network/pet_ioannou_2017_case.py

@@ -12,13 +12,14 @@
 from datetime import datetime, timedelta
 
 import numpy as np
+from matplotlib import colors
 from matplotlib import pyplot as plt
 from matplotlib.animation import FuncAnimation
 from matplotlib.ticker import FuncFormatter
-from matplotlib import colors
 
 import py_eddy_tracker.gui
 from py_eddy_tracker.appli.gui import Anim
+from py_eddy_tracker.data import get_path
 from py_eddy_tracker.observations.network import NetworkObservations
 
 
@@ -71,9 +72,7 @@ def update_axes(ax, mappable=None):
 # We know the position and the time of a specific eddy
 #
 # `n.extract_with_mask` give us the corresponding network
-n = NetworkObservations.load_file(
-    "med/Anticyclonic_seg.nc"
-)
+n = NetworkObservations.load_file(get_path("Anticyclonic_seg.nc"))
 i = np.where(
     (n.lat > 33)
     * (n.lat < 34)

examples/16_network/pet_relative.py

@@ -55,7 +55,7 @@
 
 # %%
 # Timeline by mean Effective Radius
-# -----------------------
+# ---------------------------------
 # The factor argument is applied on the chosen field
 fig = plt.figure(figsize=(15, 5))
 ax = fig.add_axes([0.04, 0.04, 0.92, 0.92])
@@ -237,7 +237,7 @@
 
 # %%
 # Get birth event
-# ------------------
+# ---------------
 # Display the starting position of non-splitted eddies
 fig = plt.figure(figsize=(15, 8))
 ax = fig.add_axes([0.04, 0.06, 0.90, 0.88], projection="full_axes")
@@ -248,7 +248,7 @@
 
 # %%
 # Get death event
-# ------------------
+# ---------------
 # Display the last position of non-merged eddies
 fig = plt.figure(figsize=(15, 8))
 ax = fig.add_axes([0.04, 0.06, 0.90, 0.88], projection="full_axes")

examples/16_network/pet_replay_segmentation.py

@@ -11,13 +11,14 @@
 from datetime import datetime, timedelta
 
 import numpy as np
+from matplotlib import colors
 from matplotlib import pyplot as plt
 from matplotlib.animation import FuncAnimation
 from matplotlib.ticker import FuncFormatter
-from matplotlib import colors
 
 import py_eddy_tracker.gui
 from py_eddy_tracker.appli.gui import Anim
+from py_eddy_tracker.data import get_path
 from py_eddy_tracker.observations.network import NetworkObservations
 from py_eddy_tracker.observations.tracking import TrackEddiesObservations
 
@@ -120,9 +121,7 @@ def get_obs(dataset):
 
 # %%
 # Get original network, we will isolate only relative at order *2*
-n = NetworkObservations.load_file(
-    "/data/adelepoulle/work/Eddies/20201217_network_build/tracking/med/Anticyclonic_seg.nc"
-)
+n = NetworkObservations.load_file(get_path("Anticyclonic_seg.nc"))
 
 n = n.extract_with_mask(n.track == n.track[get_obs(n)])
 n_ = n.relative(get_obs(n), order=2)