Skip to content

Commit 4149342

Browse files
AntSimiAntSimi
committed
Add example of simple advection
1 parent 97caffa commit 4149342

File tree

3 files changed

+331
-30
lines changed

3 files changed

+331
-30
lines changed

‎examples/06_grid_manipulation/pet_advect.py‎

Lines changed: 112 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,112 @@
1+
"""
2+
Grid advection
3+
==============
4+
5+
Dummy advection which use only static geostrophic current, which didn't resolve the complex circulation of the ocean.
6+
"""
7+
import numpy as np
8+
from matplotlib import pyplot as plt
9+
from matplotlib.animation import FuncAnimation
10+
11+
import py_eddy_tracker.gui
12+
from py_eddy_tracker import data
13+
from py_eddy_tracker.dataset.grid import RegularGridDataset
14+
from py_eddy_tracker.observations.observation import EddiesObservations
15+
16+
# %%
17+
# Load Input grid, ADT is used to detect eddies
18+
g = RegularGridDataset(
19+
data.get_path("dt_med_allsat_phy_l4_20160515_20190101.nc"), "longitude", "latitude"
20+
)
21+
# Compute u/v from height
22+
g.add_uv("adt")
23+
24+
# %%
25+
# Load detection files
26+
a = EddiesObservations.load_file(data.get_path("Anticyclonic_20160515.nc"))
27+
c = EddiesObservations.load_file(data.get_path("Cyclonic_20160515.nc"))
28+
29+
30+
# %%
31+
# Quiver from u/v with eddies
32+
fig = plt.figure(figsize=(10, 5))
33+
ax = fig.add_axes([0, 0, 1, 1], projection="full_axes")
34+
ax.set_xlim(19, 30), ax.set_ylim(31, 36.5), ax.grid()
35+
x, y = np.meshgrid(g.x_c, g.y_c)
36+
a.filled(ax, facecolors="r", alpha=0.1), c.filled(ax, facecolors="b", alpha=0.1)
37+
_ = ax.quiver(x.T, y.T, g.grid("u"), g.grid("v"), scale=20)
38+
39+
40+
# %%
41+
class VideoAnimation(FuncAnimation):
42+
def _repr_html_(self, *args, **kwargs):
43+
"""To get video in html and have a player"""
44+
return self.to_html5_video()
45+
46+
def save(self, *args, **kwargs):
47+
if args[0].endswith("gif"):
48+
# In this case gif is use to create thumbnail which are not use but consume same time than video
49+
# So we create an empty file, to save time
50+
with open(args[0], "w") as _:
51+
pass
52+
return
53+
return super().save(*args, **kwargs)
54+
55+
56+
# %%
57+
# Anim
58+
# ----
59+
# Particules positions
60+
x, y = np.meshgrid(np.arange(13, 36, 0.125), np.arange(28, 40, 0.125))
61+
x, y = x.reshape(-1), y.reshape(-1)
62+
# Remove all original position that we can't advect at first place
63+
m = ~np.isnan(g.interp("u", x, y))
64+
x, y = x[m], y[m]
65+
66+
# Movie properties
67+
kwargs = dict(frames=np.arange(51), interval=90)
68+
kw_p = dict(nb_step=2, time_step=21600)
69+
frame_t = kw_p["nb_step"] * kw_p["time_step"] / 86400.0
70+
71+
72+
def anim_ax(generator, **kw):
73+
t = 0
74+
for _ in range(4):
75+
generator.__next__()
76+
t += frame_t
77+
78+
fig = plt.figure(figsize=(10, 5), dpi=64)
79+
ax = fig.add_axes([0, 0, 1, 1], projection="full_axes")
80+
ax.set_xlim(19, 30), ax.set_ylim(31, 36.5), ax.grid()
81+
a.filled(ax, facecolors="r", alpha=0.1), c.filled(ax, facecolors="b", alpha=0.1)
82+
return fig, ax.text(21, 32.1, ""), ax.plot([], [], "k", **kw)[0], t
83+
84+
85+
def update(i_frame, t_step):
86+
global t
87+
x, y = p.__next__()
88+
t += t_step
89+
l.set_data(x, y)
90+
txt.set_text(f"T0 + {t:.1f} days")
91+
92+
93+
# %%
94+
# Filament forward
95+
# ^^^^^^^^^^^^^^^^
96+
p = g.filament(x, y, "u", "v", **kw_p, filament_size=4, rk4=True)
97+
fig, txt, l, t = anim_ax(p, lw=0.5)
98+
ani = VideoAnimation(fig, update, **kwargs, fargs=(frame_t,))
99+
100+
# %%
101+
# Filament backward
102+
# ^^^^^^^^^^^^^^^^^
103+
p = g.filament(x, y, "u", "v", **kw_p, filament_size=4, backward=True, rk4=True)
104+
fig, txt, l, t = anim_ax(p, lw=0.5)
105+
ani = VideoAnimation(fig, update, **kwargs, fargs=(-frame_t,))
106+
107+
# %%
108+
# Particule forward
109+
# ^^^^^^^^^^^^^^^^^
110+
p = g.advect(x, y, "u", "v", **kw_p, rk4=True)
111+
fig, txt, l, t = anim_ax(p, ls="", marker=".", markersize=1)
112+
ani = VideoAnimation(fig, update, **kwargs, fargs=(frame_t,))

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

