Skip to content

Commit eebbf58

Browse files
AntSimiAntSimi
committed
Add examples to display groups dependance at parameter
1 parent 7ce7bf1 commit eebbf58

File tree

4 files changed

+288
-4
lines changed

4 files changed

+288
-4
lines changed

‎CHANGELOG.md‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
88
## [Unreleased]
99

1010
### Added
11+
- Add a method to merge several indexs type for eddy obs
1112
- Acces at dataset variable like attribute, and lifetime/age are available for all observations
1213
- Add **EddyInfos** application to get general information about eddies dataset
1314
- Add method to inspect contour rejection (which are not in eddies)

‎examples/10_tracking_diagnostics/pet_groups.py‎

Lines changed: 79 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,79 @@
1+
"""
2+
Groups distribution
3+
==================
4+
5+
"""
6+
import py_eddy_tracker_sample
7+
from matplotlib import pyplot as plt
8+
from numpy import arange, ones, percentile
9+
10+
from py_eddy_tracker.observations.tracking import TrackEddiesObservations
11+
12+
# %%
13+
# Load an experimental med atlas over a period of 26 years (1993-2019)
14+
a = TrackEddiesObservations.load_file(
15+
py_eddy_tracker_sample.get_path("eddies_med_adt_allsat_dt2018/Anticyclonic.zarr")
16+
)
17+
18+
# %%
19+
# Group distribution
20+
groups = dict()
21+
bins_time = [10, 20, 30, 60, 90, 180, 360, 100000]
22+
for t0, t1 in zip(bins_time[:-1], bins_time[1:]):
23+
groups[f"lifetime_{t0}_{t1}"] = lambda dataset, t0=t0, t1=t1: (
24+
dataset.lifetime >= t0
25+
) * (dataset.lifetime < t1)
26+
bins_percentile = arange(0, 100.0001, 5)
27+
28+
29+
# %%
30+
# Function to build stats
31+
def stats_compilation(dataset, groups, field, bins, filter=None):
32+
datas = dict(ref=dataset.bins_stat(field, bins=bins, mask=filter)[1], y=dict())
33+
for k, index in groups.items():
34+
i = dataset.merge_indexs(filter, index)
35+
x, datas["y"][k] = dataset.bins_stat(field, bins=bins, mask=i)
36+
datas["x"], datas["bins"] = x, bins
37+
return datas
38+
39+
40+
def plot_stats(ax, bins, x, y, ref, box=False, cmap=None, percentiles=None, **kw):
41+
base, ref = ones(x.shape) * 100.0, ref / 100.0
42+
x = arange(bins.shape[0]).repeat(2)[1:-1] if box else x
43+
y0 = base
44+
if cmap is not None:
45+
cmap, nb_groups = plt.get_cmap(cmap), len(y)
46+
keys = tuple(y.keys())
47+
for i, k in enumerate(keys[::-1]):
48+
y1 = y0 - y[k] / ref
49+
args = (y0.repeat(2), y1.repeat(2)) if box else (y0, y1)
50+
if cmap is not None:
51+
kw["color"] = cmap(1 - i / (nb_groups - 1))
52+
ax.fill_between(x, *args, label=k, **kw)
53+
y0 = y1
54+
if percentiles:
55+
for b in bins:
56+
ax.axvline(b, **percentiles)
57+
58+
59+
# %%
60+
# Speed radius by track period
61+
stats = stats_compilation(
62+
a, groups, "radius_s", percentile(a.radius_s, bins_percentile)
63+
)
64+
fig = plt.figure()
65+
ax = fig.add_subplot(111)
66+
plot_stats(ax, **stats, cmap="magma", percentiles=dict(color="gray", ls="-.", lw=0.4))
67+
ax.set_xlabel("Speed radius (m)"), ax.set_ylabel("% of class"), ax.set_ylim(0, 100)
68+
ax.grid(), ax.legend()
69+
70+
# %%
71+
# Amplitude by track period
72+
stats = stats_compilation(
73+
a, groups, "amplitude", percentile(a.amplitude, bins_percentile)
74+
)
75+
fig = plt.figure()
76+
ax = fig.add_subplot(111)
77+
plot_stats(ax, **stats, cmap="magma")
78+
ax.set_xlabel("Amplitude (m)"), ax.set_ylabel("% of class"), ax.set_ylim(0, 100)
79+
ax.grid(), ax.legend()

