Skip to content

Commit 58d8f02

Browse files
AntSimiAntSimi
committed
Add time cube example with FSLE
1 parent 1b2207a commit 58d8f02

File tree

4 files changed

+303
-0
lines changed

4 files changed

+303
-0
lines changed

‎CHANGELOG.rst‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -34,6 +34,7 @@ Added
3434
- Network:
3535
- Add method to find relatives segments
3636
- Add method to get cloase network in an other atlas
37+
- Management of time cube data for advection
3738

3839
[3.3.0] - 2020-12-03
3940
--------------------

‎examples/07_cube_manipulation/README.rst‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,6 @@
1+
Time grid computation
2+
=====================
3+
4+
.. warning::
5+
6+
Time grid is under development, API could move quickly!

‎examples/07_cube_manipulation/pet_fsle_med.py‎

Lines changed: 152 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,152 @@
1+
"""
2+
FSLE experiment in med
3+
======================
4+
5+
Example to build FSLE, parameter values must be adapted for your case.
6+
7+
Example use a method similar to `AVISO flse`_
8+
9+
.. _AVISO flse:
10+
https://www.aviso.altimetry.fr/en/data/products/value-added-products/
11+
fsle-finite-size-lyapunov-exponents/fsle-description.html
12+
13+
"""
14+
15+
from matplotlib import pyplot as plt
16+
from numba import njit
17+
from numpy import arange, empty, isnan, log2, ma, meshgrid, zeros
18+
19+
from py_eddy_tracker import start_logger
20+
from py_eddy_tracker.data import get_path
21+
from py_eddy_tracker.dataset.grid import GridCollection, RegularGridDataset
22+
23+
start_logger().setLevel("ERROR")
24+
25+
26+
# %%
27+
# ADT in med
28+
# ---------------------
29+
c = GridCollection.from_netcdf_cube(
30+
get_path("dt_med_allsat_phy_l4_2005T2.nc"),
31+
"longitude",
32+
"latitude",
33+
"time",
34+
heigth="adt",
35+
)
36+
37+
38+
# %%
39+
# Methods to compute fsle
40+
# -----------------------
41+
@njit(cache=True, fastmath=True)
42+
def check_p(x, y, g, m, dt, dist_init=0.02, dist_max=0.6):
43+
"""
44+
Check if distance between eastern or northern particle to center particle is bigger than `dist_max`
45+
"""
46+
nb_p = x.shape[0] // 3
47+
delta = dist_max ** 2
48+
for i in range(nb_p):
49+
i0 = i * 3
50+
i_n = i0 + 1
51+
i_e = i0 + 2
52+
# If particle already set, we skip
53+
if m[i0] or m[i_n] or m[i_e]:
54+
continue
55+
# Distance with north
56+
dxn, dyn = x[i0] - x[i_n], y[i0] - y[i_n]
57+
dn = dxn ** 2 + dyn ** 2
58+
# Distance with east
59+
dxe, dye = x[i0] - x[i_e], y[i0] - y[i_e]
60+
de = dxe ** 2 + dye ** 2
61+
62+
if dn >= delta or de >= delta:
63+
s1 = dxe ** 2 + dxn ** 2 + dye ** 2 + dyn ** 2
64+
s2 = ((dxn + dye) ** 2 + (dxe - dyn) ** 2) * (
65+
(dxn - dye) ** 2 + (dxe + dyn) ** 2
66+
)
67+
g[i] = 1 / (2 * dt) * log2(1 / (2 * dist_init ** 2) * (s1 + s2 ** 0.5))
68+
m[i0], m[i_n], m[i_e] = True, True, True
69+
70+
71+
@njit(cache=True)
72+
def build_triplet(x, y, step=0.02):
73+
"""
74+
Triplet building for each position we add east and north point with defined step
75+
"""
76+
nb_x = x.shape[0]
77+
x_ = empty(nb_x * 3, dtype=x.dtype)
78+
y_ = empty(nb_x * 3, dtype=y.dtype)
79+
for i in range(nb_x):
80+
i0 = i * 3
81+
i_n, i_e = i0 + 1, i0 + 2
82+
x__, y__ = x[i], y[i]
83+
x_[i0], y_[i0] = x__, y__
84+
x_[i_n], y_[i_n] = x__, y__ + step
85+
x_[i_e], y_[i_e] = x__ + step, y__
86+
return x_, y_
87+
88+
89+
# %%
90+
# Particles
91+
# ---------
92+
step = 0.02
93+
t0 = 20268
94+
x0_, y0_ = -5, 30
95+
lon_p, lat_p = arange(x0_, x0_ + 43, step), arange(y0_, y0_ + 16, step)
96+
x0, y0 = meshgrid(lon_p, lat_p)
97+
grid_shape = x0.shape
98+
x0, y0 = x0.reshape(-1), y0.reshape(-1)
99+
# Identify all particle not on land
100+
m = ~isnan(c[t0].interp("adt", x0, y0))
101+
x0, y0 = x0[m], y0[m]
102+
103+
# %%
104+
# FSLE
105+
# ----
106+
time_step_by_days = 5
107+
# Array to compute fsle
108+
fsle = zeros(x0.shape[0], dtype="f4")
109+
x, y = build_triplet(x0, y0)
110+
used = zeros(x.shape[0], dtype="bool")
111+
112+
# advection generator
113+
kw = dict(t_init=t0, nb_step=1, backward=True, mask_particule=used)
114+
p = c.advect(x, y, "u", "v", time_step=86400 / time_step_by_days, **kw)
115+
116+
nb_days = 85
117+
# We check at each step of advection if particle distance is over `dist_max`
118+
for i in range(time_step_by_days * nb_days):
119+
t, xt, yt = p.__next__()
120+
dt = t / 86400.0 - t0
121+
check_p(xt, yt, fsle, used, dt, dist_max=0.2, dist_init=step)
122+
123+
# Get index with original_position
124+
i = ((x0 - x0_) / step).astype("i4")
125+
j = ((y0 - y0_) / step).astype("i4")
126+
fsle_ = empty(grid_shape, dtype="f4")
127+
used_ = zeros(grid_shape, dtype="bool")
128+
fsle_[j, i] = fsle
129+
used_[j, i] = used[::3]
130+
# Create a grid object
131+
fsle_custom = RegularGridDataset.with_array(
132+
coordinates=("lon", "lat"),
133+
datas=dict(
134+
fsle=ma.array(fsle_.T, mask=~used_.T),
135+
lon=lon_p,
136+
lat=lat_p,
137+
),
138+
centered=True,
139+
)
140+
141+
# %%
142+
# Display FSLE
143+
# ------------
144+
fig = plt.figure(figsize=(13, 5), dpi=150)
145+
ax = fig.add_axes([0.03, 0.03, 0.90, 0.94])
146+
ax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)
147+
ax.set_aspect("equal")
148+
ax.set_title("Finite size lyapunov exponent", weight="bold")
149+
kw = dict(cmap="viridis_r", vmin=-15, vmax=0)
150+
m = fsle_custom.display(ax, 1 / fsle_custom.grid("fsle"), **kw)
151+
ax.grid()
152+
cb = plt.colorbar(m, cax=fig.add_axes([0.94, 0.05, 0.01, 0.9]))

