Skip to content

Commit d96823c

Browse files
AntSimiAntSimi
committed
Add an example about filtering needs for detection
1 parent 146c578 commit d96823c

File tree

2 files changed

+261
-0
lines changed

2 files changed

+261
-0
lines changed

‎examples/02_eddy_identification/pet_filter_and_detection.py‎

Lines changed: 106 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,106 @@
1+
"""
2+
Eddy detection and filter
3+
=========================
4+
5+
"""
6+
from datetime import datetime
7+
from matplotlib import pyplot as plt
8+
from py_eddy_tracker.dataset.grid import RegularGridDataset
9+
from py_eddy_tracker import data
10+
from numpy import arange
11+
12+
13+
# %%
14+
def start_axes(title):
15+
fig = plt.figure(figsize=(13, 5))
16+
ax = fig.add_axes([0.03, 0.03, 0.90, 0.94])
17+
ax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)
18+
ax.set_aspect("equal")
19+
ax.set_title(title)
20+
return ax
21+
22+
23+
def update_axes(ax, mappable=None):
24+
ax.grid()
25+
if mappable:
26+
plt.colorbar(mappable, cax=ax.figure.add_axes([0.95, 0.05, 0.01, 0.9]))
27+
28+
29+
# %%
30+
# Load Input grid, ADT will be used to detect eddies
31+
32+
g = RegularGridDataset(
33+
data.get_path("dt_med_allsat_phy_l4_20160515_20190101.nc"), "longitude", "latitude",
34+
)
35+
g.add_uv("adt")
36+
g.copy("adt", "adt_high")
37+
wavelength = 400
38+
g.bessel_high_filter("adt_high", wavelength, order=3)
39+
date = datetime(2016, 5, 15)
40+
41+
# %%
42+
# Run algorithm of detection
43+
a_f, c_f = g.eddy_identification("adt_high", "u", "v", date, 0.002)
44+
merge_f = a_f.merge(c_f)
45+
a_r, c_r = g.eddy_identification("adt", "u", "v", date, 0.002)
46+
merge_r = a_r.merge(c_r)
47+
48+
# %%
49+
# Display detection
50+
ax = start_axes("Eddies detected over ADT")
51+
m = g.display(ax, "adt", vmin=-0.15, vmax=0.15)
52+
merge_f.display(ax, lw=0.5, label="Eddy from filtered grid", ref=-10, color="k")
53+
merge_r.display(ax, lw=0.5, label="Eddy from raw grid", ref=-10, color="r")
54+
ax.legend()
55+
update_axes(ax, m)
56+
57+
# %%
58+
# Parameters distribution
59+
# -----------------------
60+
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)")
63+
ax_a.grid(), ax_r.grid()
64+
ax_a.hist(
65+
merge_f["amplitude"] * 100,
66+
bins=arange(0.0005, 100, 1),
67+
label="Eddy from filtered grid",
68+
histtype="step",
69+
)
70+
ax_a.hist(
71+
merge_r["amplitude"] * 100,
72+
bins=arange(0.0005, 100, 1),
73+
label="Eddy from raw grid",
74+
histtype="step",
75+
)
76+
ax_a.set_xlim(0, 10)
77+
ax_r.hist(merge_f["radius_s"] / 1000.0, bins=arange(0, 300, 5), histtype="step")
78+
ax_r.hist(merge_r["radius_s"] / 1000.0, bins=arange(0, 300, 5), histtype="step")
79+
ax_r.set_xlim(0, 100)
80+
ax_a.legend()
81+
82+
# %%
83+
# Match detection and compare
84+
# ---------------------------
85+
86+
i, j, c = merge_f.match(merge_r)
87+
m = c > 0.1
88+
i_, j_ = i[m], j[m]
89+
fig = plt.figure(figsize=(12, 12))
90+
fig.suptitle(f"Scatter plot of speed_radius(km) ({m.sum()} matches)")
91+
92+
for i, (label, field, factor, stop) in enumerate(
93+
zip(
94+
("speed radius (km)", "outter radius (km)", "amplitude (cm)"),
95+
("radius_s", "radius_e", "amplitude"),
96+
(0.001, 0.001, 100),
97+
(80, 120, 25),
98+
)
99+
):
100+
ax = fig.add_subplot(
101+
2, 2, i + 1, xlabel="filtered grid", ylabel="raw grid", title=label
102+
)
103+
ax.plot(merge_f[field][i_] * factor, merge_r[field][j_] * factor, ".")
104+
ax.set_aspect("equal"), ax.grid()
105+
ax.plot((0, 1000), (0, 1000), "r")
106+
ax.set_xlim(0, stop), ax.set_ylim(0, stop)

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

