Skip to content

Commit 8f9e9cf

Browse files
AntSimiAntSimi
committed
Add example: how work detection on gulf stream
1 parent 55860fb commit 8f9e9cf

File tree

5 files changed

+402
-5
lines changed

5 files changed

+402
-5
lines changed

‎examples/02_eddy_identification/pet_eddy_detection_gulf_stream.py‎

Lines changed: 131 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,131 @@
1+
"""
2+
Eddy detection Gulf stream
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+
Figures will show different step to detect eddies.
8+
9+
"""
10+
from datetime import datetime
11+
from matplotlib import pyplot as plt
12+
from py_eddy_tracker.dataset.grid import RegularGridDataset
13+
from py_eddy_tracker import data
14+
from py_eddy_tracker.eddy_feature import Contours
15+
16+
17+
# %%
18+
def start_axes(title):
19+
fig = plt.figure(figsize=(13, 8))
20+
ax = fig.add_axes([0.03, 0.03, 0.90, 0.94])
21+
ax.set_xlim(279, 304), ax.set_ylim(29, 44)
22+
ax.set_aspect("equal")
23+
ax.set_title(title)
24+
return ax
25+
26+
27+
def update_axes(ax, mappable=None):
28+
ax.grid()
29+
if mappable:
30+
plt.colorbar(mappable, cax=ax.figure.add_axes([0.95, 0.05, 0.01, 0.9]))
31+
32+
33+
# %%
34+
# Load Input grid, ADT will be used to detect eddies
35+
margin = 30
36+
g = RegularGridDataset(
37+
data.get_path("nrt_global_allsat_phy_l4_20190223_20190226.nc"),
38+
"longitude",
39+
"latitude",
40+
# Manual area subset
41+
indexs=dict(
42+
longitude=slice(1116 - margin, 1216 + margin),
43+
latitude=slice(476 - margin, 536 + margin),
44+
),
45+
)
46+
47+
ax = start_axes("ADT (m)")
48+
m = g.display(ax, "adt", vmin=-0.15, vmax=1)
49+
# Draw line on the gulf stream front
50+
great_current = Contours(g.x_c, g.y_c, g.grid("adt"), levels=(0.35,), keep_unclose=True)
51+
great_current.display(ax, color="k")
52+
update_axes(ax, m)
53+
54+
# %%
55+
# Get u/v
56+
# -------
57+
# U/V are deduced from ADT, this algortihm are not usable around equator (~+- 2°)
58+
g.add_uv("adt")
59+
60+
# %%
61+
# Pre-processings
62+
# ---------------
63+
# Apply high filter to remove long scale to highlight mesoscale
64+
g.bessel_high_filter("adt", 700)
65+
ax = start_axes("ADT (m) filtered (700km)")
66+
m = g.display(ax, "adt", vmin=-0.25, vmax=0.25)
67+
great_current.display(ax, color="k")
68+
update_axes(ax, m)
69+
70+
# %%
71+
# Identification
72+
# --------------
73+
# run identification with slice of 2 mm
74+
date = datetime(2016, 5, 15)
75+
a, c = g.eddy_identification("adt", "u", "v", date, 0.002)
76+
77+
# %%
78+
# All closed contour found in this input grid (Display only 1 contour every 5)
79+
ax = start_axes("ADT closed contour (only 1 / 5 levels)")
80+
g.contours.display(ax, step=5, lw=1)
81+
great_current.display(ax, color="k")
82+
update_axes(ax)
83+
84+
# %%
85+
# Contours include in eddies
86+
ax = start_axes("ADT contour used as eddies")
87+
g.contours.display(ax, only_used=True, lw=0.25)
88+
great_current.display(ax, color="k")
89+
update_axes(ax)
90+
91+
# %%
92+
# Contours reject from several origin (shape error to high, several extremum in contour, ...)
93+
ax = start_axes("ADT contour reject")
94+
g.contours.display(ax, only_unused=True, lw=0.25)
95+
great_current.display(ax, color="k")
96+
update_axes(ax)
97+
98+
# %%
99+
# Contours closed which contains several eddies
100+
ax = start_axes("ADT contour reject but which contain eddies")
101+
g.contours.label_contour_unused_which_contain_eddies(a)
102+
g.contours.label_contour_unused_which_contain_eddies(c)
103+
g.contours.display(
104+
ax, only_contain_eddies=True, color="k", lw=1, label="Could be interaction contour"
105+
)
106+
a.display(ax, color="r", linewidth=0.5, label="Anticyclonic", ref=-10)
107+
c.display(ax, color="b", linewidth=0.5, label="Cyclonic", ref=-10)
108+
ax.legend()
109+
update_axes(ax)
110+
111+
# %%
112+
# Output
113+
# ------
114+
# Display detected eddies, dashed lines represent effective contour
115+
# and solid lines represent contour of maximum of speed. See figure 1 of https://doi.org/10.1175/JTECH-D-14-00019.1
116+
117+
ax = start_axes("Eddies detected")
118+
a.display(ax, color="r", linewidth=0.5, label="Anticyclonic", ref=-10)
119+
c.display(ax, color="b", linewidth=0.5, label="Cyclonic", ref=-10)
120+
ax.legend()
121+
great_current.display(ax, color="k")
122+
update_axes(ax)
123+
124+
125+
# %%
126+
# Display speed radius of eddies detected
127+
ax = start_axes("Eddies speed radius (km)")
128+
a.filled(ax, "radius_e", vmin=10, vmax=150, cmap="magma_r", factor=0.001, lut=14)
129+
m = c.filled(ax, "radius_e", vmin=10, vmax=150, cmap="magma_r", factor=0.001, lut=14)
130+
great_current.display(ax, color="k")
131+
update_axes(ax, m)

