Skip to content

Commit aa82bbc

Browse files
AntSimiAntSimi
committed
Add an example about LAVD
1 parent 401f1ab commit aa82bbc

File tree

4 files changed

+429
-2
lines changed

4 files changed

+429
-2
lines changed

‎examples/06_grid_manipulation/pet_lavd.py‎

Lines changed: 174 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,174 @@
1+
"""
2+
LAVD experiment
3+
===============
4+
5+
Naive method to reproduce LAVD(Lagrangian-Averaged Vorticity deviation) method with a static velocity field.
6+
In the current example we didn't remove a mean vorticity.
7+
8+
Method are described here:
9+
10+
- `Transport by Lagrangian Vortices in the Eastern Pacific <https://doi.org/10.1175/JPO-D-17-0102.1>`_
11+
- `Transport by Coherent Lagrangian Vortices`_
12+
13+
.. _Transport by Coherent Lagrangian Vortices:
14+
https://usclivar.org/sites/default/files/meetings/2019/presentations/Aberernathey_CLIVAR.pdf
15+
16+
"""
17+
import re
18+
19+
from matplotlib import pyplot as plt
20+
from matplotlib.animation import FuncAnimation
21+
from numpy import arange, meshgrid, zeros
22+
23+
import py_eddy_tracker.gui
24+
from py_eddy_tracker.data import get_path
25+
from py_eddy_tracker.dataset.grid import RegularGridDataset
26+
from py_eddy_tracker.observations.network import NetworkObservations
27+
28+
29+
# %%
30+
def start_ax(title="", dpi=90):
31+
fig = plt.figure(figsize=(16, 9), dpi=dpi)
32+
ax = fig.add_axes([0, 0, 1, 1], projection="full_axes")
33+
ax.set_xlim(0, 32), ax.set_ylim(28, 46), ax.update_env()
34+
ax.set_title(title)
35+
return fig, ax, ax.text(3, 32, "", fontsize=20)
36+
37+
38+
def update_axes(ax, mappable=None):
39+
ax.grid()
40+
if mappable:
41+
cb = plt.colorbar(
42+
mappable,
43+
cax=ax.figure.add_axes([0.05, 0.1, 0.9, 0.01]),
44+
orientation="horizontal",
45+
)
46+
cb.set_label("Vorticity integration along trajectory at initial position")
47+
48+
49+
kw_vorticity = dict(vmin=0, vmax=2e-5, cmap="viridis")
50+
51+
52+
# %%
53+
class VideoAnimation(FuncAnimation):
54+
def _repr_html_(self, *args, **kwargs):
55+
"""To get video in html and have a player"""
56+
content = self.to_html5_video()
57+
return re.sub(
58+
r'width="[0-9]*"\sheight="[0-9]*"', 'width="100%" height="100%"', content
59+
)
60+
61+
return
62+
63+
def save(self, *args, **kwargs):
64+
if args[0].endswith("gif"):
65+
# In this case gif is use to create thumbnail which are not use but consume same time than video
66+
# So we create an empty file, to save time
67+
with open(args[0], "w") as _:
68+
pass
69+
return
70+
return super().save(*args, **kwargs)
71+
72+
73+
# %%
74+
# Data
75+
# ----
76+
# To compute vorticity(:math:`\omega`) we compute u/v field with a stencil and apply the following equation with stencil
77+
# method :
78+
#
79+
# .. math::
80+
# \omega = \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y}
81+
g = RegularGridDataset(
82+
get_path("dt_med_allsat_phy_l4_20160515_20190101.nc"), "longitude", "latitude"
83+
)
84+
g.add_uv("adt")
85+
u_y = g.compute_stencil(g.grid("u"), vertical=True)
86+
v_x = g.compute_stencil(g.grid("v"))
87+
g.vars["vort"] = v_x - u_y
88+
89+
# %%
90+
# Display vorticity field
91+
fig, ax, _ = start_ax()
92+
mappable = g.display(ax, abs(g.grid("vort")), **kw_vorticity)
93+
update_axes(ax, mappable)
94+
95+
# %%
96+
# Particles
97+
# ---------
98+
# Particles specification
99+
step = 1 / 32
100+
x_g, y_g = arange(0, 36, step), arange(28, 46, step)
101+
x, y = meshgrid(x_g, y_g)
102+
original_shape = x.shape
103+
x, y = x.reshape(-1), y.reshape(-1)
104+
print(f"{len(x)} particles advected")
105+
# A frame every 8h
106+
step_by_day = 3
107+
# Compute step of advection every 4h
108+
nb_step = 2
109+
kw_p = dict(nb_step=nb_step, time_step=86400 / step_by_day / nb_step)
110+
# Start a generator which at each iteration return new position at next time step
111+
particule = g.advect(x, y, "u", "v", **kw_p, rk4=True)
112+
113+
# %%
114+
# LAVD
115+
# ----
116+
lavd = zeros(original_shape)
117+
# Advection time
118+
nb_days = 8
119+
# Nb frame
120+
nb_time = step_by_day * nb_days
121+
i = 0.0
122+
123+
124+
# %%
125+
# Anim
126+
# ^^^^
127+
# Movie of LAVD integration at each integration time step.
128+
def update(i_frame):
129+
global lavd, i
130+
i += 1
131+
x, y = particule.__next__()
132+
# Interp vorticity on new_position
133+
lavd += abs(g.interp("vort", x, y).reshape(original_shape) * 1 / nb_time)
134+
txt.set_text(f"T0 + {i / step_by_day:.2f} days of advection")
135+
pcolormesh.set_array(lavd / i * nb_time)
136+
return pcolormesh, txt
137+
138+
139+
kw_video = dict(frames=arange(nb_time), interval=1000.0 / step_by_day / 2, blit=True)
140+
fig, ax, txt = start_ax(dpi=60)
141+
x_g_, y_g_ = arange(0 - step / 2, 36 + step / 2, step), arange(
142+
28 - step / 2, 46 + step / 2, step
143+
)
144+
# pcolorfast will be faster than pcolormesh, we could use pcolorfast due to x and y are regular
145+
pcolormesh = ax.pcolorfast(x_g_, y_g_, lavd, **kw_vorticity)
146+
update_axes(ax, pcolormesh)
147+
ani = VideoAnimation(ax.figure, update, **kw_video)
148+
149+
# %%
150+
# Final LAVD
151+
# ^^^^^^^^^^
152+
153+
# %%
154+
# Format LAVD data
155+
lavd = RegularGridDataset.with_array(
156+
coordinates=("lon", "lat"),
157+
datas=dict(
158+
lavd=lavd.T,
159+
lon=x_g,
160+
lat=y_g,
161+
),
162+
centered=True,
163+
)
164+
165+
# %%
166+
# Display final LAVD with py eddy tracker detection.
167+
# Period used for LAVD integration(8 days) is too short for a real use, but choose for example efficiency.
168+
fig, ax, _ = start_ax()
169+
mappable = lavd.display(ax, "lavd", **kw_vorticity)
170+
NetworkObservations.load_file(get_path("Anticyclonic_20160515.nc")).display(
171+
ax, color="k"
172+
)
173+
NetworkObservations.load_file(get_path("Cyclonic_20160515.nc")).display(ax, color="k")
174+
update_axes(ax, mappable)