‎notebooks/python_module/10_tracking_diagnostics/pet_groups.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# Groups distribution\n"
19+
]
20+
},
21+
{
22+
"cell_type": "code",
23+
"execution_count": null,
24+
"metadata": {
25+
"collapsed": false
26+
},
27+
"outputs": [],
28+
"source": [
29+
"import py_eddy_tracker_sample\nfrom matplotlib import pyplot as plt\nfrom numpy import arange, ones, percentile\n\nfrom py_eddy_tracker.observations.tracking import TrackEddiesObservations"
30+
]
31+
},
32+
{
33+
"cell_type": "markdown",
34+
"metadata": {},
35+
"source": [
36+
"Load an experimental med atlas over a period of 26 years (1993-2019)\n\n"
37+
]
38+
},
39+
{
40+
"cell_type": "code",
41+
"execution_count": null,
42+
"metadata": {
43+
"collapsed": false
44+
},
45+
"outputs": [],
46+
"source": [
47+
"a = TrackEddiesObservations.load_file(\n py_eddy_tracker_sample.get_path(\"eddies_med_adt_allsat_dt2018/Anticyclonic.zarr\")\n)"
48+
]
49+
},
50+
{
51+
"cell_type": "markdown",
52+
"metadata": {},
53+
"source": [
54+
"Group distribution\n\n"
55+
]
56+
},
57+
{
58+
"cell_type": "code",
59+
"execution_count": null,
60+
"metadata": {
61+
"collapsed": false
62+
},
63+
"outputs": [],
64+
"source": [
65+
"groups = dict()\nbins_time = [10, 20, 30, 60, 90, 180, 360, 100000]\nfor t0, t1 in zip(bins_time[:-1], bins_time[1:]):\n groups[f\"lifetime_{t0}_{t1}\"] = lambda dataset, t0=t0, t1=t1: (\n dataset.lifetime >= t0\n ) * (dataset.lifetime < t1)\nbins_percentile = arange(0, 100.0001, 5)"
66+
]
67+
},
68+
{
69+
"cell_type": "markdown",
70+
"metadata": {},
71+
"source": [
72+
"Function to build stats\n\n"
73+
]
74+
},
75+
{
76+
"cell_type": "code",
77+
"execution_count": null,
78+
"metadata": {
79+
"collapsed": false
80+
},
81+
"outputs": [],
82+
"source": [
83+
"def stats_compilation(dataset, groups, field, bins, filter=None):\n datas = dict(ref=dataset.bins_stat(field, bins=bins, mask=filter)[1], y=dict())\n for k, index in groups.items():\n i = dataset.merge_indexs(filter, index)\n x, datas[\"y\"][k] = dataset.bins_stat(field, bins=bins, mask=i)\n datas[\"x\"], datas[\"bins\"] = x, bins\n return datas\n\n\ndef plot_stats(ax, bins, x, y, ref, box=False, cmap=None, percentiles=None, **kw):\n base, ref = ones(x.shape) * 100.0, ref / 100.0\n x = arange(bins.shape[0]).repeat(2)[1:-1] if box else x\n y0 = base\n if cmap is not None:\n cmap, nb_groups = plt.get_cmap(cmap), len(y)\n keys = tuple(y.keys())\n for i, k in enumerate(keys[::-1]):\n y1 = y0 - y[k] / ref\n args = (y0.repeat(2), y1.repeat(2)) if box else (y0, y1)\n if cmap is not None:\n kw[\"color\"] = cmap(1 - i / (nb_groups - 1))\n ax.fill_between(x, *args, label=k, **kw)\n y0 = y1\n if percentiles:\n for b in bins:\n ax.axvline(b, **percentiles)"
84+
]
85+
},
86+
{
87+
"cell_type": "markdown",
88+
"metadata": {},
89+
"source": [
90+
"Speed radius by track period\n\n"
91+
]
92+
},
93+
{
94+
"cell_type": "code",
95+
"execution_count": null,
96+
"metadata": {
97+
"collapsed": false
98+
},
99+
"outputs": [],
100+
"source": [
101+
"stats = stats_compilation(\n a, groups, \"radius_s\", percentile(a.radius_s, bins_percentile)\n)\nfig = plt.figure()\nax = fig.add_subplot(111)\nplot_stats(ax, **stats, cmap=\"magma\", percentiles=dict(color=\"gray\", ls=\"-.\", lw=0.4))\nax.set_xlabel(\"Speed radius (m)\"), ax.set_ylabel(\"% of class\"), ax.set_ylim(0, 100)\nax.grid(), ax.legend()"
102+
]
103+
},
104+
{
105+
"cell_type": "markdown",
106+
"metadata": {},
107+
"source": [
108+
"Amplitude by track period\n\n"
109+
]
110+
},
111+
{
112+
"cell_type": "code",
113+
"execution_count": null,
114+
"metadata": {
115+
"collapsed": false
116+
},
117+
"outputs": [],
118+
"source": [
119+
"stats = stats_compilation(\n a, groups, \"amplitude\", percentile(a.amplitude, bins_percentile)\n)\nfig = plt.figure()\nax = fig.add_subplot(111)\nplot_stats(ax, **stats, cmap=\"magma\")\nax.set_xlabel(\"Amplitude (m)\"), ax.set_ylabel(\"% of class\"), ax.set_ylim(0, 100)\nax.grid(), ax.legend()"
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+
}

‎src/py_eddy_tracker/observations/observation.py‎

Lines changed: 64 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -28,6 +28,7 @@
2828
floor,
2929
histogram,
3030
histogram2d,
31+
in1d,
3132
isnan,
3233
linspace,
3334
ma,
@@ -1685,6 +1686,66 @@ def filled(
16851686
c.norm = Normalize(vmin=vmin, vmax=vmax)
16861687
return c
16871688

1689+
def merge_indexs(self, filter, index=None):
1690+
"""
1691+
Compute an intersectin between indexs and filters after to evaluate each of them
1692+
1693+
:param callable,None,slice,array[int],array[bool] filter:
1694+
:param callable,None,slice,array[int],array[bool] index:
1695+
1696+
:return: Return applicable object to numpy.array
1697+
:rtype: slice, index, mask
1698+
"""
1699+
# If filter is a function we apply on dataset
1700+
if callable(filter):
1701+
filter = filter(self)
1702+
# Solve indexs case
1703+
if index is not None:
1704+
index = self.merge_indexs(index)
1705+
# Merge indexs and filter
1706+
if index is None and filter is None:
1707+
return slice(None)
1708+
if index is None:
1709+
return filter
1710+
if filter is None:
1711+
return index
1712+
if isinstance(index, slice):
1713+
reject = ones(len(self), dtype="bool")
1714+
reject[index] = False
1715+
if isinstance(filter, slice):
1716+
reject[filter] = False
1717+
return ~reject
1718+
# Mask case
1719+
elif filter.dtype == bool:
1720+
return ~reject * filter
1721+
# index case
1722+
else:
1723+
return filter[~reject[filter]]
1724+
# mask case
1725+
elif index.dtype == bool:
1726+
if isinstance(filter, slice):
1727+
select = zeros(len(self), dtype="bool")
1728+
select[filter] = True
1729+
return select * index
1730+
# Mask case
1731+
elif filter.dtype == bool:
1732+
return filter * index
1733+
# index case
1734+
else:
1735+
return filter[index[filter]]
1736+
# index case
1737+
else:
1738+
if isinstance(filter, slice):
1739+
select = zeros(len(self), dtype="bool")
1740+
select[filter] = True
1741+
return index[select[index]]
1742+
# Mask case
1743+
elif filter.dtype == bool:
1744+
return index[filter[index]]
1745+
# index case
1746+
else:
1747+
return index[in1d(index, filter)]
1748+
16881749
def bins_stat(self, xname, bins=None, yname=None, method=None, mask=None):
16891750
"""
16901751
:param str,array xname: variable to compute stats on
@@ -1698,16 +1759,15 @@ def bins_stat(self, xname, bins=None, yname=None, method=None, mask=None):
16981759
.. minigallery:: py_eddy_tracker.EddiesObservations.bins_stat
16991760
"""
17001761
v = self[xname] if isinstance(xname, str) else xname
1701-
if mask is not None:
1702-
v = v[mask]
1762+
mask = self.merge_indexs(mask)
1763+
v = v[mask]
17031764
if bins is None:
17041765
bins = arange(v.min(), v.max() + 2)
17051766
y, x = hist_numba(v, bins=bins)
17061767
x = (x[1:] + x[:-1]) / 2
17071768
if method == "mean":
17081769
y_v = self[yname] if isinstance(yname, str) else yname
1709-
if mask is not None:
1710-
y_v = y_v[mask]
1770+
y_v = y_v[mask]
17111771
y_, _ = histogram(v, bins=bins, weights=y_v)
17121772
y = y_ / y
17131773
return x, y

0 commit comments

Comments
 (0)