Add another LAVD example

Author: AntSimi

examples/07_cube_manipulation/pet_lavd_detection.py

@@ -0,0 +1,218 @@
+"""
+LAVD detection and geometric detection
+======================================
+
+Naive method to reproduce LAVD(Lagrangian-Averaged Vorticity deviation).
+In the current example we didn't remove a mean vorticity.
+
+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
+
+.. _Transport by Coherent Lagrangian Vortices:
+    https://usclivar.org/sites/default/files/meetings/2019/presentations/Aberernathey_CLIVAR.pdf
+
+"""
+from datetime import datetime
+
+from matplotlib import pyplot as plt
+from numpy import arange, isnan, ma, meshgrid, zeros
+
+import py_eddy_tracker.gui
+from py_eddy_tracker import start_logger
+from py_eddy_tracker.data import get_path
+from py_eddy_tracker.dataset.grid import GridCollection, RegularGridDataset
+
+start_logger().setLevel("ERROR")
+
+
+# %%
+class LAVDGrid(RegularGridDataset):
+    def init_speed_coef(self, uname="u", vname="v"):
+        """Hack to be able to identify eddy with LAVD field"""
+        self._speed_ev = self.grid("lavd")
+
+    @classmethod
+    def from_(cls, x, y, z):
+        z.mask += isnan(z.data)
+        datas = dict(lavd=z, lon=x, lat=y)
+        return cls.with_array(coordinates=("lon", "lat"), datas=datas, centered=True)
+
+
+# %%
+def start_ax(title="", dpi=90):
+    fig = plt.figure(figsize=(12, 5), dpi=dpi)
+    ax = fig.add_axes([0.05, 0.08, 0.9, 0.9], projection="full_axes")
+    ax.set_xlim(-6, 36), ax.set_ylim(31, 45)
+    ax.set_title(title)
+    return fig, ax, ax.text(3, 32, "", fontsize=20)
+
+
+def update_axes(ax, mappable=None):
+    ax.grid()
+    if mappable:
+        cb = plt.colorbar(
+            mappable,
+            cax=ax.figure.add_axes([0.05, 0.1, 0.9, 0.01]),
+            orientation="horizontal",
+        )
+        cb.set_label("LAVD at initial position")
+        return cb
+
+
+kw_lavd = dict(vmin=0, vmax=2e-5, cmap="viridis")
+
+# %%
+# Data
+# ----
+
+# Load data cube of 3 month
+c = GridCollection.from_netcdf_cube(
+    get_path("dt_med_allsat_phy_l4_2005T2.nc"),
+    "longitude",
+    "latitude",
+    "time",
+    heigth="adt",
+)
+
+# Add vorticity at each time step
+for g in c:
+    u_y = g.compute_stencil(g.grid("u"), vertical=True)
+    v_x = g.compute_stencil(g.grid("v"))
+    g.vars["vort"] = v_x - u_y
+
+# %%
+# Particles
+# ---------
+
+# Time properties, for example with advection only 25 days
+nb_days, step_by_day = 25, 6
+nb_time = step_by_day * nb_days
+kw_p = dict(nb_step=1, time_step=86400 / step_by_day)
+t0 = 20236
+t0_grid = c[t0]
+# Geographic properties, we use a coarser resolution for time consuming reasons
+step = 1 / 32.0
+x_g, y_g = arange(-6, 36, step), arange(30, 46, step)
+x0, y0 = meshgrid(x_g, y_g)
+original_shape = x0.shape
+x0, y0 = x0.reshape(-1), y0.reshape(-1)
+# Get all particles in defined area
+m = ~isnan(t0_grid.interp("vort", x0, y0))
+x0, y0 = x0[m], y0[m]
+print(f"{x0.size} particles advected")
+# Gridded mask
+m = m.reshape(original_shape)
+
+# %%
+# LAVD forward (dynamic field)
+# ----------------------------
+lavd = zeros(original_shape)
+lavd_ = lavd[m]
+p = c.advect(x0.copy(), y0.copy(), "u", "v", t_init=t0, **kw_p)
+for _ in range(nb_time):
+    t, x, y = p.__next__()
+    lavd_ += abs(c.interp("vort", t / 86400.0, x, y))
+lavd[m] = lavd_ / nb_time
+# Put LAVD result in a standard py eddy tracker grid
+lavd_forward = LAVDGrid.from_(x_g, y_g, ma.array(lavd, mask=~m).T)
+# Display
+fig, ax, _ = start_ax("LAVD with a forward advection")
+mappable = lavd_forward.display(ax, "lavd", **kw_lavd)
+_ = update_axes(ax, mappable)
+
+# %%
+# LAVD backward (dynamic field)
+# -----------------------------
+lavd = zeros(original_shape)
+lavd_ = lavd[m]
+p = c.advect(x0.copy(), y0.copy(), "u", "v", t_init=t0, backward=True, **kw_p)
+for i in range(nb_time):
+    t, x, y = p.__next__()
+    lavd_ += abs(c.interp("vort", t / 86400.0, x, y))
+lavd[m] = lavd_ / nb_time
+# Put LAVD result in a standard py eddy tracker grid
+lavd_backward = LAVDGrid.from_(x_g, y_g, ma.array(lavd, mask=~m).T)
+# Display
+fig, ax, _ = start_ax("LAVD with a backward advection")
+mappable = lavd_backward.display(ax, "lavd", **kw_lavd)
+_ = update_axes(ax, mappable)
+
+# %%
+# LAVD forward (static field)
+# ---------------------------
+lavd = zeros(original_shape)
+lavd_ = lavd[m]
+p = t0_grid.advect(x0.copy(), y0.copy(), "u", "v", **kw_p)
+for _ in range(nb_time):
+    x, y = p.__next__()
+    lavd_ += abs(t0_grid.interp("vort", x, y))
+lavd[m] = lavd_ / nb_time
+# Put LAVD result in a standard py eddy tracker grid
+lavd_forward_static = LAVDGrid.from_(x_g, y_g, ma.array(lavd, mask=~m).T)
+# Display
+fig, ax, _ = start_ax("LAVD with a forward advection on a static velocity field")
+mappable = lavd_forward_static.display(ax, "lavd", **kw_lavd)
+_ = update_axes(ax, mappable)
+
+# %%
+# LAVD backward (static field)
+# ----------------------------
+lavd = zeros(original_shape)
+lavd_ = lavd[m]
+p = t0_grid.advect(x0.copy(), y0.copy(), "u", "v", backward=True, **kw_p)
+for i in range(nb_time):
+    x, y = p.__next__()
+    lavd_ += abs(t0_grid.interp("vort", x, y))
+lavd[m] = lavd_ / nb_time
+# Put LAVD result in a standard py eddy tracker grid
+lavd_backward_static = LAVDGrid.from_(x_g, y_g, ma.array(lavd, mask=~m).T)
+# Display
+fig, ax, _ = start_ax("LAVD with a backward advection on a static velocity field")
+mappable = lavd_backward_static.display(ax, "lavd", **kw_lavd)
+_ = update_axes(ax, mappable)
+
+# %%
+# Contourdetection
+# ----------------
+# To extract contour from LAVD grid, we will used method design for SSH, with some hacks and adapted options.
+# It will produce false amplitude and speed.
+kw_ident = dict(
+    force_speed_unit="m/s",
+    force_height_unit="m",
+    pixel_limit=(40, 200000),
+    date=datetime(2005, 5, 18),
+    uname=None,
+    vname=None,
+    grid_height="lavd",
+    shape_error=70,
+    step=1e-6,
+)
+fig, ax, _ = start_ax("Detection of eddies with several method")
+t0_grid.bessel_high_filter("adt", 700)
+a, c = t0_grid.eddy_identification(
+    "adt", "u", "v", kw_ident["date"], step=0.002, shape_error=70
+)
+kw_ed = dict(ax=ax, intern=True, ref=-10)
+a.filled(
+    facecolors="#FFEFCD", label="Anticyclonic SSH detection {nb_obs} eddies", **kw_ed
+)
+c.filled(facecolors="#DEDEDE", label="Cyclonic SSH detection {nb_obs} eddies", **kw_ed)
+kw_cont = dict(ax=ax, extern_only=True, ls="-", ref=-10)
+forward, _ = lavd_forward.eddy_identification(**kw_ident)
+forward.display(label="LAVD forward {nb_obs} eddies", color="g", **kw_cont)
+backward, _ = lavd_backward.eddy_identification(**kw_ident)
+backward.display(label="LAVD backward {nb_obs} eddies", color="r", **kw_cont)
+forward, _ = lavd_forward_static.eddy_identification(**kw_ident)
+forward.display(label="LAVD forward static {nb_obs} eddies", color="cyan", **kw_cont)
+backward, _ = lavd_backward_static.eddy_identification(**kw_ident)
+backward.display(
+    label="LAVD backward static {nb_obs} eddies", color="orange", **kw_cont
+)
+ax.legend()
+update_axes(ax)

