Skip to content

Commit 7977d6d

Browse files
AntSimiAntSimi
committed
Add new example for network
1 parent 9eea120 commit 7977d6d

File tree

2 files changed

+318
-0
lines changed

2 files changed

+318
-0
lines changed

‎examples/16_network/pet_atlas.py‎

Lines changed: 144 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,144 @@
1+
"""
2+
Network analyse
3+
===============
4+
"""
5+
from matplotlib import pyplot as plt
6+
from numpy import ma
7+
8+
import py_eddy_tracker.gui
9+
from py_eddy_tracker.observations.network import NetworkObservations
10+
11+
n = NetworkObservations.load_file(
12+
"/data/adelepoulle/work/Eddies/20201217_network_build/tracking/med/Anticyclonic_seg.nc"
13+
)
14+
15+
# %%
16+
# Parameters
17+
step = 1 / 10.0
18+
bins = ((-10, 37, step), (30, 46, step))
19+
kw_time = dict(cmap="terrain_r", factor=100.0 / n.nb_days, name="count")
20+
kw_ratio = dict(cmap=plt.get_cmap("coolwarm_r", 10))
21+
22+
# %%
23+
# Function
24+
def start_axes(title):
25+
fig = plt.figure(figsize=(13, 5))
26+
ax = fig.add_axes([0.03, 0.03, 0.90, 0.94], projection="full_axes")
27+
ax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)
28+
ax.set_aspect("equal")
29+
ax.set_title(title, weight="bold")
30+
return ax
31+
32+
def update_axes(ax, mappable=None):
33+
ax.grid()
34+
ax.update_env()
35+
if mappable:
36+
return plt.colorbar(mappable, cax=ax.figure.add_axes([0.94, 0.05, 0.01, 0.9]))
37+
38+
# %%
39+
# All
40+
# ---
41+
ax = start_axes("")
42+
g_all = n.grid_count(bins)
43+
m = g_all.display(ax, **kw_time, vmin=0, vmax=75)
44+
update_axes(ax, m).set_label("Pixel used in % of time")
45+
46+
# %%
47+
# Network longer than 10 days
48+
# ---------------------------
49+
ax = start_axes("")
50+
n10 = n.longer_than(10)
51+
g_10 = n10.grid_count(bins)
52+
m = g_10.display(ax, **kw_time, vmin=0, vmax=75)
53+
update_axes(ax, m).set_label("Pixel used in % of time")
54+
55+
ax = start_axes("")
56+
m = g_10.display(
57+
ax,
58+
**kw_ratio,
59+
vmin=50,
60+
vmax=100,
61+
name=g_10.vars["count"] * 100.0 / g_all.vars["count"]
62+
)
63+
update_axes(ax, m).set_label("Pixel used in % all atlas")
64+
# %%
65+
# Network longer than 20 days
66+
# ---------------------------
67+
ax = start_axes("")
68+
n20 = n.longer_than(20)
69+
g_20 = n20.grid_count(bins)
70+
m = g_20.display(ax, **kw_time, vmin=0, vmax=75)
71+
update_axes(ax, m).set_label("Pixel used in % of time")
72+
73+
ax = start_axes("")
74+
m = g_20.display(
75+
ax,
76+
**kw_ratio,
77+
vmin=50,
78+
vmax=100,
79+
name=g_20.vars["count"] * 100.0 / g_all.vars["count"]
80+
)
81+
update_axes(ax, m).set_label("Pixel used in % all atlas")
82+
ax = start_axes("")
83+
m = g_20.display(
84+
ax,
85+
**kw_ratio,
86+
vmin=50,
87+
vmax=100,
88+
name=ma.array(
89+
g_20.vars["count"] * 100.0 / g_all.vars["count"], mask=g_all.vars["count"] < 365
90+
)
91+
)
92+
update_axes(ax, m).set_label("Pixel used in % all atlas")
93+
ax = start_axes("")
94+
m = g_20.display(
95+
ax,
96+
**kw_ratio,
97+
vmin=50,
98+
vmax=100,
99+
name=ma.array(
100+
g_20.vars["count"] * 100.0 / g_all.vars["count"],
101+
mask=g_all.vars["count"] >= 365,
102+
)
103+
)
104+
update_axes(ax, m).set_label("Pixel used in % all atlas")
105+
106+
107+
# %%
108+
# All merging
109+
# -----------
110+
ax = start_axes("")
111+
g_all_merging = n.merging_event().grid_count(bins)
112+
m = g_all_merging.display(ax, **kw_time, vmin=0, vmax=0.2)
113+
update_axes(ax, m).set_label("Pixel used in % of time")
114+
# %%
115+
ax = start_axes("")
116+
m = g_all_merging.display(
117+
ax,
118+
**kw_ratio,
119+
vmin=0,
120+
vmax=1,
121+
name=g_all_merging.vars["count"] * 100.0 / g_all.vars["count"]
122+
)
123+
update_axes(ax, m).set_label("Pixel used in % all atlas")
124+
125+
# %%
126+
# Merging in network longer than 10 days
127+
# --------------------------------------
128+
ax = start_axes("")
129+
g_10_merging = n10.merging_event().grid_count(bins)
130+
m = g_10_merging.display(ax, **kw_time, vmin=0, vmax=0.2)
131+
update_axes(ax, m).set_label("Pixel used in % of time")
132+
# %%
133+
ax = start_axes("")
134+
m = g_10_merging.display(
135+
ax,
136+
**kw_ratio,
137+
vmin=0,
138+
vmax=0.5,
139+
name=ma.array(
140+
g_10_merging.vars["count"] * 100.0 / g_10.vars["count"],
141+
mask=g_10.vars["count"] < 365,
142+
)
143+
)
144+
update_axes(ax, m).set_label("Pixel used in % all atlas")