‎notebooks/python_module/06_grid_manipulation/pet_lavd.ipynb‎

Lines changed: 209 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,209 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "code",
5+
"execution_count": null,
6+
"metadata": {
7+
"collapsed": false
8+
},
9+
"outputs": [],
10+
"source": [
11+
"%matplotlib inline"
12+
]
13+
},
14+
{
15+
"cell_type": "markdown",
16+
"metadata": {},
17+
"source": [
18+
"\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 - `Transport by Lagrangian Vortices in the Eastern Pacific <https://doi.org/10.1175/JPO-D-17-0102.1>`_\n - `Transport by Coherent Lagrangian Vortices`_\n\n https://usclivar.org/sites/default/files/meetings/2019/presentations/Aberernathey_CLIVAR.pdf\n"
19+
]
20+
},
21+
{
22+
"cell_type": "code",
23+
"execution_count": null,
24+
"metadata": {
25+
"collapsed": false
26+
},
27+
"outputs": [],
28+
"source": [
29+
"import re\n\nfrom matplotlib import pyplot as plt\nfrom matplotlib.animation import FuncAnimation\nfrom numpy import arange, meshgrid, zeros\n\nimport py_eddy_tracker.gui\nfrom py_eddy_tracker.data import get_path\nfrom py_eddy_tracker.dataset.grid import RegularGridDataset\nfrom py_eddy_tracker.observations.network import NetworkObservations"
30+
]
31+
},
32+
{
33+
"cell_type": "code",
34+
"execution_count": null,
35+
"metadata": {
36+
"collapsed": false
37+
},
38+
"outputs": [],
39+
"source": [
40+
"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=\"full_axes\")\n ax.set_xlim(0, 32), ax.set_ylim(28, 46), ax.update_env()\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\n\nkw_vorticity = dict(vmin=0, vmax=2e-5, cmap=\"viridis\")"
41+
]
42+
},
43+
{
44+
"cell_type": "code",
45+
"execution_count": null,
46+
"metadata": {
47+
"collapsed": false
48+
},
49+
"outputs": [],
50+
"source": [
51+
"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 return\n\n def save(self, *args, **kwargs):\n if args[0].endswith(\"gif\"):\n # In this case gif is use to create thumbnail which are not use 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)"
52+
]
53+
},
54+
{
55+
"cell_type": "markdown",
56+
"metadata": {},
57+
"source": [
58+
"## 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"
59+
]
60+
},
61+
{
62+
"cell_type": "code",
63+
"execution_count": null,
64+
"metadata": {
65+
"collapsed": false
66+
},
67+
"outputs": [],
68+
"source": [
69+
"g = RegularGridDataset(\n get_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"
70+
]
71+
},
72+
{
73+
"cell_type": "markdown",
74+
"metadata": {},
75+
"source": [
76+
"Display vorticity field\n\n"
77+
]
78+
},
79+
{
80+
"cell_type": "code",
81+
"execution_count": null,
82+
"metadata": {
83+
"collapsed": false
84+
},
85+
"outputs": [],
86+
"source": [
87+
"fig, ax, _ = start_ax()\nmappable = g.display(ax, abs(g.grid(\"vort\")), **kw_vorticity)\nupdate_axes(ax, mappable)"
88+
]
89+
},
90+
{
91+
"cell_type": "markdown",
92+
"metadata": {},
93+
"source": [
94+
"## Particles\nParticles specification\n\n"
95+
]
96+
},
97+
{
98+
"cell_type": "code",
99+
"execution_count": null,
100+
"metadata": {
101+
"collapsed": false
102+
},
103+
"outputs": [],
104+
"source": [
105+
"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)\n# Start a generator which at each iteration return new position at next time step\nparticule = g.advect(x, y, \"u\", \"v\", **kw_p, rk4=True)"
106+
]
107+
},
108+
{
109+
"cell_type": "markdown",
110+
"metadata": {},
111+
"source": [
112+
"## LAVD\n\n"
113+
]
114+
},
115+
{
116+
"cell_type": "code",
117+
"execution_count": null,
118+
"metadata": {
119+
"collapsed": false
120+
},
121+
"outputs": [],
122+
"source": [
123+
"lavd = zeros(original_shape)\n# Advection time\nnb_days = 8\n# Nb frame\nnb_time = step_by_day * nb_days\ni = 0.0"
124+
]
125+
},
126+
{
127+
"cell_type": "markdown",
128+
"metadata": {},
129+
"source": [
130+
"### Anim\nMovie of LAVD integration at each integration time step.\n\n"
131+
]
132+
},
133+
{
134+
"cell_type": "code",
135+
"execution_count": null,
136+
"metadata": {
137+
"collapsed": false
138+
},
139+
"outputs": [],
140+
"source": [
141+
"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_ = arange(0 - step / 2, 36 + step / 2, step), arange(\n 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)\nani = VideoAnimation(ax.figure, update, **kw_video)"
142+
]
143+
},
144+
{
145+
"cell_type": "markdown",
146+
"metadata": {},
147+
"source": [
148+
"### Final LAVD\n\n"
149+
]
150+
},
151+
{
152+
"cell_type": "markdown",
153+
"metadata": {},
154+
"source": [
155+
"Format LAVD data\n\n"
156+
]
157+
},
158+
{
159+
"cell_type": "code",
160+
"execution_count": null,
161+
"metadata": {
162+
"collapsed": false
163+
},
164+
"outputs": [],
165+
"source": [
166+
"lavd = RegularGridDataset.with_array(\n coordinates=(\"lon\", \"lat\"),\n datas=dict(\n lavd=lavd.T,\n lon=x_g,\n lat=y_g,\n ),\n centered=True,\n)"
167+
]
168+
},
169+
{
170+
"cell_type": "markdown",
171+
"metadata": {},
172+
"source": [
173+
"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"
174+
]
175+
},
176+
{
177+
"cell_type": "code",
178+
"execution_count": null,
179+
"metadata": {
180+
"collapsed": false
181+
},
182+
"outputs": [],
183+
"source": [
184+
"fig, ax, _ = start_ax()\nmappable = lavd.display(ax, \"lavd\", **kw_vorticity)\nNetworkObservations.load_file(get_path(\"Anticyclonic_20160515.nc\")).display(\n ax, color=\"k\"\n)\nNetworkObservations.load_file(get_path(\"Cyclonic_20160515.nc\")).display(ax, color=\"k\")\nupdate_axes(ax, mappable)"
185+
]
186+
}
187+
],
188+
"metadata": {
189+
"kernelspec": {
190+
"display_name": "Python 3",
191+
"language": "python",
192+
"name": "python3"
193+
},
194+
"language_info": {
195+
"codemirror_mode": {
196+
"name": "ipython",
197+
"version": 3
198+
},
199+
"file_extension": ".py",
200+
"mimetype": "text/x-python",
201+
"name": "python",
202+
"nbconvert_exporter": "python",
203+
"pygments_lexer": "ipython3",
204+
"version": "3.7.7"
205+
}
206+
},
207+
"nbformat": 4,
208+
"nbformat_minor": 0
209+
}

0 commit comments

Comments
 (0)