py-eddy-tracker/examples/07_cube_manipulation/pet_cube.py at master · Yuxuan-PO/py-eddy-tracker · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
"""
Time advection
==============

Example which use CMEMS surface current with a Runge-Kutta 4 algorithm to advect particles.
"""
# sphinx_gallery_thumbnail_number = 2
import re
from datetime import datetime, timedelta

from matplotlib import pyplot as plt
from matplotlib.animation import FuncAnimation
from numpy import arange, isnan, meshgrid, ones

import py_eddy_tracker.gui
from py_eddy_tracker import start_logger
from py_eddy_tracker.data import get_path
from py_eddy_tracker.dataset.grid import GridCollection

start_logger().setLevel("ERROR")


# %%
class VideoAnimation(FuncAnimation):
    def _repr_html_(self, *args, **kwargs):
        """To get video in html and have a player"""
        content = self.to_html5_video()
        return re.sub(
            r'width="[0-9]*"\sheight="[0-9]*"', 'width="100%" height="100%"', content
        )

    def save(self, *args, **kwargs):
        if args[0].endswith("gif"):
            # In this case gif is use to create thumbnail which are not use but consume same time than video
            # So we create an empty file, to save time
            with open(args[0], "w") as _:
                pass
            return
        return super().save(*args, **kwargs)


# %%
# Data
# ----
# Load Input time grid ADT
c = GridCollection.from_netcdf_cube(
    get_path("dt_med_allsat_phy_l4_2005T2.nc"),
    "longitude",
    "latitude",
    "time",
    # To create U/V variable
    heigth="adt",
)

# %%
# Anim
# ----
# Particles setup
step_p = 1 / 8
x, y = meshgrid(arange(13, 36, step_p), arange(28, 40, step_p))
x, y = x.reshape(-1), y.reshape(-1)
# Remove all original position that we can't advect at first place
t0 = 20181
m = ~isnan(c[t0].interp("u", x, y))
x0, y0 = x[m], y[m]
x, y = x0.copy(), y0.copy()


# %%
# Function
def anim_ax(**kw):
    fig = plt.figure(figsize=(10, 5), dpi=55)
    axes = fig.add_axes([0, 0, 1, 1], projection="full_axes")
    axes.set_xlim(19, 30), axes.set_ylim(31, 36.5), axes.grid()
    line = axes.plot([], [], "k", **kw)[0]
    return fig, axes.text(21, 32.1, ""), line


def update(_):
    tt, xt, yt = f.__next__()
    mappable.set_data(xt, yt)
    d = timedelta(tt / 86400.0) + datetime(1950, 1, 1)
    txt.set_text(f"{d:%Y/%m/%d-%H}")


# %%
f = c.filament(x, y, "u", "v", t_init=t0, nb_step=2, time_step=21600, filament_size=3)
fig, txt, mappable = anim_ax(lw=0.5)
ani = VideoAnimation(fig, update, frames=arange(160), interval=100)


# %%
# Particules stat
# ---------------
# Time_step settings
# ^^^^^^^^^^^^^^^^^^
# Dummy experiment to test advection precision, we run particles 50 days forward and backward with different time step
# and we measure distance between new positions and original positions.
fig = plt.figure()
ax = fig.add_subplot(111)
kw = dict(
    bins=arange(0, 50, 0.002),
    cumulative=True,
    weights=ones(x0.shape) / x0.shape[0] * 100.0,
    histtype="step",
)
kw_p = dict(u_name="u", v_name="v", nb_step=1)
for time_step in (10800, 21600, 43200, 86400):
    x, y = x0.copy(), y0.copy()
    nb = int(30 * 86400 / time_step)
    # Go forward
    p = c.advect(x, y, time_step=time_step, t_init=20181.5, **kw_p)
    for i in range(nb):
        t_, _, _ = p.__next__()
    # Go backward
    p = c.advect(x, y, time_step=time_step, backward=True, t_init=t_ / 86400.0, **kw_p)
    for i in range(nb):
        t_, _, _ = p.__next__()
    d = ((x - x0) ** 2 + (y - y0) ** 2) ** 0.5
    ax.hist(d, **kw, label=f"{86400. / time_step:.0f} time step by day")
ax.set_xlim(0, 0.25), ax.set_ylim(0, 100), ax.legend(loc="lower right"), ax.grid()
ax.set_title("Distance after 50 days forward and 50 days backward")
ax.set_xlabel("Distance between original position and final position (in degrees)")
_ = ax.set_ylabel("Percent of particles with distance lesser than")

# %%
# Time duration
# ^^^^^^^^^^^^^
# We keep same time_step but change time duration
fig = plt.figure()
ax = fig.add_subplot(111)
time_step = 10800
for duration in (10, 40, 80):
    x, y = x0.copy(), y0.copy()
    nb = int(duration * 86400 / time_step)
    # Go forward
    p = c.advect(x, y, time_step=time_step, t_init=20181.5, **kw_p)
    for i in range(nb):
        t_, _, _ = p.__next__()
    # Go backward
    p = c.advect(x, y, time_step=time_step, backward=True, t_init=t_ / 86400.0, **kw_p)
    for i in range(nb):
        t_, _, _ = p.__next__()
    d = ((x - x0) ** 2 + (y - y0) ** 2) ** 0.5
    ax.hist(d, **kw, label=f"Time duration {duration} days")
ax.set_xlim(0, 0.25), ax.set_ylim(0, 100), ax.legend(loc="lower right"), ax.grid()
ax.set_title(
    "Distance after N days forward and N days backward\nwith a time step of 1/8 days"
)
ax.set_xlabel("Distance between original position and final position (in degrees)")
_ = ax.set_ylabel("Percent of particles with distance lesser than ")