Skip to content

Commit cb87f3d

Browse files
AntSimiAntSimi
committed
Add circum polar example
1 parent 7565ca1 commit cb87f3d

File tree

2 files changed

+316
-0
lines changed

2 files changed

+316
-0
lines changed

‎examples/02_eddy_identification/pet_eddy_detection_ACC.py‎

Lines changed: 168 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,168 @@
1+
"""
2+
Eddy detection : Antartic circum polar
3+
======================================
4+
5+
Script will detect eddies on adt field, and compute u,v with method add_uv(which could use, only if equator is avoid)
6+
7+
Two ones with filtering adt and another without
8+
9+
"""
10+
from datetime import datetime
11+
12+
from matplotlib import pyplot as plt
13+
14+
from py_eddy_tracker import data
15+
from py_eddy_tracker.dataset.grid import RegularGridDataset
16+
17+
18+
def quad_axes(title):
19+
fig = plt.figure(figsize=(13, 8.5))
20+
fig.suptitle(title, weight="bold")
21+
axes = list()
22+
for position in (
23+
[0.05, 0.53, 0.44, 0.44],
24+
[0.53, 0.53, 0.44, 0.44],
25+
[0.05, 0.03, 0.44, 0.44],
26+
[0.53, 0.03, 0.44, 0.44],
27+
):
28+
ax = fig.add_axes(position)
29+
ax.set_xlim(5, 45), ax.set_ylim(-60, -37)
30+
ax.set_aspect("equal"), ax.grid(True)
31+
axes.append(ax)
32+
return axes
33+
34+
35+
# %%
36+
# Load Input grid, ADT is used to detect eddies
37+
margin = 30
38+
39+
kw_data = dict(
40+
filename=data.get_path("nrt_global_allsat_phy_l4_20190223_20190226.nc"),
41+
x_name="longitude",
42+
y_name="latitude",
43+
# Manual area subset
44+
indexs=dict(
45+
latitude=slice(100 - margin, 220 + margin),
46+
longitude=slice(0, 230 + margin),
47+
),
48+
)
49+
g_raw = RegularGridDataset(**kw_data)
50+
g_raw.add_uv("adt")
51+
g = RegularGridDataset(**kw_data)
52+
g.copy("adt", "adt_low")
53+
g.bessel_high_filter("adt", 700)
54+
g.bessel_low_filter("adt_low", 700)
55+
g.add_uv("adt")
56+
57+
# %%
58+
# Identification
59+
# ^^^^^^^^^^^^^^
60+
# Run the identification step with slices of 2 mm
61+
date = datetime(2016, 5, 15)
62+
kw_ident = dict(
63+
date=date, step=0.002, shape_error=70, sampling=30, uname="u", vname="v"
64+
)
65+
a, c = g.eddy_identification("adt", **kw_ident)
66+
a_, c_ = g_raw.eddy_identification("adt", **kw_ident)
67+
68+
69+
# %%
70+
# Figures
71+
# -------
72+
axs = quad_axes("General properties field")
73+
m = g_raw.display(axs[0], "adt", vmin=-1, vmax=1, cmap="RdBu_r")
74+
axs[0].set_title("ADT(m)")
75+
m = g.display(axs[1], "adt_low", vmin=-1, vmax=1, cmap="RdBu_r")
76+
axs[1].set_title("ADT (m) large scale with cut at 700 km")
77+
m = g.display(axs[2], "adt", vmin=-1, vmax=1, cmap="RdBu_r")
78+
axs[2].set_title("ADT (m) high scale with cut at 700 km")
79+
cb = plt.colorbar(
80+
m, cax=axs[0].figure.add_axes([0.03, 0.51, 0.94, 0.01]), orientation="horizontal"
81+
)
82+
cb.set_label("ADT(m)", labelpad=-2)
83+
84+
# %%
85+
axs = quad_axes("")
86+
axs[0].set_title("Without filter")
87+
axs[0].set_ylabel("Contours used in eddies")
88+
axs[1].set_title("With filter")
89+
axs[2].set_ylabel("Closed contours but not used")
90+
g_raw.contours.display(axs[0], lw=0.5, only_used=True)
91+
g.contours.display(axs[1], lw=0.5, only_used=True)
92+
g_raw.contours.display(axs[2], lw=0.5, only_unused=True)
93+
g.contours.display(axs[3], lw=0.5, only_unused=True)
94+
95+
# %%
96+
kw = dict(ref=-10, linewidth=0.75)
97+
kw_a = dict(color="r", label="Anticyclonic ({nb_obs} eddies)")
98+
kw_c = dict(color="b", label="Cyclonic ({nb_obs} eddies)")
99+
kw_filled = dict(vmin=0, vmax=100, cmap="Spectral_r", lut=20, intern=True, factor=100)
100+
axs = quad_axes("Comparison between two detection")
101+
# Match with intern/inner contour
102+
i_a, j_a, s_a = a_.match(a, intern=True, cmin=0.15)
103+
i_c, j_c, s_c = c_.match(c, intern=True, cmin=0.15)
104+
105+
a_.index(i_a).filled(axs[0], s_a, **kw_filled)
106+
a.index(j_a).filled(axs[1], s_a, **kw_filled)
107+
c_.index(i_c).filled(axs[0], s_c, **kw_filled)
108+
m = c.index(j_c).filled(axs[1], s_c, **kw_filled)
109+
110+
cb = plt.colorbar(
111+
m, cax=axs[0].figure.add_axes([0.03, 0.51, 0.94, 0.01]), orientation="horizontal"
112+
)
113+
cb.set_label("Similarity index", labelpad=-5)
114+
a_.display(axs[0], **kw, **kw_a), c_.display(axs[0], **kw, **kw_c)
115+
a.display(axs[1], **kw, **kw_a), c.display(axs[1], **kw, **kw_c)
116+
117+
axs[0].set_title("Without filter")
118+
axs[0].set_ylabel("Detection")
119+
axs[1].set_title("With filter")
120+
axs[2].set_ylabel("Contours' rejection criteria")
121+
122+
g_raw.contours.display(axs[2], lw=0.5, only_unused=True, display_criterion=True)
123+
g.contours.display(axs[3], lw=0.5, only_unused=True, display_criterion=True)
124+
125+
for ax in axs:
126+
ax.legend()
127+
# %%
128+
# Criteria for rejecting a contour :
129+
# 0. Accepted (green)
130+
# 1. Rejection for shape error (red)
131+
# 2. Masked value within contour (blue)
132+
# 3. Under or over the pixel limit bounds (black)
133+
# 4. Amplitude criterion (yellow)
134+
135+
# %%
136+
i_a, j_a = i_a[s_a >= 0.4], j_a[s_a >= 0.4]
137+
i_c, j_c = i_c[s_c >= 0.4], j_c[s_c >= 0.4]
138+
fig = plt.figure(figsize=(12, 12))
139+
fig.suptitle(f"Scatter plot (A : {i_a.shape[0]}, C : {i_c.shape[0]} matches)")
140+
141+
for i, (label, field, factor, stop) in enumerate(
142+
(
143+
("speed radius (km)", "radius_s", 0.001, 120),
144+
("outter radius (km)", "radius_e", 0.001, 120),
145+
("amplitude (cm)", "amplitude", 100, 25),
146+
("speed max (cm/s)", "speed_average", 100, 25),
147+
)
148+
):
149+
ax = fig.add_subplot(2, 2, i + 1, title=label)
150+
ax.set_xlabel("Without filter")
151+
ax.set_ylabel("With filter")
152+
153+
ax.plot(
154+
a_[field][i_a] * factor,
155+
a[field][j_a] * factor,
156+
"r.",
157+
label="Anticyclonic",
158+
)
159+
ax.plot(
160+
c_[field][i_c] * factor,
161+
c[field][j_c] * factor,
162+
"b.",
163+
label="Cyclonic",
164+
)
165+
ax.set_aspect("equal"), ax.grid()
166+
ax.plot((0, 1000), (0, 1000), "g")
167+
ax.set_xlim(0, stop), ax.set_ylim(0, stop)
168+
ax.legend()

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