examples/16_network/pet_relative.py

@@ -136,54 +136,46 @@
 ax.set_title(f"Clean network ({n_clean.infos()})")
 _ = n_clean.display_timeline(ax)
 
+# %%
+# For further figure we will use clean path
+n = n_clean
 
 # %%
-# change splittint-merging events
-# ------------------
+# Change splitting-merging events
+# -------------------------------
 # change event where seg A split to B, then A merge into B, to A split to B then B merge into A
-fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 12), dpi=120)
-
-ax1.set_title(f"Clean network ({n_clean.infos()})")
-n_clean.display_timeline(ax1)
+fig = plt.figure(figsize=(15, 12))
+ax = fig.add_axes([0.04, 0.54, 0.90, 0.40])
+ax.set_title(f"Clean network ({n.infos()})")
+n.display_timeline(ax)
 
-clean_modified = n_clean.copy()
-clean_modified.correct_close_events(100)
-ax2.set_title(f"resplitted network ({clean_modified.infos()})")
-_ = clean_modified.display_timeline(ax2)
+clean_modified = n.copy()
+# If it's happen in less than 40 days
+clean_modified.correct_close_events(40)
 
+ax = fig.add_axes([0.04, 0.04, 0.90, 0.40])
+ax.set_title(f"resplitted network ({clean_modified.infos()})")
+_ = clean_modified.display_timeline(ax)
 
 # %%