Lines changed: 155 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,155 @@
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+
"\nEddy detection and filter\n=========================\n"
19+
]
20+
},
21+
{
22+
"cell_type": "code",
23+
"execution_count": null,
24+
"metadata": {
25+
"collapsed": false
26+
},
27+
"outputs": [],
28+
"source": [
29+
"from datetime import datetime\nfrom matplotlib import pyplot as plt\nfrom py_eddy_tracker.dataset.grid import RegularGridDataset\nfrom py_eddy_tracker import data\nfrom numpy import arange"
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+
"Load Input grid, ADT will be used to detect eddies\n\n"
48+
]
49+
},
50+
{
51+
"cell_type": "code",
52+
"execution_count": null,
53+
"metadata": {
54+
"collapsed": false
55+
},
56+
"outputs": [],
57+
"source": [
58+
"g = RegularGridDataset(\n data.get_path(\"dt_med_allsat_phy_l4_20160515_20190101.nc\"), \"longitude\", \"latitude\",\n)\ng.add_uv(\"adt\")\ng.copy(\"adt\", \"adt_high\")\nwavelength = 400\ng.bessel_high_filter(\"adt_high\", wavelength, order=3)\ndate = datetime(2016, 5, 15)"
59+
]
60+
},
61+
{
62+
"cell_type": "markdown",
63+
"metadata": {},
64+
"source": [
65+
"Run algorithm of detection\n\n"
66+
]
67+
},
68+
{
69+
"cell_type": "code",
70+
"execution_count": null,
71+
"metadata": {
72+
"collapsed": false
73+
},
74+
"outputs": [],
75+
"source": [
76+
"a_f, c_f = g.eddy_identification(\"adt_high\", \"u\", \"v\", date, 0.002)\nmerge_f = a_f.merge(c_f)\na_r, c_r = g.eddy_identification(\"adt\", \"u\", \"v\", date, 0.002)\nmerge_r = a_r.merge(c_r)"
77+
]
78+
},
79+
{
80+
"cell_type": "markdown",
81+
"metadata": {},
82+
"source": [
83+
"Display detection\n\n"
84+
]
85+
},
86+
{
87+
"cell_type": "code",
88+
"execution_count": null,
89+
"metadata": {
90+
"collapsed": false
91+
},
92+
"outputs": [],
93+
"source": [
94+
"ax = start_axes(\"Eddies detected over ADT\")\nm = g.display(ax, \"adt\", vmin=-0.15, vmax=0.15)\nmerge_f.display(ax, lw=0.5, label=\"Eddy from filtered grid\", ref=-10, color=\"k\")\nmerge_r.display(ax, lw=0.5, label=\"Eddy from raw grid\", ref=-10, color=\"r\")\nax.legend()\nupdate_axes(ax, m)"
95+
]
96+
},
97+
{
98+
"cell_type": "markdown",
99+
"metadata": {},
100+
"source": [
101+
"Parameters distribution\n-----------------------\n\n"
102+
]
103+
},
104+
{
105+
"cell_type": "code",
106+
"execution_count": null,
107+
"metadata": {
108+
"collapsed": false
109+
},
110+
"outputs": [],
111+
"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.grid(), ax_r.grid()\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()"
113+
]
114+
},
115+
{
116+
"cell_type": "markdown",
117+
"metadata": {},
118+
"source": [
119+
"Match detection and compare\n---------------------------\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, j, c = merge_f.match(merge_r)\nm = c > 0.1\ni_, j_ = i[m], j[m]\nfig = plt.figure(figsize=(12, 12))\nfig.suptitle(f\"Scatter plot of speed_radius(km) ({m.sum()} matches)\")\n\nfor i, (label, field, factor, stop) in enumerate(\n zip(\n (\"speed radius (km)\", \"outter radius (km)\", \"amplitude (cm)\"),\n (\"radius_s\", \"radius_e\", \"amplitude\"),\n (0.001, 0.001, 100),\n (80, 120, 25),\n )\n):\n ax = fig.add_subplot(\n 2, 2, i + 1, xlabel=\"filtered grid\", ylabel=\"raw grid\", title=label\n )\n ax.plot(merge_f[field][i_] * factor, merge_r[field][j_] * factor, \".\")\n ax.set_aspect(\"equal\"), ax.grid()\n ax.plot((0, 1000), (0, 1000), \"r\")\n ax.set_xlim(0, stop), ax.set_ylim(0, stop)"
131+
]
132+
}
133+
],
134+
"metadata": {
135+
"kernelspec": {
136+
"display_name": "Python 3",
137+
"language": "python",
138+
"name": "python3"
139+
},
140+
"language_info": {
141+
"codemirror_mode": {
142+
"name": "ipython",
143+
"version": 3
144+
},
145+
"file_extension": ".py",
146+
"mimetype": "text/x-python",
147+
"name": "python",
148+
"nbconvert_exporter": "python",
149+
"pygments_lexer": "ipython3",
150+
"version": "3.7.7"
151+
}
152+
},
153+
"nbformat": 4,
154+
"nbformat_minor": 0
155+
}

0 commit comments

Comments
 (0)