Lines changed: 148 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,148 @@
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# Eddy detection : Antartic circum polar\n\nScript will detect eddies on adt field, and compute u,v with method add_uv(which could use, only if equator is avoid)\n\nTwo ones with filtering adt and another without\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\n\nfrom matplotlib import pyplot as plt\n\nfrom py_eddy_tracker import data\nfrom py_eddy_tracker.dataset.grid import RegularGridDataset\n\n\ndef quad_axes(title):\n fig = plt.figure(figsize=(13, 8.5))\n fig.suptitle(title, weight=\"bold\")\n axes = list()\n for position in (\n [0.05, 0.53, 0.44, 0.44],\n [0.53, 0.53, 0.44, 0.44],\n [0.05, 0.03, 0.44, 0.44],\n [0.53, 0.03, 0.44, 0.44],\n ):\n ax = fig.add_axes(position)\n ax.set_xlim(5, 45), ax.set_ylim(-60, -37)\n ax.set_aspect(\"equal\"), ax.grid(True)\n axes.append(ax)\n return axes"
30+
]
31+
},
32+
{
33+
"cell_type": "markdown",
34+
"metadata": {},
35+
"source": [
36+
"Load Input grid, ADT is used to detect eddies\n\n"
37+
]
38+
},
39+
{
40+
"cell_type": "code",
41+
"execution_count": null,
42+
"metadata": {
43+
"collapsed": false
44+
},
45+
"outputs": [],
46+
"source": [
47+
"margin = 30\n\nkw_data = dict(\n filename=data.get_path(\"nrt_global_allsat_phy_l4_20190223_20190226.nc\"),\n x_name=\"longitude\",\n y_name=\"latitude\",\n # Manual area subset\n indexs=dict(\n latitude=slice(100 - margin, 220 + margin),\n longitude=slice(0, 230 + margin),\n ),\n)\ng_raw = RegularGridDataset(**kw_data)\ng_raw.add_uv(\"adt\")\ng = RegularGridDataset(**kw_data)\ng.copy(\"adt\", \"adt_low\")\ng.bessel_high_filter(\"adt\", 700)\ng.bessel_low_filter(\"adt_low\", 700)\ng.add_uv(\"adt\")"
48+
]
49+
},
50+
{
51+
"cell_type": "markdown",
52+
"metadata": {},
53+
"source": [
54+
"## Identification\nRun the identification step with slices of 2 mm\n\n"
55+
]
56+
},
57+
{
58+
"cell_type": "code",
59+
"execution_count": null,
60+
"metadata": {
61+
"collapsed": false
62+
},
63+
"outputs": [],
64+
"source": [
65+
"date = datetime(2016, 5, 15)\nkw_ident = dict(\n date=date, step=0.002, shape_error=70, sampling=30, uname=\"u\", vname=\"v\"\n)\na, c = g.eddy_identification(\"adt\", **kw_ident)\na_, c_ = g_raw.eddy_identification(\"adt\", **kw_ident)"
66+
]
67+
},
68+
{
69+
"cell_type": "markdown",
70+
"metadata": {},
71+
"source": [
72+
"### Figures\n\n"
73+
]
74+
},
75+
{
76+
"cell_type": "code",
77+
"execution_count": null,
78+
"metadata": {
79+
"collapsed": false
80+
},
81+
"outputs": [],
82+
"source": [
83+
"axs = quad_axes(\"General properties field\")\nm = g_raw.display(axs[0], \"adt\", vmin=-1, vmax=1, cmap=\"RdBu_r\")\naxs[0].set_title(\"ADT(m)\")\nm = g.display(axs[1], \"adt_low\", vmin=-1, vmax=1, cmap=\"RdBu_r\")\naxs[1].set_title(\"ADT (m) large scale with cut at 700 km\")\nm = g.display(axs[2], \"adt\", vmin=-1, vmax=1, cmap=\"RdBu_r\")\naxs[2].set_title(\"ADT (m) high scale with cut at 700 km\")\ncb = plt.colorbar(\n m, cax=axs[0].figure.add_axes([0.03, 0.51, 0.94, 0.01]), orientation=\"horizontal\"\n)\ncb.set_label(\"ADT(m)\", labelpad=-2)"
84+
]
85+
},
86+
{
87+
"cell_type": "code",
88+
"execution_count": null,
89+
"metadata": {
90+
"collapsed": false
91+
},
92+
"outputs": [],
93+
"source": [
94+
"axs = quad_axes(\"\")\naxs[0].set_title(\"Without filter\")\naxs[0].set_ylabel(\"Contours used in eddies\")\naxs[1].set_title(\"With filter\")\naxs[2].set_ylabel(\"Closed contours but not used\")\ng_raw.contours.display(axs[0], lw=0.5, only_used=True)\ng.contours.display(axs[1], lw=0.5, only_used=True)\ng_raw.contours.display(axs[2], lw=0.5, only_unused=True)\ng.contours.display(axs[3], lw=0.5, only_unused=True)"
95+
]
96+
},
97+
{
98+
"cell_type": "code",
99+
"execution_count": null,
100+
"metadata": {
101+
"collapsed": false
102+
},
103+
"outputs": [],
104+
"source": [
105+
"kw = dict(ref=-10, linewidth=0.75)\nkw_a = dict(color=\"r\", label=\"Anticyclonic ({nb_obs} eddies)\")\nkw_c = dict(color=\"b\", label=\"Cyclonic ({nb_obs} eddies)\")\nkw_filled = dict(vmin=0, vmax=100, cmap=\"Spectral_r\", lut=20, intern=True, factor=100)\naxs = quad_axes(\"Comparison between two detection\")\n# Match with intern/inner contour\ni_a, j_a, s_a = a_.match(a, intern=True, cmin=0.15)\ni_c, j_c, s_c = c_.match(c, intern=True, cmin=0.15)\n\na_.index(i_a).filled(axs[0], s_a, **kw_filled)\na.index(j_a).filled(axs[1], s_a, **kw_filled)\nc_.index(i_c).filled(axs[0], s_c, **kw_filled)\nm = c.index(j_c).filled(axs[1], s_c, **kw_filled)\n\ncb = plt.colorbar(\n m, cax=axs[0].figure.add_axes([0.03, 0.51, 0.94, 0.01]), orientation=\"horizontal\"\n)\ncb.set_label(\"Similarity index\", labelpad=-5)\na_.display(axs[0], **kw, **kw_a), c_.display(axs[0], **kw, **kw_c)\na.display(axs[1], **kw, **kw_a), c.display(axs[1], **kw, **kw_c)\n\naxs[0].set_title(\"Without filter\")\naxs[0].set_ylabel(\"Detection\")\naxs[1].set_title(\"With filter\")\naxs[2].set_ylabel(\"Contours' rejection criteria\")\n\ng_raw.contours.display(axs[2], lw=0.5, only_unused=True, display_criterion=True)\ng.contours.display(axs[3], lw=0.5, only_unused=True, display_criterion=True)\n\nfor ax in axs:\n ax.legend()"
106+
]
107+
},
108+
{
109+
"cell_type": "markdown",
110+
"metadata": {},
111+
"source": [
112+
"Criteria for rejecting a contour :\n 0. Accepted (green)\n 1. Rejection for shape error (red)\n 2. Masked value within contour (blue)\n 3. Under or over the pixel limit bounds (black)\n 4. Amplitude criterion (yellow)\n\n"
113+
]
114+
},
115+
{
116+
"cell_type": "code",
117+
"execution_count": null,
118+
"metadata": {
119+
"collapsed": false
120+
},
121+
"outputs": [],
122+
"source": [
123+
"i_a, j_a = i_a[s_a >= 0.4], j_a[s_a >= 0.4]\ni_c, j_c = i_c[s_c >= 0.4], j_c[s_c >= 0.4]\nfig = plt.figure(figsize=(12, 12))\nfig.suptitle(f\"Scatter plot (A : {i_a.shape[0]}, C : {i_c.shape[0]} matches)\")\n\nfor i, (label, field, factor, stop) in enumerate(\n (\n (\"speed radius (km)\", \"radius_s\", 0.001, 120),\n (\"outter radius (km)\", \"radius_e\", 0.001, 120),\n (\"amplitude (cm)\", \"amplitude\", 100, 25),\n (\"speed max (cm/s)\", \"speed_average\", 100, 25),\n )\n):\n ax = fig.add_subplot(2, 2, i + 1, title=label)\n ax.set_xlabel(\"Without filter\")\n ax.set_ylabel(\"With filter\")\n\n ax.plot(\n a_[field][i_a] * factor,\n a[field][j_a] * factor,\n \"r.\",\n label=\"Anticyclonic\",\n )\n ax.plot(\n c_[field][i_c] * factor,\n c[field][j_c] * factor,\n \"b.\",\n label=\"Cyclonic\",\n )\n ax.set_aspect(\"equal\"), ax.grid()\n ax.plot((0, 1000), (0, 1000), \"g\")\n ax.set_xlim(0, stop), ax.set_ylim(0, stop)\n ax.legend()"
124+
]
125+
}
126+
],
127+
"metadata": {
128+
"kernelspec": {
129+
"display_name": "Python 3",
130+
"language": "python",
131+
"name": "python3"
132+
},
133+
"language_info": {
134+
"codemirror_mode": {
135+
"name": "ipython",
136+
"version": 3
137+
},
138+
"file_extension": ".py",
139+
"mimetype": "text/x-python",
140+
"name": "python",
141+
"nbconvert_exporter": "python",
142+
"pygments_lexer": "ipython3",
143+
"version": "3.7.7"
144+
}
145+
},
146+
"nbformat": 4,
147+
"nbformat_minor": 0
148+
}

0 commit comments

Comments
 (0)