.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "python_module/07_cube_manipulation/pet_lavd_detection.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. or to run this example in your browser via Binder .. rst-class:: sphx-glr-example-title .. _sphx_glr_python_module_07_cube_manipulation_pet_lavd_detection.py: 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 .. GENERATED FROM PYTHON SOURCE LINES 21-34 .. code-block:: Python from datetime import datetime from matplotlib import pyplot as plt from numpy import arange, isnan, ma, meshgrid, zeros from py_eddy_tracker import start_logger from py_eddy_tracker.data import get_demo_path from py_eddy_tracker.dataset.grid import GridCollection, RegularGridDataset from py_eddy_tracker.gui import GUI_AXES start_logger().setLevel("ERROR") .. GENERATED FROM PYTHON SOURCE LINES 35-47 .. code-block:: Python 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) .. GENERATED FROM PYTHON SOURCE LINES 48-70 .. code-block:: Python 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=GUI_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") .. GENERATED FROM PYTHON SOURCE LINES 71-73 Data ---- .. GENERATED FROM PYTHON SOURCE LINES 73-89 .. code-block:: Python # Load data cube of 3 month c = GridCollection.from_netcdf_cube( get_demo_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 .. rst-class:: sphx-glr-script-out .. code-block:: none /home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/conda/latest/lib/python3.12/site-packages/numpy/ma/core.py:1020: RuntimeWarning: overflow encountered in multiply result = self.f(da, db, *args, **kwargs) .. GENERATED FROM PYTHON SOURCE LINES 90-92 Particles --------- .. GENERATED FROM PYTHON SOURCE LINES 92-112 .. code-block:: Python # 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, u_name="u", v_name="v") 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) .. rst-class:: sphx-glr-script-out .. code-block:: none 252256 particles advected .. GENERATED FROM PYTHON SOURCE LINES 113-115 LAVD forward (dynamic field) ---------------------------- .. GENERATED FROM PYTHON SOURCE LINES 115-129 .. code-block:: Python lavd = zeros(original_shape) lavd_ = lavd[m] p = c.advect(x0.copy(), y0.copy(), 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) .. image-sg:: /python_module/07_cube_manipulation/images/sphx_glr_pet_lavd_detection_001.png :alt: LAVD with a forward advection :srcset: /python_module/07_cube_manipulation/images/sphx_glr_pet_lavd_detection_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/conda/latest/lib/python3.12/site-packages/numpy/ma/core.py:4365: RuntimeWarning: invalid value encountered in multiply self._data.__imul__(other_data) .. GENERATED FROM PYTHON SOURCE LINES 130-132 LAVD backward (dynamic field) ----------------------------- .. GENERATED FROM PYTHON SOURCE LINES 132-146 .. code-block:: Python lavd = zeros(original_shape) lavd_ = lavd[m] p = c.advect(x0.copy(), y0.copy(), 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) .. image-sg:: /python_module/07_cube_manipulation/images/sphx_glr_pet_lavd_detection_002.png :alt: LAVD with a backward advection :srcset: /python_module/07_cube_manipulation/images/sphx_glr_pet_lavd_detection_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 147-149 LAVD forward (static field) --------------------------- .. GENERATED FROM PYTHON SOURCE LINES 149-163 .. code-block:: Python lavd = zeros(original_shape) lavd_ = lavd[m] p = t0_grid.advect(x0.copy(), y0.copy(), **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) .. image-sg:: /python_module/07_cube_manipulation/images/sphx_glr_pet_lavd_detection_003.png :alt: LAVD with a forward advection on a static velocity field :srcset: /python_module/07_cube_manipulation/images/sphx_glr_pet_lavd_detection_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 164-166 LAVD backward (static field) ---------------------------- .. GENERATED FROM PYTHON SOURCE LINES 166-180 .. code-block:: Python lavd = zeros(original_shape) lavd_ = lavd[m] p = t0_grid.advect(x0.copy(), y0.copy(), 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) .. image-sg:: /python_module/07_cube_manipulation/images/sphx_glr_pet_lavd_detection_004.png :alt: LAVD with a backward advection on a static velocity field :srcset: /python_module/07_cube_manipulation/images/sphx_glr_pet_lavd_detection_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 181-185 Contour detection ----------------- 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. .. GENERATED FROM PYTHON SOURCE LINES 185-219 .. code-block:: Python 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) .. image-sg:: /python_module/07_cube_manipulation/images/sphx_glr_pet_lavd_detection_005.png :alt: Detection of eddies with several method :srcset: /python_module/07_cube_manipulation/images/sphx_glr_pet_lavd_detection_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/conda/latest/lib/python3.12/site-packages/py_eddy_tracker/dataset/grid.py:1915: RuntimeWarning: invalid value encountered in sqrt self._speed_ev = sqrt(u * u + v * v) /home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/conda/latest/lib/python3.12/site-packages/numpy/lib/function_base.py:4824: UserWarning: Warning: 'partition' will ignore the 'mask' of the MaskedArray. arr.partition( /home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/conda/latest/lib/python3.12/site-packages/numpy/lib/function_base.py:4824: UserWarning: Warning: 'partition' will ignore the 'mask' of the MaskedArray. arr.partition( .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 15.071 seconds) .. _sphx_glr_download_python_module_07_cube_manipulation_pet_lavd_detection.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: binder-badge .. image:: images/binder_badge_logo.svg :target: https://mybinder.org/v2/gh/AntSimi/py-eddy-tracker/master?urlpath=lab/tree/notebooks/python_module/07_cube_manipulation/pet_lavd_detection.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: pet_lavd_detection.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: pet_lavd_detection.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: pet_lavd_detection.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_