Skip to content

Commit 50885c2

Browse files
committed
Add example to compare radius and area
1 parent 23cc80b commit 50885c2

File tree

5 files changed

+174
-8
lines changed

5 files changed

+174
-8
lines changed

examples/02_eddy_identification/pet_filter_and_detection.py

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -58,8 +58,8 @@ def update_axes(ax, mappable=None):
5858
# Parameters distribution
5959
# -----------------------
6060
fig = plt.figure(figsize=(12, 5))
61-
ax_a = plt.subplot(121, xlabel="amplitdue(cm)")
62-
ax_r = plt.subplot(122, xlabel="speed radius (km)")
61+
ax_a = fig.add_subplot(121, xlabel="amplitdue(cm)")
62+
ax_r = fig.add_subplot(122, xlabel="speed radius (km)")
6363
ax_a.hist(
6464
merge_f["amplitude"] * 100,
6565
bins=arange(0.0005, 100, 1),
Lines changed: 49 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,49 @@
1+
"""
2+
Radius vs area
3+
==============
4+
5+
"""
6+
from matplotlib import pyplot as plt
7+
from py_eddy_tracker.observations.observation import EddiesObservations
8+
from py_eddy_tracker import data
9+
10+
from py_eddy_tracker.poly import poly_area
11+
from py_eddy_tracker.generic import coordinates_to_local
12+
from numpy import pi, array
13+
14+
15+
# %%
16+
# Load detection files
17+
a = EddiesObservations.load_file(data.get_path("Anticyclonic_20190223.nc"))
18+
areas = list()
19+
# For each contour area will be compute in local reference
20+
for i in a:
21+
x, y = coordinates_to_local(
22+
i["contour_lon_s"], i["contour_lat_s"], i["lon"], i["lat"]
23+
)
24+
areas.append(poly_area(x, y))
25+
areas = array(areas)
26+
27+
# %%
28+
# Radius provided by eddy detection is computed with :func:`~py_eddy_tracker.poly.fit_circle` method.
29+
# This radius will be compared with an equivalent radius deduced from polygon area.
30+
ax = plt.subplot(111)
31+
ax.set_aspect("equal")
32+
ax.grid()
33+
ax.set_xlabel("Speed radius computed with fit_circle")
34+
ax.set_ylabel("Radius deduced from area\nof contour_lon_s/contour_lat_s")
35+
ax.set_title('Area vs radius')
36+
ax.plot(a["radius_s"] / 1000.0, (areas / pi) ** 0.5 / 1000.0, ".")
37+
ax.plot((0, 250), (0, 250), "r")
38+
39+
# %%
40+
# Fit circle give a radius bigger than polygon area
41+
42+
# %%
43+
# When error is tiny, radius are very close.
44+
ax = plt.subplot(111)
45+
ax.grid()
46+
ax.set_xlabel("Radius ratio")
47+
ax.set_ylabel("Shape error")
48+
ax.set_title("err = f(radius_ratio)")
49+
ax.plot(a["radius_s"] / (areas / pi) ** .5 , a['shape_error_s'], ".")

notebooks/python_module/02_eddy_identification/pet_filter_and_detection.ipynb

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,7 @@
1515
"cell_type": "markdown",
1616
"metadata": {},
1717
"source": [
18-
"\nEddy detection and filter\n=========================\n"
18+
"\n# Eddy detection and filter\n"
1919
]
2020
},
2121
{
@@ -98,7 +98,7 @@
9898
"cell_type": "markdown",
9999
"metadata": {},
100100
"source": [
101-
"Parameters distribution\n-----------------------\n\n"
101+
"## Parameters distribution\n\n"
102102
]
103103
},
104104
{
@@ -109,14 +109,14 @@
109109
},
110110
"outputs": [],
111111
"source": [
112-
"fig = plt.figure(figsize=(12, 5))\nax_a = plt.subplot(121, xlabel=\"amplitdue(cm)\")\nax_r = plt.subplot(122, xlabel=\"speed radius (km)\")\nax_a.hist(\n merge_f[\"amplitude\"] * 100,\n bins=arange(0.0005, 100, 1),\n label=\"Eddy from filtered grid\",\n histtype=\"step\",\n)\nax_a.hist(\n merge_r[\"amplitude\"] * 100,\n bins=arange(0.0005, 100, 1),\n label=\"Eddy from raw grid\",\n histtype=\"step\",\n)\nax_a.set_xlim(0, 10)\nax_r.hist(merge_f[\"radius_s\"] / 1000.0, bins=arange(0, 300, 5), histtype=\"step\")\nax_r.hist(merge_r[\"radius_s\"] / 1000.0, bins=arange(0, 300, 5), histtype=\"step\")\nax_r.set_xlim(0, 100)\nax_a.legend()"
112+
"fig = plt.figure(figsize=(12, 5))\nax_a = fig.add_subplot(121, xlabel=\"amplitdue(cm)\")\nax_r = fig.add_subplot(122, xlabel=\"speed radius (km)\")\nax_a.hist(\n merge_f[\"amplitude\"] * 100,\n bins=arange(0.0005, 100, 1),\n label=\"Eddy from filtered grid\",\n histtype=\"step\",\n)\nax_a.hist(\n merge_r[\"amplitude\"] * 100,\n bins=arange(0.0005, 100, 1),\n label=\"Eddy from raw grid\",\n histtype=\"step\",\n)\nax_a.set_xlim(0, 10)\nax_r.hist(merge_f[\"radius_s\"] / 1000.0, bins=arange(0, 300, 5), histtype=\"step\")\nax_r.hist(merge_r[\"radius_s\"] / 1000.0, bins=arange(0, 300, 5), histtype=\"step\")\nax_r.set_xlim(0, 100)\nax_a.legend()"
113113
]
114114
},
115115
{
116116
"cell_type": "markdown",
117117
"metadata": {},
118118
"source": [
119-
"Match detection and compare\n---------------------------\n\n"
119+
"## Match detection and compare\n\n"
120120
]
121121
},
122122
{
@@ -127,7 +127,7 @@
127127
},
128128
"outputs": [],
129129
"source": [
130-
"i_, j_, c = merge_f.match(merge_r, 0.1)"
130+
"i_, j_, c = merge_f.match(merge_r, cmin=0.1)"
131131
]
132132
},
133133
{
@@ -176,7 +176,7 @@
176176
"name": "python",
177177
"nbconvert_exporter": "python",
178178
"pygments_lexer": "ipython3",
179-
"version": "3.7.7"
179+
"version": "3.7.0"
180180
}
181181
},
182182
"nbformat": 4,
Lines changed: 115 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,115 @@
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# Radius vs area\n"
19+
]
20+
},
21+
{
22+
"cell_type": "code",
23+
"execution_count": null,
24+
"metadata": {
25+
"collapsed": false
26+
},
27+
"outputs": [],
28+
"source": [
29+
"from matplotlib import pyplot as plt\nfrom py_eddy_tracker.observations.observation import EddiesObservations\nfrom py_eddy_tracker import data\n\nfrom py_eddy_tracker.poly import poly_area\nfrom py_eddy_tracker.generic import coordinates_to_local\nfrom numpy import pi, array"
30+
]
31+
},
32+
{
33+
"cell_type": "markdown",
34+
"metadata": {},
35+
"source": [
36+
"Load detection files\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 = EddiesObservations.load_file(data.get_path(\"Anticyclonic_20190223.nc\"))\nareas = list()\n# For each contour area will be compute in local reference\nfor i in a:\n x, y = coordinates_to_local(\n i[\"contour_lon_s\"], i[\"contour_lat_s\"], i[\"lon\"], i[\"lat\"]\n )\n areas.append(poly_area(x, y))\nareas = array(areas)"
48+
]
49+
},
50+
{
51+
"cell_type": "markdown",
52+
"metadata": {},
53+
"source": [
54+
"Radius provided by eddy detection is computed with :func:`~py_eddy_tracker.poly.fit_circle` method.\nThis radius will be compared with an equivalent radius deduced from polygon area.\n\n"
55+
]
56+
},
57+
{
58+
"cell_type": "code",
59+
"execution_count": null,
60+
"metadata": {
61+
"collapsed": false
62+
},
63+
"outputs": [],
64+
"source": [
65+
"ax = plt.subplot(111)\nax.set_aspect(\"equal\")\nax.grid()\nax.set_xlabel(\"Speed radius computed with fit_circle\")\nax.set_ylabel(\"Radius deduced from area\\nof contour_lon_s/contour_lat_s\")\nax.set_title('Area vs radius')\nax.plot(a[\"radius_s\"] / 1000.0, (areas / pi) ** 0.5 / 1000.0, \".\")\nax.plot((0, 250), (0, 250), \"r\")"
66+
]
67+
},
68+
{
69+
"cell_type": "markdown",
70+
"metadata": {},
71+
"source": [
72+
"Fit circle give a radius bigger than polygon area\n\n"
73+
]
74+
},
75+
{
76+
"cell_type": "markdown",
77+
"metadata": {},
78+
"source": [
79+
"When error is tiny, radius are very close.\n\n"
80+
]
81+
},
82+
{
83+
"cell_type": "code",
84+
"execution_count": null,
85+
"metadata": {
86+
"collapsed": false
87+
},
88+
"outputs": [],
89+
"source": [
90+
"ax = plt.subplot(111)\nax.grid()\nax.set_xlabel(\"Radius ratio\")\nax.set_ylabel(\"Shape error\")\nax.set_title(\"err = f(radius_ratio)\")\nax.plot(a[\"radius_s\"] / (areas / pi) ** .5 , a['shape_error_s'], \".\")"
91+
]
92+
}
93+
],
94+
"metadata": {
95+
"kernelspec": {
96+
"display_name": "Python 3",
97+
"language": "python",
98+
"name": "python3"
99+
},
100+
"language_info": {
101+
"codemirror_mode": {
102+
"name": "ipython",
103+
"version": 3
104+
},
105+
"file_extension": ".py",
106+
"mimetype": "text/x-python",
107+
"name": "python",
108+
"nbconvert_exporter": "python",
109+
"pygments_lexer": "ipython3",
110+
"version": "3.7.0"
111+
}
112+
},
113+
"nbformat": 4,
114+
"nbformat_minor": 0
115+
}

src/py_eddy_tracker/poly.py

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -86,6 +86,8 @@ def poly_area_vertice(v):
8686
@njit(cache=True)
8787
def poly_area(x, y):
8888
"""
89+
Must be call with local coordinates (in m, to get an area in m²)
90+
8991
:param array x:
9092
:param array y:
9193
:return: area of polygon in coordinates unit

0 commit comments

Comments
 (0)