‎notebooks/python_module/07_cube_manipulation/pet_fsle_med.ipynb‎

Lines changed: 144 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,144 @@
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# FSLE experiment in med\n\nExample to build FSLE, parameter values must be adapted for your case.\n\nExample use a method similar to `AVISO flse`_\n\n https://www.aviso.altimetry.fr/en/data/products/value-added-products/\n fsle-finite-size-lyapunov-exponents/fsle-description.html\n"
19+
]
20+
},
21+
{
22+
"cell_type": "code",
23+
"execution_count": null,
24+
"metadata": {
25+
"collapsed": false
26+
},
27+
"outputs": [],
28+
"source": [
29+
"from matplotlib import pyplot as plt\nfrom numba import njit\nfrom numpy import arange, empty, isnan, log2, ma, meshgrid, zeros\n\nfrom py_eddy_tracker import start_logger\nfrom py_eddy_tracker.data import get_path\nfrom py_eddy_tracker.dataset.grid import GridCollection, RegularGridDataset\n\nstart_logger().setLevel(\"ERROR\")"
30+
]
31+
},
32+
{
33+
"cell_type": "markdown",
34+
"metadata": {},
35+
"source": [
36+
"## ADT in med\n\n"
37+
]
38+
},
39+
{
40+
"cell_type": "code",
41+
"execution_count": null,
42+
"metadata": {
43+
"collapsed": false
44+
},
45+
"outputs": [],
46+
"source": [
47+
"c = GridCollection.from_netcdf_cube(\n get_path(\"dt_med_allsat_phy_l4_2005T2.nc\"),\n \"longitude\",\n \"latitude\",\n \"time\",\n heigth=\"adt\",\n)"
48+
]
49+
},
50+
{
51+
"cell_type": "markdown",
52+
"metadata": {},
53+
"source": [
54+
"## Methods to compute fsle\n\n"
55+
]
56+
},
57+
{
58+
"cell_type": "code",
59+
"execution_count": null,
60+
"metadata": {
61+
"collapsed": false
62+
},
63+
"outputs": [],
64+
"source": [
65+
"@njit(cache=True, fastmath=True)\ndef check_p(x, y, g, m, dt, dist_init=0.02, dist_max=0.6):\n \"\"\"\n Check if distance between eastern or northern particle to center particle is bigger than `dist_max`\n \"\"\"\n nb_p = x.shape[0] // 3\n delta = dist_max ** 2\n for i in range(nb_p):\n i0 = i * 3\n i_n = i0 + 1\n i_e = i0 + 2\n # If particle already set, we skip\n if m[i0] or m[i_n] or m[i_e]:\n continue\n # Distance with north\n dxn, dyn = x[i0] - x[i_n], y[i0] - y[i_n]\n dn = dxn ** 2 + dyn ** 2\n # Distance with east\n dxe, dye = x[i0] - x[i_e], y[i0] - y[i_e]\n de = dxe ** 2 + dye ** 2\n\n if dn >= delta or de >= delta:\n s1 = dxe ** 2 + dxn ** 2 + dye ** 2 + dyn ** 2\n s2 = ((dxn + dye) ** 2 + (dxe - dyn) ** 2) * (\n (dxn - dye) ** 2 + (dxe + dyn) ** 2\n )\n g[i] = 1 / (2 * dt) * log2(1 / (2 * dist_init ** 2) * (s1 + s2 ** 0.5))\n m[i0], m[i_n], m[i_e] = True, True, True\n\n\n@njit(cache=True)\ndef build_triplet(x, y, step=0.02):\n \"\"\"\n Triplet building for each position we add east and north point with defined step\n \"\"\"\n nb_x = x.shape[0]\n x_ = empty(nb_x * 3, dtype=x.dtype)\n y_ = empty(nb_x * 3, dtype=y.dtype)\n for i in range(nb_x):\n i0 = i * 3\n i_n, i_e = i0 + 1, i0 + 2\n x__, y__ = x[i], y[i]\n x_[i0], y_[i0] = x__, y__\n x_[i_n], y_[i_n] = x__, y__ + step\n x_[i_e], y_[i_e] = x__ + step, y__\n return x_, y_"
66+
]
67+
},
68+
{
69+
"cell_type": "markdown",
70+
"metadata": {},
71+
"source": [
72+
"## Particles\n\n"
73+
]
74+
},
75+
{
76+
"cell_type": "code",
77+
"execution_count": null,
78+
"metadata": {
79+
"collapsed": false
80+
},
81+
"outputs": [],
82+
"source": [
83+
"step = 0.02\nt0 = 20268\nx0_, y0_ = -5, 30\nlon_p, lat_p = arange(x0_, x0_ + 43, step), arange(y0_, y0_ + 16, step)\nx0, y0 = meshgrid(lon_p, lat_p)\ngrid_shape = x0.shape\nx0, y0 = x0.reshape(-1), y0.reshape(-1)\n# Identify all particle not on land\nm = ~isnan(c[t0].interp(\"adt\", x0, y0))\nx0, y0 = x0[m], y0[m]"
84+
]
85+
},
86+
{
87+
"cell_type": "markdown",
88+
"metadata": {},
89+
"source": [
90+
"## FSLE\n\n"
91+
]
92+
},
93+
{
94+
"cell_type": "code",
95+
"execution_count": null,
96+
"metadata": {
97+
"collapsed": false
98+
},
99+
"outputs": [],
100+
"source": [
101+
"time_step_by_days = 5\n# Array to compute fsle\nfsle = zeros(x0.shape[0], dtype=\"f4\")\nx, y = build_triplet(x0, y0)\nused = zeros(x.shape[0], dtype=\"bool\")\n\n# advection generator\nkw = dict(t_init=t0, nb_step=1, backward=True, mask_particule=used)\np = c.advect(x, y, \"u\", \"v\", time_step=86400 / time_step_by_days, **kw)\n\nnb_days = 85\n# We check at each step of advection if particle distance is over `dist_max`\nfor i in range(time_step_by_days * nb_days):\n t, xt, yt = p.__next__()\n dt = t / 86400.0 - t0\n check_p(xt, yt, fsle, used, dt, dist_max=0.2, dist_init=step)\n\n# Get index with original_position\ni = ((x0 - x0_) / step).astype(\"i4\")\nj = ((y0 - y0_) / step).astype(\"i4\")\nfsle_ = empty(grid_shape, dtype=\"f4\")\nused_ = zeros(grid_shape, dtype=\"bool\")\nfsle_[j, i] = fsle\nused_[j, i] = used[::3]\n# Create a grid object\nfsle_custom = RegularGridDataset.with_array(\n coordinates=(\"lon\", \"lat\"),\n datas=dict(\n fsle=ma.array(fsle_.T, mask=~used_.T),\n lon=lon_p,\n lat=lat_p,\n ),\n centered=True,\n)"
102+
]
103+
},
104+
{
105+
"cell_type": "markdown",
106+
"metadata": {},
107+
"source": [
108+
"## Display FSLE\n\n"
109+
]
110+
},
111+
{
112+
"cell_type": "code",
113+
"execution_count": null,
114+
"metadata": {
115+
"collapsed": false
116+
},
117+
"outputs": [],
118+
"source": [
119+
"fig = plt.figure(figsize=(13, 5), dpi=150)\nax = fig.add_axes([0.03, 0.03, 0.90, 0.94])\nax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)\nax.set_aspect(\"equal\")\nax.set_title(\"Finite size lyapunov exponent\", weight=\"bold\")\nkw = dict(cmap=\"viridis_r\", vmin=-15, vmax=0)\nm = fsle_custom.display(ax, 1 / fsle_custom.grid(\"fsle\"), **kw)\nax.grid()\ncb = plt.colorbar(m, cax=fig.add_axes([0.94, 0.05, 0.01, 0.9]))"
120+
]
121+
}
122+
],
123+
"metadata": {
124+
"kernelspec": {
125+
"display_name": "Python 3",
126+
"language": "python",
127+
"name": "python3"
128+
},
129+
"language_info": {
130+
"codemirror_mode": {
131+
"name": "ipython",
132+
"version": 3
133+
},
134+
"file_extension": ".py",
135+
"mimetype": "text/x-python",
136+
"name": "python",
137+
"nbconvert_exporter": "python",
138+
"pygments_lexer": "ipython3",
139+
"version": "3.7.7"
140+
}
141+
},
142+
"nbformat": 4,
143+
"nbformat_minor": 0
144+
}

0 commit comments

Comments
 (0)