Lines changed: 191 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,191 @@
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# Grid advection\n\nDummy advection which use only static geostrophic current, which didn't resolve the complex circulation of the ocean.\n"
19+
]
20+
},
21+
{
22+
"cell_type": "code",
23+
"execution_count": null,
24+
"metadata": {
25+
"collapsed": false
26+
},
27+
"outputs": [],
28+
"source": [
29+
"import numpy as np\nfrom matplotlib import pyplot as plt\nfrom matplotlib.animation import FuncAnimation\n\nimport py_eddy_tracker.gui\nfrom py_eddy_tracker import data\nfrom py_eddy_tracker.dataset.grid import RegularGridDataset\nfrom py_eddy_tracker.observations.observation import EddiesObservations"
30+
]
31+
},
32+
{
33+
"cell_type": "markdown",
34+
"metadata": {},
35+
"source": [
36+
"Load Input grid, ADT is used to detect eddies\n\n"
37+
]
38+
},
39+
{
40+
"cell_type": "code",
41+
"execution_count": null,
42+
"metadata": {
43+
"collapsed": false
44+
},
45+
"outputs": [],
46+
"source": [
47+
"g = RegularGridDataset(\n data.get_path(\"dt_med_allsat_phy_l4_20160515_20190101.nc\"), \"longitude\", \"latitude\"\n)\n# Compute u/v from height\ng.add_uv(\"adt\")"
48+
]
49+
},
50+
{
51+
"cell_type": "markdown",
52+
"metadata": {},
53+
"source": [
54+
"Load detection files\n\n"
55+
]
56+
},
57+
{
58+
"cell_type": "code",
59+
"execution_count": null,
60+
"metadata": {
61+
"collapsed": false
62+
},
63+
"outputs": [],
64+
"source": [
65+
"a = EddiesObservations.load_file(data.get_path(\"Anticyclonic_20160515.nc\"))\nc = EddiesObservations.load_file(data.get_path(\"Cyclonic_20160515.nc\"))"
66+
]
67+
},
68+
{
69+
"cell_type": "markdown",
70+
"metadata": {},
71+
"source": [
72+
"Quiver from u/v with eddies\n\n"
73+
]
74+
},
75+
{
76+
"cell_type": "code",
77+
"execution_count": null,
78+
"metadata": {
79+
"collapsed": false
80+
},
81+
"outputs": [],
82+
"source": [
83+
"fig = plt.figure(figsize=(10, 5))\nax = fig.add_axes([0, 0, 1, 1], projection=\"full_axes\")\nax.set_xlim(19, 30), ax.set_ylim(31, 36.5), ax.grid()\nx, y = np.meshgrid(g.x_c, g.y_c)\na.filled(ax, facecolors=\"r\", alpha=0.1), c.filled(ax, facecolors=\"b\", alpha=0.1)\n_ = ax.quiver(x.T, y.T, g.grid(\"u\"), g.grid(\"v\"), scale=20)"
84+
]
85+
},
86+
{
87+
"cell_type": "code",
88+
"execution_count": null,
89+
"metadata": {
90+
"collapsed": false
91+
},
92+
"outputs": [],
93+
"source": [
94+
"class VideoAnimation(FuncAnimation):\n def _repr_html_(self, *args, **kwargs):\n \"\"\"To get video in html and have a player\"\"\"\n return self.to_html5_video()\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)"
95+
]
96+
},
97+
{
98+
"cell_type": "markdown",
99+
"metadata": {},
100+
"source": [
101+
"## Anim\nParticules positions\n\n"
102+
]
103+
},
104+
{
105+
"cell_type": "code",
106+
"execution_count": null,
107+
"metadata": {
108+
"collapsed": false
109+
},
110+
"outputs": [],
111+
"source": [
112+
"x, y = np.meshgrid(np.arange(13, 36, 0.125), np.arange(28, 40, 0.125))\nx, y = x.reshape(-1), y.reshape(-1)\n# Remove all original position that we can't advect at first place\nm = ~np.isnan(g.interp(\"u\", x, y))\nx, y = x[m], y[m]\n\n# Movie properties\nkwargs = dict(frames=np.arange(51), interval=90)\nkw_p = dict(nb_step=2, time_step=21600)\nframe_t = kw_p[\"nb_step\"] * kw_p[\"time_step\"] / 86400.0\n\n\ndef anim_ax(generator, **kw):\n t = 0\n for _ in range(4):\n generator.__next__()\n t += frame_t\n\n fig = plt.figure(figsize=(10, 5), dpi=64)\n ax = fig.add_axes([0, 0, 1, 1], projection=\"full_axes\")\n ax.set_xlim(19, 30), ax.set_ylim(31, 36.5), ax.grid()\n a.filled(ax, facecolors=\"r\", alpha=0.1), c.filled(ax, facecolors=\"b\", alpha=0.1)\n return fig, ax.text(21, 32.1, \"\"), ax.plot([], [], \"k\", **kw)[0], t\n\n\ndef update(i_frame, t_step):\n global t\n x, y = p.__next__()\n t += t_step\n l.set_data(x, y)\n txt.set_text(f\"T0 + {t:.1f} days\")"
113+
]
114+
},
115+
{
116+
"cell_type": "markdown",
117+
"metadata": {},
118+
"source": [
119+
"### Filament forward\n\n"
120+
]
121+
},
122+
{
123+
"cell_type": "code",
124+
"execution_count": null,
125+
"metadata": {
126+
"collapsed": false
127+
},
128+
"outputs": [],
129+
"source": [
130+
"p = g.filament(x, y, \"u\", \"v\", **kw_p, filament_size=4, rk4=True)\nfig, txt, l, t = anim_ax(p, lw=0.5)\nani = VideoAnimation(fig, update, **kwargs, fargs=(frame_t,))"
131+
]
132+
},
133+
{
134+
"cell_type": "markdown",
135+
"metadata": {},
136+
"source": [
137+
"### Filament backward\n\n"
138+
]
139+
},
140+
{
141+
"cell_type": "code",
142+
"execution_count": null,
143+
"metadata": {
144+
"collapsed": false
145+
},
146+
"outputs": [],
147+
"source": [
148+
"p = g.filament(x, y, \"u\", \"v\", **kw_p, filament_size=4, backward=True, rk4=True)\nfig, txt, l, t = anim_ax(p, lw=0.5)\nani = VideoAnimation(fig, update, **kwargs, fargs=(-frame_t,))"
149+
]
150+
},
151+
{
152+
"cell_type": "markdown",
153+
"metadata": {},
154+
"source": [
155+
"### Particule forward\n\n"
156+
]
157+
},
158+
{
159+
"cell_type": "code",
160+
"execution_count": null,
161+
"metadata": {
162+
"collapsed": false
163+
},
164+
"outputs": [],
165+
"source": [
166+
"p = g.advect(x, y, \"u\", \"v\", **kw_p, rk4=True)\nfig, txt, l, t = anim_ax(p, ls=\"\", marker=\".\", markersize=1)\nani = VideoAnimation(fig, update, **kwargs, fargs=(frame_t,))"
167+
]
168+
}
169+
],
170+
"metadata": {
171+
"kernelspec": {
172+
"display_name": "Python 3",
173+
"language": "python",
174+
"name": "python3"
175+
},
176+
"language_info": {
177+
"codemirror_mode": {
178+
"name": "ipython",
179+
"version": 3
180+
},
181+
"file_extension": ".py",
182+
"mimetype": "text/x-python",
183+
"name": "python",
184+
"nbconvert_exporter": "python",
185+
"pygments_lexer": "ipython3",
186+
"version": "3.7.7"
187+
}
188+
},
189+
"nbformat": 4,
190+
"nbformat_minor": 0
191+
}

0 commit comments

Comments
 (0)