‎notebooks/README.md‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,4 @@
11
python setup.py install
2+
# rm build/sphinx/ doc/python_module/ doc/gen_modules/ -rf
23
python setup.py build_sphinx
34
rsync -vrltp doc/python_module notebooks/. --include '*/' --include '*.ipynb' --exclude '*' --prune-empty-dirs

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

Lines changed: 245 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,245 @@
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 Gulf stream\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\nFigures will show different step to detect eddies.\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 py_eddy_tracker.eddy_feature import Contours"
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, 8))\n ax = fig.add_axes([0.03, 0.03, 0.90, 0.94])\n ax.set_xlim(279, 304), ax.set_ylim(29, 44)\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+
"margin = 30\ng = RegularGridDataset(\n data.get_path(\"nrt_global_allsat_phy_l4_20190223_20190226.nc\"),\n \"longitude\",\n \"latitude\",\n # Manual area subset\n indexs=dict(\n longitude=slice(1116 - margin, 1216 + margin),\n latitude=slice(476 - margin, 536 + margin),\n ),\n)\n\nax = start_axes(\"ADT (m)\")\nm = g.display(ax, \"adt\", vmin=-0.15, vmax=1)\n# Draw line on the gulf stream front\ngreat_current = Contours(g.x_c, g.y_c, g.grid(\"adt\"), levels=(0.35,), keep_unclose=True)\ngreat_current.display(ax, color=\"k\")\nupdate_axes(ax, m)"
59+
]
60+
},
61+
{
62+
"cell_type": "markdown",
63+
"metadata": {},
64+
"source": [
65+
"## Get u/v\nU/V are deduced from ADT, this algortihm are not usable around equator (~+- 2\u00b0)\n\n"
66+
]
67+
},
68+
{
69+
"cell_type": "code",
70+
"execution_count": null,
71+
"metadata": {
72+
"collapsed": false
73+
},
74+
"outputs": [],
75+
"source": [
76+
"g.add_uv(\"adt\")"
77+
]
78+
},
79+
{
80+
"cell_type": "markdown",
81+
"metadata": {},
82+
"source": [
83+
"## Pre-processings\nApply high filter to remove long scale to highlight mesoscale\n\n"
84+
]
85+
},
86+
{
87+
"cell_type": "code",
88+
"execution_count": null,
89+
"metadata": {
90+
"collapsed": false
91+
},
92+
"outputs": [],
93+
"source": [
94+
"g.bessel_high_filter(\"adt\", 700)\nax = start_axes(\"ADT (m) filtered (700km)\")\nm = g.display(ax, \"adt\", vmin=-0.25, vmax=0.25)\ngreat_current.display(ax, color=\"k\")\nupdate_axes(ax, m)"
95+
]
96+
},
97+
{
98+
"cell_type": "markdown",
99+
"metadata": {},
100+
"source": [
101+
"## Identification\nrun identification with slice of 2 mm\n\n"
102+
]
103+
},
104+
{
105+
"cell_type": "code",
106+
"execution_count": null,
107+
"metadata": {
108+
"collapsed": false
109+
},
110+
"outputs": [],
111+
"source": [
112+
"date = datetime(2016, 5, 15)\na, c = g.eddy_identification(\"adt\", \"u\", \"v\", date, 0.002)"
113+
]
114+
},
115+
{
116+
"cell_type": "markdown",
117+
"metadata": {},
118+
"source": [
119+
"All closed contour found in this input grid (Display only 1 contour every 5)\n\n"
120+
]
121+
},
122+
{
123+
"cell_type": "code",
124+
"execution_count": null,
125+
"metadata": {
126+
"collapsed": false
127+
},
128+
"outputs": [],
129+
"source": [
130+
"ax = start_axes(\"ADT closed contour (only 1 / 5 levels)\")\ng.contours.display(ax, step=5, lw=1)\ngreat_current.display(ax, color=\"k\")\nupdate_axes(ax)"
131+
]
132+
},
133+
{
134+
"cell_type": "markdown",
135+
"metadata": {},
136+
"source": [
137+
"Contours include in eddies\n\n"
138+
]
139+
},
140+
{
141+
"cell_type": "code",
142+
"execution_count": null,
143+
"metadata": {
144+
"collapsed": false
145+
},
146+
"outputs": [],
147+
"source": [
148+
"ax = start_axes(\"ADT contour used as eddies\")\ng.contours.display(ax, only_used=True, lw=0.25)\ngreat_current.display(ax, color=\"k\")\nupdate_axes(ax)"
149+
]
150+
},
151+
{
152+
"cell_type": "markdown",
153+
"metadata": {},
154+
"source": [
155+
"Contours reject from several origin (shape error to high, several extremum in contour, ...)\n\n"
156+
]
157+
},
158+
{
159+
"cell_type": "code",
160+
"execution_count": null,
161+
"metadata": {
162+
"collapsed": false
163+
},
164+
"outputs": [],
165+
"source": [
166+
"ax = start_axes(\"ADT contour reject\")\ng.contours.display(ax, only_unused=True, lw=0.25)\ngreat_current.display(ax, color=\"k\")\nupdate_axes(ax)"
167+
]
168+
},
169+
{
170+
"cell_type": "markdown",
171+
"metadata": {},
172+
"source": [
173+
"Contours closed which contains several eddies\n\n"
174+
]
175+
},
176+
{
177+
"cell_type": "code",
178+
"execution_count": null,
179+
"metadata": {
180+
"collapsed": false
181+
},
182+
"outputs": [],
183+
"source": [
184+
"ax = start_axes(\"ADT contour reject but which contain eddies\")\ng.contours.label_contour_unused_which_contain_eddies(a)\ng.contours.label_contour_unused_which_contain_eddies(c)\ng.contours.display(\n ax, only_contain_eddies=True, color=\"k\", lw=1, label=\"Could be interaction contour\"\n)\na.display(ax, color=\"r\", linewidth=0.5, label=\"Anticyclonic\", ref=-10)\nc.display(ax, color=\"b\", linewidth=0.5, label=\"Cyclonic\", ref=-10)\nax.legend()\nupdate_axes(ax)"
185+
]
186+
},
187+
{
188+
"cell_type": "markdown",
189+
"metadata": {},
190+
"source": [
191+
"## Output\nDisplay detected eddies, dashed lines represent effective contour\nand solid lines represent contour of maximum of speed. See figure 1 of https://doi.org/10.1175/JTECH-D-14-00019.1\n\n"
192+
]
193+
},
194+
{
195+
"cell_type": "code",
196+
"execution_count": null,
197+
"metadata": {
198+
"collapsed": false
199+
},
200+
"outputs": [],
201+
"source": [
202+
"ax = start_axes(\"Eddies detected\")\na.display(ax, color=\"r\", linewidth=0.5, label=\"Anticyclonic\", ref=-10)\nc.display(ax, color=\"b\", linewidth=0.5, label=\"Cyclonic\", ref=-10)\nax.legend()\ngreat_current.display(ax, color=\"k\")\nupdate_axes(ax)"
203+
]
204+
},
205+
{
206+
"cell_type": "markdown",
207+
"metadata": {},
208+
"source": [
209+
"Display speed radius of eddies detected\n\n"
210+
]
211+
},
212+
{
213+
"cell_type": "code",
214+
"execution_count": null,
215+
"metadata": {
216+
"collapsed": false
217+
},
218+
"outputs": [],
219+
"source": [
220+
"ax = start_axes(\"Eddies speed radius (km)\")\na.filled(ax, \"radius_e\", vmin=10, vmax=150, cmap=\"magma_r\", factor=0.001, lut=14)\nm = c.filled(ax, \"radius_e\", vmin=10, vmax=150, cmap=\"magma_r\", factor=0.001, lut=14)\ngreat_current.display(ax, color=\"k\")\nupdate_axes(ax, m)"
221+
]
222+
}
223+
],
224+
"metadata": {
225+
"kernelspec": {
226+
"display_name": "Python 3",
227+
"language": "python",
228+
"name": "python3"
229+
},
230+
"language_info": {
231+
"codemirror_mode": {
232+
"name": "ipython",
233+
"version": 3
234+
},
235+
"file_extension": ".py",
236+
"mimetype": "text/x-python",
237+
"name": "python",
238+
"nbconvert_exporter": "python",
239+
"pygments_lexer": "ipython3",
240+
"version": "3.7.7"
241+
}
242+
},
243+
"nbformat": 4,
244+
"nbformat_minor": 0
245+
}

0 commit comments

Comments
 (0)