Skip to content

Commit b9a801e

Browse files
AntSimiAntSimi
authored
Merge pull request AntSimi#127 from AntSimi/raw_stats
example to compute statistics on raw identification
2 parents ab67c55 + cc540c0 commit b9a801e

File tree

6 files changed

+330
-5
lines changed

6 files changed

+330
-5
lines changed

‎README.md‎

Lines changed: 12 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -17,15 +17,17 @@ Method was described in :
1717
### Use case ###
1818

1919
Method is used in :
20-
20+
2121
[Mason, E., A. Pascual, P. Gaube, S.Ruiz, J. Pelegrí, A. Delepoulle, 2017: Subregional characterization of mesoscale eddies across the Brazil-Malvinas Confluence](https://doi.org/10.1002/2016JC012611)
2222

2323
### How do I get set up? ###
2424

2525
#### Short story ####
26+
2627
```bash
2728
pip install pyeddytracker
2829
```
30+
2931
#### Long story ####
3032

3133
To avoid problems with installation, use of the virtualenv Python virtual environment is recommended.
@@ -36,12 +38,20 @@ Then use pip to install all dependencies (numpy, scipy, matplotlib, netCDF4, ...
3638
pip install numpy scipy netCDF4 matplotlib opencv-python pyyaml pint polygon3
3739
```
3840

39-
Then run the following to install the eddy tracker:
41+
Clone :
42+
43+
```bash
44+
git clone https://github.com/AntSimi/py-eddy-tracker
45+
```
46+
47+
Then run the following to install the eddy tracker :
4048

4149
```bash
4250
python setup.py install
4351
```
52+
4453
### Tools gallery ###
54+
4555
Several examples based on py eddy tracker module are [here](https://py-eddy-tracker.readthedocs.io/en/latest/python_module/index.html).
4656

4757
[![](https://py-eddy-tracker.readthedocs.io/en/latest/_static/logo.png)](https://py-eddy-tracker.readthedocs.io/en/latest/python_module/index.html)

‎examples/02_eddy_identification/pet_statistics_on_identification.py‎

Lines changed: 105 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,105 @@
1+
"""
2+
Stastics on identification files
3+
================================
4+
5+
Some statistics on raw identification without any tracking
6+
"""
7+
import numpy as np
8+
from matplotlib import pyplot as plt
9+
from matplotlib.dates import date2num
10+
11+
from py_eddy_tracker import start_logger
12+
from py_eddy_tracker.data import get_remote_demo_sample
13+
from py_eddy_tracker.observations.observation import EddiesObservations
14+
15+
start_logger().setLevel("ERROR")
16+
17+
18+
# %%
19+
def start_axes(title):
20+
fig = plt.figure(figsize=(13, 5))
21+
ax = fig.add_axes([0.03, 0.03, 0.90, 0.94])
22+
ax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)
23+
ax.set_aspect("equal")
24+
ax.set_title(title)
25+
return ax
26+
27+
28+
def update_axes(ax, mappable=None):
29+
ax.grid()
30+
if mappable:
31+
plt.colorbar(mappable, cax=ax.figure.add_axes([0.95, 0.05, 0.01, 0.9]))
32+
33+
34+
# %%
35+
# We load demo sample and take only first year.
36+
#
37+
# Replace by a list of filename to apply on your own dataset.
38+
file_objects = get_remote_demo_sample(
39+
"eddies_med_adt_allsat_dt2018/Anticyclonic_2010_2011_2012"
40+
)[:365]
41+
42+
# %%
43+
# Merge all identification datasets in one object
44+
all_a = EddiesObservations.concatenate(
45+
[EddiesObservations.load_file(i) for i in file_objects]
46+
)
47+
48+
# %%
49+
# We define polygon bound
50+
x0, x1, y0, y1 = 15, 20, 33, 38
51+
xs = np.array([[x0, x1, x1, x0, x0]], dtype="f8")
52+
ys = np.array([[y0, y0, y1, y1, y0]], dtype="f8")
53+
# Polygon object created for the match function use.
54+
polygon = dict(contour_lon_e=xs, contour_lat_e=ys, contour_lon_s=xs, contour_lat_s=ys)
55+
56+
# %%
57+
# Geographic frequency of eddies
58+
step = 0.125
59+
ax = start_axes("")
60+
# Count pixel encompassed in each contour
61+
g_a = all_a.grid_count(bins=((-10, 37, step), (30, 46, step)), intern=True)
62+
m = g_a.display(
63+
ax, cmap="terrain_r", vmin=0, vmax=0.75, factor=1 / all_a.nb_days, name="count"
64+
)
65+
ax.plot(polygon["contour_lon_e"][0], polygon["contour_lat_e"][0], "r")
66+
update_axes(ax, m)
67+
68+
# %%
69+
# We use the match function to count the number of eddies that intersect the polygon defined previously
70+
# `p1_area` option allow to get in c_e/c_s output, precentage of area occupy by eddies in the polygon.
71+
i_e, j_e, c_e = all_a.match(polygon, p1_area=True, intern=False)
72+
i_s, j_s, c_s = all_a.match(polygon, p1_area=True, intern=True)
73+
74+
# %%
75+
dt = np.datetime64("1970-01-01") - np.datetime64("1950-01-01")
76+
kw_hist = dict(
77+
bins=date2num(np.arange(21900, 22300).astype("datetime64") - dt), histtype="step"
78+
)
79+
# translate julian day in datetime64
80+
t = all_a.time.astype("datetime64") - dt
81+
# %%
82+
# Number of eddies within a polygon
83+
ax = plt.figure(figsize=(12, 6)).add_subplot(111)
84+
ax.set_title("Different ways to count eddies within a polygon")
85+
ax.set_ylabel("Count")
86+
m = all_a.mask_from_polygons(((xs, ys),))
87+
ax.hist(t[m], label="Eddy Center in polygon", **kw_hist)
88+
ax.hist(t[i_s[c_s > 0]], label="Intersection Speed contour and polygon", **kw_hist)
89+
ax.hist(t[i_e[c_e > 0]], label="Intersection Effective contour and polygon", **kw_hist)
90+
ax.legend()
91+
ax.set_xlim(np.datetime64("2010"), np.datetime64("2011"))
92+
ax.grid()
93+
94+
# %%
95+
# Percent of the area of interest occupied by eddies.
96+
ax = plt.figure(figsize=(12, 6)).add_subplot(111)
97+
ax.set_title("Percent of polygon occupied by an anticyclonic eddy")
98+
ax.set_ylabel("Percent of polygon")
99+
ax.hist(t[i_s[c_s > 0]], weights=c_s[c_s > 0] * 100.0, label="speed contour", **kw_hist)
100+
ax.hist(
101+
t[i_e[c_e > 0]], weights=c_e[c_e > 0] * 100.0, label="effective contour", **kw_hist
102+
)
103+
ax.legend(), ax.set_ylim(0, 25)
104+
ax.set_xlim(np.datetime64("2010"), np.datetime64("2011"))
105+
ax.grid()

‎notebooks/python_module/02_eddy_identification/pet_statistics_on_identification.ipynb‎

Lines changed: 202 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,202 @@
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# Stastics on identification files\n\nSome statistics on raw identification without any tracking\n"
19+
]
20+
},
21+
{
22+
"cell_type": "code",
23+
"execution_count": null,
24+
"metadata": {
25+
"collapsed": false
26+
},
27+
"outputs": [],
28+
"source": [
29+
"import numpy as np\nfrom matplotlib import pyplot as plt\nfrom matplotlib.dates import date2num\n\nfrom py_eddy_tracker import start_logger\nfrom py_eddy_tracker.data import get_remote_demo_sample\nfrom py_eddy_tracker.observations.observation import EddiesObservations\n\nstart_logger().setLevel(\"ERROR\")"
30+
]
31+
},
32+
{
33+
"cell_type": "code",
34+
"execution_count": null,
35+
"metadata": {
36+
"collapsed": false
37+
},
38+
"outputs": [],
39+
"source": [
40+
"def start_axes(title):\n fig = plt.figure(figsize=(13, 5))\n ax = fig.add_axes([0.03, 0.03, 0.90, 0.94])\n ax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)\n ax.set_aspect(\"equal\")\n ax.set_title(title)\n return ax\n\n\ndef update_axes(ax, mappable=None):\n ax.grid()\n if mappable:\n plt.colorbar(mappable, cax=ax.figure.add_axes([0.95, 0.05, 0.01, 0.9]))"
41+
]
42+
},
43+
{
44+
"cell_type": "markdown",
45+
"metadata": {},
46+
"source": [
47+
"We load demo sample and take only first year.\n\nReplace by a list of filename to apply on your own dataset.\n\n"
48+
]
49+
},
50+
{
51+
"cell_type": "code",
52+
"execution_count": null,
53+
"metadata": {
54+
"collapsed": false
55+
},
56+
"outputs": [],
57+
"source": [
58+
"file_objects = get_remote_demo_sample(\n \"eddies_med_adt_allsat_dt2018/Anticyclonic_2010_2011_2012\"\n)[:365]"
59+
]
60+
},
61+
{
62+
"cell_type": "markdown",
63+
"metadata": {},
64+
"source": [
65+
"Merge all identification dataset in one object\n\n"
66+
]
67+
},
68+
{
69+
"cell_type": "code",
70+
"execution_count": null,
71+
"metadata": {
72+
"collapsed": false
73+
},
74+
"outputs": [],
75+
"source": [
76+
"all_a = EddiesObservations.concatenate(\n [EddiesObservations.load_file(i) for i in file_objects]\n)"
77+
]
78+
},
79+
{
80+
"cell_type": "markdown",
81+
"metadata": {},
82+
"source": [
83+
"We define polygon bound\n\n"
84+
]
85+
},
86+
{
87+
"cell_type": "code",
88+
"execution_count": null,
89+
"metadata": {
90+
"collapsed": false
91+
},
92+
"outputs": [],
93+
"source": [
94+
"x0, x1, y0, y1 = 15, 20, 33, 38\nxs = np.array([[x0, x1, x1, x0, x0]], dtype=\"f8\")\nys = np.array([[y0, y0, y1, y1, y0]], dtype=\"f8\")\n# Polygon object is create to be usable by match function.\npolygon = dict(contour_lon_e=xs, contour_lat_e=ys, contour_lon_s=xs, contour_lat_s=ys)"
95+
]
96+
},
97+
{
98+
"cell_type": "markdown",
99+
"metadata": {},
100+
"source": [
101+
"Geographic frequency of eddies\n\n"
102+
]
103+
},
104+
{
105+
"cell_type": "code",
106+
"execution_count": null,
107+
"metadata": {
108+
"collapsed": false
109+
},
110+
"outputs": [],
111+
"source": [
112+
"step = 0.125\nax = start_axes(\"\")\n# Count pixel used for each contour\ng_a = all_a.grid_count(bins=((-10, 37, step), (30, 46, step)), intern=True)\nm = g_a.display(\n ax, cmap=\"terrain_r\", vmin=0, vmax=0.75, factor=1 / all_a.nb_days, name=\"count\"\n)\nax.plot(polygon[\"contour_lon_e\"][0], polygon[\"contour_lat_e\"][0], \"r\")\nupdate_axes(ax, m)"
113+
]
114+
},
115+
{
116+
"cell_type": "markdown",
117+
"metadata": {},
118+
"source": [
119+
"We use match function to count number of eddies which intersect the polygon defined previously.\n`p1_area` option allow to get in c_e/c_s output, precentage of area occupy by eddies in the polygon.\n\n"
120+
]
121+
},
122+
{
123+
"cell_type": "code",
124+
"execution_count": null,
125+
"metadata": {
126+
"collapsed": false
127+
},
128+
"outputs": [],
129+
"source": [
130+
"i_e, j_e, c_e = all_a.match(polygon, p1_area=True, intern=False)\ni_s, j_s, c_s = all_a.match(polygon, p1_area=True, intern=True)"
131+
]
132+
},
133+
{
134+
"cell_type": "code",
135+
"execution_count": null,
136+
"metadata": {
137+
"collapsed": false
138+
},
139+
"outputs": [],
140+
"source": [
141+
"dt = np.datetime64(\"1970-01-01\") - np.datetime64(\"1950-01-01\")\nkw_hist = dict(\n bins=date2num(np.arange(21900, 22300).astype(\"datetime64\") - dt), histtype=\"step\"\n)\n# translate julian day in datetime64\nt = all_a.time.astype(\"datetime64\") - dt"
142+
]
143+
},
144+
{
145+
"cell_type": "markdown",
146+
"metadata": {},
147+
"source": [
148+
"Count how many are in polygon\n\n"
149+
]
150+
},
151+
{
152+
"cell_type": "code",
153+
"execution_count": null,
154+
"metadata": {
155+
"collapsed": false
156+
},
157+
"outputs": [],
158+
"source": [
159+
"ax = plt.figure(figsize=(12, 6)).add_subplot(111)\nax.set_title(\"Different way to count eddies presence in a polygon\")\nax.set_ylabel(\"Count\")\nm = all_a.mask_from_polygons(((xs, ys),))\nax.hist(t[m], label=\"center in polygon\", **kw_hist)\nax.hist(t[i_s[c_s > 0]], label=\"intersect speed contour with polygon\", **kw_hist)\nax.hist(t[i_e[c_e > 0]], label=\"intersect extern contour with polygon\", **kw_hist)\nax.legend()\nax.set_xlim(np.datetime64(\"2010\"), np.datetime64(\"2011\"))\nax.grid()"
160+
]
161+
},
162+
{
163+
"cell_type": "markdown",
164+
"metadata": {},
165+
"source": [
166+
"Percent of are of interest occupy by eddies\n\n"
167+
]
168+
},
169+
{
170+
"cell_type": "code",
171+
"execution_count": null,
172+
"metadata": {
173+
"collapsed": false
174+
},
175+
"outputs": [],
176+
"source": [
177+
"ax = plt.figure(figsize=(12, 6)).add_subplot(111)\nax.set_title(\"Percent of polygon occupy by an anticyclonic eddy\")\nax.set_ylabel(\"Percent of polygon\")\nax.hist(t[i_s[c_s > 0]], weights=c_s[c_s > 0] * 100.0, label=\"speed contour\", **kw_hist)\nax.hist(t[i_e[c_e > 0]], weights=c_e[c_e > 0] * 100.0, label=\"effective contour\", **kw_hist)\nax.legend(), ax.set_ylim(0, 25)\nax.set_xlim(np.datetime64(\"2010\"), np.datetime64(\"2011\"))\nax.grid()"
178+
]
179+
}
180+
],
181+
"metadata": {
182+
"kernelspec": {
183+
"display_name": "Python 3",
184+
"language": "python",
185+
"name": "python3"
186+
},
187+
"language_info": {
188+
"codemirror_mode": {
189+
"name": "ipython",
190+
"version": 3
191+
},
192+
"file_extension": ".py",
193+
"mimetype": "text/x-python",
194+
"name": "python",
195+
"nbconvert_exporter": "python",
196+
"pygments_lexer": "ipython3",
197+
"version": "3.7.7"
198+
}
199+
},
200+
"nbformat": 4,
201+
"nbformat_minor": 0
202+
}

‎src/py_eddy_tracker/dataset/grid.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -807,7 +807,7 @@ def eddy_identification(
807807
else:
808808
centi = reset_centroid[0]
809809
centj = reset_centroid[1]
810-
# To move in regular and unregular grid
810+
# FIXME : To move in regular and unregular grid
811811
if len(x.shape) == 1:
812812
centlon_e = x[centi]
813813
centlat_e = y[centj]

‎src/py_eddy_tracker/observations/observation.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@
88
from tarfile import ExFileObject
99
from tokenize import TokenError
1010

11-
import packaging
11+
import packaging.version
1212
import zarr
1313
from matplotlib.cm import get_cmap
1414
from matplotlib.collections import PolyCollection

0 commit comments

Comments
 (0)