py-eddy-tracker/examples/10_tracking_diagnostics/pet_center_count.py at ROMSgrid · HelloKaito/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
"""
Count center
============

Do Geo stat with center and compare with frequency method
show: :ref:`sphx_glr_python_module_10_tracking_diagnostics_pet_pixel_used.py`
"""
import py_eddy_tracker_sample
from matplotlib import pyplot as plt
from matplotlib.colors import LogNorm

from py_eddy_tracker.observations.tracking import TrackEddiesObservations

# %%
# Load an experimental med atlas over a period of 26 years (1993-2019)
a = TrackEddiesObservations.load_file(
    py_eddy_tracker_sample.get_demo_path(
        "eddies_med_adt_allsat_dt2018/Anticyclonic.zarr"
    )
)
c = TrackEddiesObservations.load_file(
    py_eddy_tracker_sample.get_demo_path("eddies_med_adt_allsat_dt2018/Cyclonic.zarr")
)

# %%
# Parameters
step = 0.125
bins = ((-10, 37, step), (30, 46, step))
kwargs_pcolormesh = dict(
    cmap="terrain_r", vmin=0, vmax=2, factor=1 / (a.nb_days * step ** 2), name="count"
)


# %%
# Plot
fig = plt.figure(figsize=(12, 18.5))
ax_a = fig.add_axes([0.03, 0.75, 0.90, 0.25])
ax_a.set_title("Anticyclonic center frequency")
ax_c = fig.add_axes([0.03, 0.5, 0.90, 0.25])
ax_c.set_title("Cyclonic center frequency")
ax_all = fig.add_axes([0.03, 0.25, 0.90, 0.25])
ax_all.set_title("All eddies center frequency")
ax_ratio = fig.add_axes([0.03, 0.0, 0.90, 0.25])
ax_ratio.set_title("Ratio cyclonic / Anticyclonic")

# Count pixel used for each center
g_a = a.grid_count(bins, intern=True, center=True)
g_a.display(ax_a, **kwargs_pcolormesh)
g_c = c.grid_count(bins, intern=True, center=True)
g_c.display(ax_c, **kwargs_pcolormesh)
# Compute a ratio Cyclonic / Anticyclonic
ratio = g_c.vars["count"] / g_a.vars["count"]

# Mask manipulation to be able to sum the 2 grids
m_c = g_c.vars["count"].mask
m = m_c & g_a.vars["count"].mask
g_c.vars["count"][m_c] = 0
g_c.vars["count"] += g_a.vars["count"]
g_c.vars["count"].mask = m

m = g_c.display(ax_all, **kwargs_pcolormesh)
cb = plt.colorbar(m, cax=fig.add_axes([0.94, 0.27, 0.01, 0.7]))
cb.set_label("Eddies by 1°^2 by day")

g_c.vars["count"] = ratio
m = g_c.display(
    ax_ratio, name="count", norm=LogNorm(vmin=0.1, vmax=10), cmap="coolwarm_r"
)
plt.colorbar(m, cax=fig.add_axes([0.94, 0.02, 0.01, 0.2]))

for ax in (ax_a, ax_c, ax_all, ax_ratio):
    ax.set_aspect("equal")
    ax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)
    ax.grid()

# %%
# Count Anticyclones as a function of lifetime
# --------------------------------------------
# Count at the center's position

fig = plt.figure(figsize=(12, 10))
mask = a.lifetime >= 60
ax_long = fig.add_axes([0.03, 0.53, 0.90, 0.45])
g_a = a.grid_count(bins, center=True, filter=mask)
g_a.display(ax_long, **kwargs_pcolormesh)
ax_long.set_title(f"Anticyclones with lifetime >= 60 days ({mask.sum()} Obs)")
ax_short = fig.add_axes([0.03, 0.03, 0.90, 0.45])
g_a = a.grid_count(bins, center=True, filter=~mask)
m = g_a.display(ax_short, **kwargs_pcolormesh)
ax_short.set_title(f"Anticyclones with lifetime < 60 days ({(~mask).sum()} Obs)")
for ax in (ax_short, ax_long):
    ax.set_aspect("equal"), ax.grid()
    ax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)
cb = plt.colorbar(m, cax=fig.add_axes([0.94, 0.05, 0.015, 0.9]))