{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# LAVD experiment\n\nNaive method to reproduce LAVD(Lagrangian-Averaged Vorticity deviation) method with a static velocity field.\nIn the current example we didn't remove a mean vorticity.\n\nMethod are described here:\n\n - Abernathey, Ryan, and George Haller. \"Transport by Lagrangian Vortices in the Eastern Pacific\",\n Journal of Physical Oceanography 48, 3 (2018): 667-685, accessed Feb 16, 2021,\n https://doi.org/10.1175/JPO-D-17-0102.1\n - `Transport by Coherent Lagrangian Vortices`_,\n R. Abernathey, Sinha A., Tarshish N., Liu T., Zhang C., Haller G., 2019,\n Talk a t the Sources and Sinks of Ocean Mesoscale Eddy Energy CLIVAR Workshop\n\n https://usclivar.org/sites/default/files/meetings/2019/presentations/Aberernathey_CLIVAR.pdf\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import re\n\nfrom matplotlib import pyplot as plt\nfrom matplotlib.animation import FuncAnimation\nfrom numpy import arange, meshgrid, zeros\n\nfrom py_eddy_tracker.data import get_demo_path\nfrom py_eddy_tracker.dataset.grid import RegularGridDataset\nfrom py_eddy_tracker.gui import GUI_AXES\nfrom py_eddy_tracker.observations.observation import EddiesObservations" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def start_ax(title=\"\", dpi=90):\n fig = plt.figure(figsize=(16, 9), dpi=dpi)\n ax = fig.add_axes([0, 0, 1, 1], projection=GUI_AXES)\n ax.set_xlim(0, 32), ax.set_ylim(28, 46)\n ax.set_title(title)\n return fig, ax, ax.text(3, 32, \"\", fontsize=20)\n\n\ndef update_axes(ax, mappable=None):\n ax.grid()\n if mappable:\n cb = plt.colorbar(\n mappable,\n cax=ax.figure.add_axes([0.05, 0.1, 0.9, 0.01]),\n orientation=\"horizontal\",\n )\n cb.set_label(\"Vorticity integration along trajectory at initial position\")\n return cb\n\n\nkw_vorticity = dict(vmin=0, vmax=2e-5, cmap=\"viridis\")" ] }, { "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 used to create thumbnail which is not used 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)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Data\nTo compute vorticity ($\\omega$) we compute u/v field with a stencil and apply the following equation with stencil\nmethod :\n\n\\begin{align}\\omega = \\frac{\\partial v}{\\partial x} - \\frac{\\partial u}{\\partial y}\\end{align}\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "g = RegularGridDataset(\n get_demo_path(\"dt_med_allsat_phy_l4_20160515_20190101.nc\"), \"longitude\", \"latitude\"\n)\ng.add_uv(\"adt\")\nu_y = g.compute_stencil(g.grid(\"u\"), vertical=True)\nv_x = g.compute_stencil(g.grid(\"v\"))\ng.vars[\"vort\"] = v_x - u_y" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Display vorticity field\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig, ax, _ = start_ax()\nmappable = g.display(ax, abs(g.grid(\"vort\")), **kw_vorticity)\ncb = update_axes(ax, mappable)\ncb.set_label(\"Vorticity\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Particles\nParticles specification\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "step = 1 / 32\nx_g, y_g = arange(0, 36, step), arange(28, 46, step)\nx, y = meshgrid(x_g, y_g)\noriginal_shape = x.shape\nx, y = x.reshape(-1), y.reshape(-1)\nprint(f\"{len(x)} particles advected\")\n# A frame every 8h\nstep_by_day = 3\n# Compute step of advection every 4h\nnb_step = 2\nkw_p = dict(nb_step=nb_step, time_step=86400 / step_by_day / nb_step, u_name=\"u\", v_name=\"v\")\n# Start a generator which at each iteration return new position at next time step\nparticule = g.advect(x, y, **kw_p, rk4=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## LAVD\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "lavd = zeros(original_shape)\n# Advection time\nnb_days = 8\n# Nb frame\nnb_time = step_by_day * nb_days\ni = 0.0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Anim\nMovie of LAVD integration at each integration time step.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def update(i_frame):\n global lavd, i\n i += 1\n x, y = particule.__next__()\n # Interp vorticity on new_position\n lavd += abs(g.interp(\"vort\", x, y).reshape(original_shape) * 1 / nb_time)\n txt.set_text(f\"T0 + {i / step_by_day:.2f} days of advection\")\n pcolormesh.set_array(lavd / i * nb_time)\n return pcolormesh, txt\n\n\nkw_video = dict(frames=arange(nb_time), interval=1000.0 / step_by_day / 2, blit=True)\nfig, ax, txt = start_ax(dpi=60)\nx_g_, y_g_ = (\n arange(0 - step / 2, 36 + step / 2, step),\n arange(28 - step / 2, 46 + step / 2, step),\n)\n# pcolorfast will be faster than pcolormesh, we could use pcolorfast due to x and y are regular\npcolormesh = ax.pcolorfast(x_g_, y_g_, lavd, **kw_vorticity)\nupdate_axes(ax, pcolormesh)\n_ = VideoAnimation(ax.figure, update, **kw_video)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Final LAVD\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Format LAVD data\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "lavd = RegularGridDataset.with_array(\n coordinates=(\"lon\", \"lat\"), datas=dict(lavd=lavd.T, lon=x_g, lat=y_g), centered=True\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Display final LAVD with py eddy tracker detection.\nPeriod used for LAVD integration (8 days) is too short for a real use, but choose for example efficiency.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig, ax, _ = start_ax()\nmappable = lavd.display(ax, \"lavd\", **kw_vorticity)\nEddiesObservations.load_file(get_demo_path(\"Anticyclonic_20160515.nc\")).display(\n ax, color=\"k\"\n)\nEddiesObservations.load_file(get_demo_path(\"Cyclonic_20160515.nc\")).display(\n ax, color=\"k\"\n)\n_ = update_axes(ax, mappable)" ] } ], "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.10.6" } }, "nbformat": 4, "nbformat_minor": 0 }