‎examples/16_network/pet_ioannou_2017_case.py‎

Lines changed: 174 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,174 @@
1+
"""
2+
Ioannou case
3+
============
4+
Figure 10 from https://doi.org/10.1002/2017JC013158
5+
6+
"""
7+
from datetime import datetime, timedelta
8+
9+
import numpy as np
10+
from matplotlib import pyplot as plt
11+
from matplotlib.animation import FuncAnimation
12+
from matplotlib.ticker import FuncFormatter
13+
14+
import py_eddy_tracker.gui
15+
from py_eddy_tracker.appli.gui import Anim
16+
from py_eddy_tracker.observations.network import NetworkObservations
17+
18+
19+
# %%
20+
class VideoAnimation(FuncAnimation):
21+
def _repr_html_(self, *args, **kwargs):
22+
return self.to_html5_video()
23+
24+
def save(self, *args, **kwargs):
25+
if args[0].endswith("gif"):
26+
return
27+
return super().save(*args, **kwargs)
28+
29+
30+
@FuncFormatter
31+
def formatter(x, pos):
32+
return (timedelta(x) + datetime(1950, 1, 1)).strftime("%d/%m/%Y")
33+
34+
35+
def start_axes(title=""):
36+
fig = plt.figure(figsize=(13, 6))
37+
ax = fig.add_axes([0.03, 0.03, 0.90, 0.94], projection="full_axes")
38+
ax.set_xlim(19, 29), ax.set_ylim(31, 35.5)
39+
ax.set_aspect("equal")
40+
ax.set_title(title, weight="bold")
41+
ax.update_env()
42+
return ax
43+
44+
45+
def timeline_axes(title=""):
46+
fig = plt.figure(figsize=(15, 5))
47+
ax = fig.add_axes([0.03, 0.06, 0.90, 0.88])
48+
ax.set_title(title, weight="bold")
49+
ax.xaxis.set_major_formatter(formatter), ax.grid()
50+
return ax
51+
52+
53+
def update_axes(ax, mappable=None):
54+
ax.grid(True)
55+
if mappable:
56+
return plt.colorbar(mappable, cax=ax.figure.add_axes([0.94, 0.05, 0.01, 0.9]))
57+
58+
59+
# %%
60+
n = NetworkObservations.load_file(
61+
"med/Anticyclonic_seg.nc"
62+
)
63+
i = np.where(
64+
(n.lat > 33)
65+
* (n.lat < 34)
66+
* (n.lon > 22)
67+
* (n.lon < 23)
68+
* (n.time > 20630)
69+
* (n.time < 20650)
70+
)[0][0]
71+
ioannou_case = n.extract_with_mask(n.track == n.track[i])
72+
print(ioannou_case.infos())
73+
74+
# %%
75+
ax = start_axes()
76+
ioannou_case.plot(ax)
77+
update_axes(ax)
78+
79+
# %%
80+
# Full Timeline
81+
# -------------
82+
fig = plt.figure(figsize=(15, 5))
83+
ax = fig.add_axes([0.04, 0.05, 0.92, 0.92])
84+
ax.xaxis.set_major_formatter(formatter), ax.grid()
85+
_ = ioannou_case.display_timeline(ax)
86+
87+
88+
# %%
89+
# Sub network and new numbering
90+
# -----------------------------
91+
i = np.where(
92+
(ioannou_case.lat > 33)
93+
* (ioannou_case.lat < 34)
94+
* (ioannou_case.lon > 22)
95+
* (ioannou_case.lon < 23)
96+
* (ioannou_case.time > 20630)
97+
* (ioannou_case.time < 20650)
98+
)[0][0]
99+
close_to_i3 = ioannou_case.relative(i, order=3)
100+
close_to_i3.numbering_segment()
101+
102+
# %%
103+
# Anim
104+
# ----
105+
a = Anim(
106+
close_to_i3,
107+
figsize=(12, 4),
108+
cmap="Spectral_r",
109+
nb_step=7,
110+
dpi=55,
111+
field_color="segment",
112+
field_txt="segment",
113+
)
114+
a.ax.set_xlim(19, 30), a.ax.set_ylim(32, 35.25)
115+
a.ax.update_env()
116+
a.txt.set_position((21.5, 32.7))
117+
kwargs = dict(frames=np.arange(*a.period), interval=100, blit=True)
118+
ani = VideoAnimation(a.fig, a.func_animation, **kwargs)
119+
120+
# %%
121+
# Classic display
122+
# ---------------
123+
ax = timeline_axes()
124+
_ = close_to_i3.display_timeline(ax)
125+
126+
# %%
127+
ax = start_axes("")
128+
n_copy = close_to_i3.copy()
129+
n_copy.position_filter(2, 4)
130+
n_copy.plot(ax)
131+
update_axes(ax)
132+
133+
# %%
134+
# Local Timeline
135+
# --------------
136+
ax = timeline_axes(f"Close segments ({close_to_i3.infos()})")
137+
n_copy = close_to_i3.copy()
138+
n_copy.median_filter(15, "time", "latitude")
139+
_ = n_copy.display_timeline(ax, field="lat", method="all")
140+
141+
# %%
142+
# Local radius timeline
143+
# ---------------------
144+
n_copy.median_filter(2, "time", "radius_e")
145+
n_copy.median_filter(2, "time", "radius_s")
146+
for b0, b1 in [
147+
(datetime(i, 1, 1), datetime(i, 12, 31)) for i in (2004, 2005, 2006, 2007)
148+
]:
149+
b0_, b1_ = (b0 - datetime(1950, 1, 1)).days - 10, (
150+
b1 - datetime(1950, 1, 1)
151+
).days + 10
152+
ax = timeline_axes()
153+
ax.set_xlim(b0_, b1_)
154+
ax.set_ylim(10, 115)
155+
n_copy.display_timeline(
156+
ax, field="radius_e", method="all", lw=4, markersize=8, factor=1e-3
157+
)
158+
n_copy.display_timeline(
159+
ax, field="radius_s", method="all", lw=1, markersize=3, factor=1e-3
160+
)
161+
162+
# %%
163+
# Parameter timeline
164+
# ------------------
165+
kw = dict(s=35, cmap=plt.get_cmap("Spectral_r", 8), zorder=10)
166+
ax = timeline_axes()
167+
m = close_to_i3.scatter_timeline(ax, "radius_e", factor=1e-3, vmin=20, vmax=100, **kw)
168+
cb = update_axes(ax, m["scatter"])
169+
cb.set_label("Effective radius (km)")
170+
# %%
171+
ax = timeline_axes()
172+
m = close_to_i3.scatter_timeline(ax, "shape_error_e", vmin=14, vmax=70, **kw)
173+
cb = update_axes(ax, m["scatter"])
174+
cb.set_label("Effective shape error")

0 commit comments

Comments
 (0)