-# keep only observations where water could propagate from an observation
+# Keep only observations where water could propagate from an observation
 # ----------------------------------------------------------------------
-
-fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 12), dpi=120)
 i_observation = 600
-only_linked = n_clean.find_link(i_observation)
-
-for ax, dataset in zip([ax1, ax2], [n_clean, only_linked]):
-    dataset.display_timeline(
-        ax, field="segment", marker="+", lw=2, markersize=5, colors_mode="y"
-    )
-    ax.scatter(
-        n_clean.time[i_observation],
-        n_clean.segment[i_observation],
-        marker="s",
-        s=50,
-        color="black",
-        zorder=200,
-        label="observation start",
-        alpha=0.6,
-    )
+only_linked = n.find_link(i_observation)
+
+fig = plt.figure(figsize=(15, 12))
+ax1 = fig.add_axes([0.04, 0.54, 0.90, 0.40])
+ax2 = fig.add_axes([0.04, 0.04, 0.90, 0.40])
+
+kw = dict(marker="s", s=300, color="black", zorder=200, label="observation start")
+for ax, dataset in zip([ax1, ax2], [n, only_linked]):
+    dataset.display_timeline(ax, field="segment", lw=2, markersize=5, colors_mode="y")
+    ax.scatter(n.time[i_observation], n.segment[i_observation], **kw)
     ax.legend()
 
-ax1.set_title(f"full example ({n_clean.infos()})")
+ax1.set_title(f"full example ({n.infos()})")
 ax2.set_title(f"only linked observations ({only_linked.infos()})")
-ax2.set_xlim(ax1.get_xlim())
-ax2.set_ylim(ax1.get_ylim())
-
-# %%
-# For further figure we will use clean path
-n = n_clean
+_ = ax2.set_xlim(ax1.get_xlim()), ax2.set_ylim(ax1.get_ylim())
 
 # %%
 # Keep close relative
@@ -299,7 +291,6 @@
 ax.set_xlim(17.5, 27.5), ax.set_ylim(31, 36), ax.grid()
 m1
 
-
 # %%
 # Get spliting event
 # ------------------

notebooks/README.md

@@ -1,4 +1,3 @@
-# rm build/sphinx/ doc/python_module/ doc/gen_modules/ -rf
-# rm doc/_autosummary/ doc/gen_modules -rf
+# rm build/sphinx/ doc/python_module/ doc/gen_modules/ doc/_autosummary/ -rf
 python setup.py build_sphinx
 rsync -vrltp doc/python_module notebooks/.   --include '*/' --include '*.ipynb' --exclude '*' --prune-empty-dirs