Skip to content

Commit 310972a

Browse files
AntSimiAntSimi
committed
Add okubo weiss example
1 parent 2b3200c commit 310972a

File tree

5 files changed

+390
-6
lines changed

5 files changed

+390
-6
lines changed

‎doc/conf.py‎

Lines changed: 0 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -59,8 +59,6 @@
5959
"doc_module": ("py_eddy_tracker",),
6060
"reference_url": {
6161
"py_eddy_tracker": None,
62-
"matplotlib": "https://matplotlib.org/",
63-
"numpy": "https://docs.scipy.org/doc/numpy/",
6462
},
6563
"line_numbers": False,
6664
"filename_pattern": "/pet",

‎examples/06_grid_manipulation/pet_okubo_weiss.py‎

Lines changed: 155 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,155 @@
1+
r"""
2+
Get Okubo Weis
3+
=====================
4+
5+
.. math:: OW = S_n^2 + S_s^2 + \omega^2
6+
7+
with normal strain (:math:`S_n`), shear strain (:math:`S_s`) and vorticity (:math:`\omega`)
8+
9+
.. math::
10+
S_n = \frac{\partial u}{\partial x} - \frac{\partial v}{\partial y},
11+
S_s = \frac{\partial v}{\partial x} + \frac{\partial u}{\partial y},
12+
\omega = \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y}
13+
14+
"""
15+
from matplotlib import pyplot as plt
16+
from py_eddy_tracker.dataset.grid import RegularGridDataset
17+
from py_eddy_tracker.observations.observation import EddiesObservations
18+
from py_eddy_tracker import data
19+
from numpy import arange, ma, where
20+
21+
22+
# %%
23+
def start_axes(title, zoom=False):
24+
fig = plt.figure(figsize=(12, 6))
25+
axes = fig.add_axes([0.03, 0.03, 0.90, 0.94])
26+
axes.set_xlim(0, 360), axes.set_ylim(-80, 80)
27+
if zoom:
28+
axes.set_xlim(270, 340), axes.set_ylim(20, 50)
29+
axes.set_aspect("equal")
30+
axes.set_title(title)
31+
return axes
32+
33+
34+
def update_axes(axes, mappable=None):
35+
axes.grid()
36+
if mappable:
37+
plt.colorbar(mappable, cax=axes.figure.add_axes([0.94, 0.05, 0.01, 0.9]))
38+
39+
40+
# %%
41+
# Load detection files
42+
a = EddiesObservations.load_file(data.get_path("Anticyclonic_20190223.nc"))
43+
c = EddiesObservations.load_file(data.get_path("Cyclonic_20190223.nc"))
44+
45+
# %%
46+
# Load Input grid, ADT will be used to detect eddies
47+
g = RegularGridDataset(
48+
data.get_path("nrt_global_allsat_phy_l4_20190223_20190226.nc"),
49+
"longitude",
50+
"latitude",
51+
)
52+
53+
ax = start_axes("ADT (cm)")
54+
m = g.display(ax, "adt", vmin=-120, vmax=120, factor=100)
55+
update_axes(ax, m)
56+
57+
# %%
58+
# Get parameter for ow
59+
u_x = g.compute_stencil(g.grid("ugos"))
60+
u_y = g.compute_stencil(g.grid("ugos"), vertical=True)
61+
v_x = g.compute_stencil(g.grid("vgos"))
62+
v_y = g.compute_stencil(g.grid("vgos"), vertical=True)
63+
ow = g.vars["ow"] = (u_x - v_y) ** 2 + (v_x + u_y) ** 2 - (v_x - u_y) ** 2
64+
65+
66+
ax = start_axes("Okubo weis")
67+
m = g.display(ax, "ow", vmin=-1e-10, vmax=1e-10, cmap="bwr")
68+
update_axes(ax, m)
69+
70+
# %%
71+
# Gulf stream zoom
72+
ax = start_axes("Okubo weis, Gulf stream", zoom=True)
73+
m = g.display(ax, "ow", vmin=-1e-10, vmax=1e-10, cmap="bwr")
74+
kw_ed = dict(intern_only=True, color="k", lw=1)
75+
a.display(ax, **kw_ed), c.display(ax, **kw_ed)
76+
update_axes(ax, m)
77+
78+
# %%
79+
# only negative OW
80+
ax = start_axes("Okubo weis, Gulf stream", zoom=True)
81+
threshold = ow.std() * -0.2
82+
ow = ma.array(ow, mask=ow > threshold)
83+
m = g.display(ax, ow, vmin=-1e-10, vmax=1e-10, cmap="bwr")
84+
a.display(ax, **kw_ed), c.display(ax, **kw_ed)
85+
update_axes(ax, m)
86+
87+
88+
# %%
89+
# Get okubo-weiss mean/min/center in eddies
90+
plt.figure(figsize=(8, 6))
91+
ax = plt.subplot(111)
92+
ax.set_xlabel("Okubo-Weiss parameter")
93+
kw_hist = dict(bins=arange(-20e-10, 20e-10, 50e-12), histtype="step")
94+
for method in ("mean", "center", "min"):
95+
kw_interp = dict(grid_object=g, varname="ow", method=method, intern=True)
96+
_, _, m = ax.hist(
97+
a.interp_grid(**kw_interp), label=f"Anticyclonic - OW {method}", **kw_hist
98+
)
99+
ax.hist(
100+
c.interp_grid(**kw_interp),
101+
label=f"Cyclonic - OW {method}",
102+
color=m[0].get_edgecolor(),
103+
ls="--",
104+
**kw_hist,
105+
)
106+
ax.axvline(threshold, color="r")
107+
ax.set_yscale("log")
108+
ax.grid()
109+
ax.set_ylim(1, 1e4)
110+
ax.set_xlim(-15e-10, 15e-10)
111+
ax.legend()
112+
113+
# %%
114+
# Catch eddies with bad OW
115+
ax = start_axes("Eddies with a min OW in speed contour over threshold")
116+
ow_min = a.interp_grid(**kw_interp)
117+
a_bad_ow = a.index(where(ow_min > threshold)[0])
118+
a_bad_ow.display(ax, color="r", label="Anticyclonic")
119+
ow_min = c.interp_grid(**kw_interp)
120+
c_bad_ow = c.index(where(ow_min > threshold)[0])
121+
c_bad_ow.display(ax, color="b", label="Cyclonic")
122+
ax.legend()
123+
124+
# %%
125+
# Display Radius and amplitude of eddies
126+
fig = plt.figure(figsize=(12, 5))
127+
fig.suptitle(
128+
"Parameter distribution (solid line) and cumulative distribution (dashed line)"
129+
)
130+
ax_amp, ax_rad = fig.add_subplot(121), fig.add_subplot(122)
131+
ax_amp_c, ax_rad_c = ax_amp.twinx(), ax_rad.twinx()
132+
ax_amp_c.set_ylim(0, 1), ax_rad_c.set_ylim(0, 1)
133+
kw_a = dict(xname="amplitude", bins=arange(0, 2, 0.002).astype("f4"))
134+
kw_r = dict(xname="radius_s", bins=arange(0, 500e6, 2e3).astype("f4"))
135+
for d, label, color in (
136+
(a, "Anticyclonic all", "r"),
137+
(a_bad_ow, "Anticyclonic bad OW", "orange"),
138+
(c, "Cyclonic all", "blue"),
139+
(c_bad_ow, "Cyclonic bad OW", "lightblue"),
140+
):
141+
x, y = d.bins_stat(**kw_a)
142+
ax_amp.plot(x * 100, y, label=label, color=color)
143+
ax_amp_c.plot(
144+
x * 100, y.cumsum() / y.sum(), label=label, color=color, ls="-.", lw=0.5
145+
)
146+
x, y = d.bins_stat(**kw_r)
147+
ax_rad.plot(x * 1e-3, y, label=label, color=color)
148+
ax_rad_c.plot(
149+
x * 1e-3, y.cumsum() / y.sum(), label=label, color=color, ls="-.", lw=0.5
150+
)
151+
152+
ax_amp.set_xlim(0, 12.5), ax_amp.grid(), ax_amp.set_ylim(0), ax_amp.legend()
153+
ax_rad.set_xlim(0, 120), ax_rad.grid(), ax_rad.set_ylim(0)
154+
ax_amp.set_xlabel("Amplitude (cm)"), ax_amp.set_ylabel("Nb eddies")
155+
ax_rad.set_xlabel("Speed radius (km)")

‎notebooks/python_module/06_grid_manipulation/pet_okubo_weiss.ipynb‎

Lines changed: 209 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,209 @@
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# Get Okubo Weis\n\n\\begin{align}OW = S_n^2 + S_s^2 + \\omega^2\\end{align}\n\nwith normal strain ($S_n$), shear strain ($S_s$) and vorticity ($\\omega$)\n\n\\begin{align}S_n = \\frac{\\partial u}{\\partial x} - \\frac{\\partial v}{\\partial y},\n S_s = \\frac{\\partial v}{\\partial x} + \\frac{\\partial u}{\\partial y},\n \\omega = \\frac{\\partial v}{\\partial x} - \\frac{\\partial u}{\\partial y}\\end{align}\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.dataset.grid import RegularGridDataset\nfrom py_eddy_tracker.observations.observation import EddiesObservations\nfrom py_eddy_tracker import data\nfrom numpy import arange, ma, where"
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, zoom=False):\n fig = plt.figure(figsize=(12, 6))\n axes = fig.add_axes([0.03, 0.03, 0.90, 0.94])\n axes.set_xlim(0, 360), axes.set_ylim(-80, 80)\n if zoom:\n axes.set_xlim(270, 340), axes.set_ylim(20, 50)\n axes.set_aspect(\"equal\")\n axes.set_title(title)\n return axes\n\n\ndef update_axes(axes, mappable=None):\n axes.grid()\n if mappable:\n plt.colorbar(mappable, cax=axes.figure.add_axes([0.94, 0.05, 0.01, 0.9]))"
41+
]
42+
},
43+
{
44+
"cell_type": "markdown",
45+
"metadata": {},
46+
"source": [
47+
"Load detection files\n\n"
48+
]
49+
},
50+
{
51+
"cell_type": "code",
52+
"execution_count": null,
53+
"metadata": {
54+
"collapsed": false
55+
},
56+
"outputs": [],
57+
"source": [
58+
"a = EddiesObservations.load_file(data.get_path(\"Anticyclonic_20190223.nc\"))\nc = EddiesObservations.load_file(data.get_path(\"Cyclonic_20190223.nc\"))"
59+
]
60+
},
61+
{
62+
"cell_type": "markdown",
63+
"metadata": {},
64+
"source": [
65+
"Load Input grid, ADT will be used to detect eddies\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 = RegularGridDataset(\n data.get_path(\"nrt_global_allsat_phy_l4_20190223_20190226.nc\"),\n \"longitude\",\n \"latitude\",\n)\n\nax = start_axes(\"ADT (cm)\")\nm = g.display(ax, \"adt\", vmin=-120, vmax=120, factor=100)\nupdate_axes(ax, m)"
77+
]
78+
},
79+
{
80+
"cell_type": "markdown",
81+
"metadata": {},
82+
"source": [
83+
"Get parameter for ow\n\n"
84+
]
85+
},
86+
{
87+
"cell_type": "code",
88+
"execution_count": null,
89+
"metadata": {
90+
"collapsed": false
91+
},
92+
"outputs": [],
93+
"source": [
94+
"u_x = g.compute_stencil(g.grid(\"ugos\"))\nu_y = g.compute_stencil(g.grid(\"ugos\"), vertical=True)\nv_x = g.compute_stencil(g.grid(\"vgos\"))\nv_y = g.compute_stencil(g.grid(\"vgos\"), vertical=True)\now = g.vars[\"ow\"] = (u_x - v_y) ** 2 + (v_x + u_y) ** 2 - (v_x - u_y) ** 2\n\n\nax = start_axes(\"Okubo weis\")\nm = g.display(ax, \"ow\", vmin=-1e-10, vmax=1e-10, cmap=\"bwr\")\nupdate_axes(ax, m)"
95+
]
96+
},
97+
{
98+
"cell_type": "markdown",
99+
"metadata": {},
100+
"source": [
101+
"Gulf stream zoom\n\n"
102+
]
103+
},
104+
{
105+
"cell_type": "code",
106+
"execution_count": null,
107+
"metadata": {
108+
"collapsed": false
109+
},
110+
"outputs": [],
111+
"source": [
112+
"ax = start_axes(\"Okubo weis, Gulf stream\", zoom=True)\nm = g.display(ax, \"ow\", vmin=-1e-10, vmax=1e-10, cmap=\"bwr\")\nkw_ed = dict(intern_only=True, color=\"k\", lw=1)\na.display(ax, **kw_ed), c.display(ax, **kw_ed)\nupdate_axes(ax, m)"
113+
]
114+
},
115+
{
116+
"cell_type": "markdown",
117+
"metadata": {},
118+
"source": [
119+
"only negative OW\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(\"Okubo weis, Gulf stream\", zoom=True)\nthreshold = ow.std() * -0.2\now = ma.array(ow, mask=ow > threshold)\nm = g.display(ax, ow, vmin=-1e-10, vmax=1e-10, cmap=\"bwr\")\na.display(ax, **kw_ed), c.display(ax, **kw_ed)\nupdate_axes(ax, m)"
131+
]
132+
},
133+
{
134+
"cell_type": "markdown",
135+
"metadata": {},
136+
"source": [
137+
"Get okubo-weiss mean/min/center 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+
"plt.figure(figsize=(8, 6))\nax = plt.subplot(111)\nax.set_xlabel(\"Okubo-Weiss parameter\")\nkw_hist = dict(bins=arange(-20e-10, 20e-10, 50e-12), histtype=\"step\")\nfor method in (\"mean\", \"center\", \"min\"):\n kw_interp = dict(grid_object=g, varname=\"ow\", method=method, intern=True)\n _, _, m = ax.hist(\n a.interp_grid(**kw_interp), label=f\"Anticyclonic - OW {method}\", **kw_hist\n )\n ax.hist(\n c.interp_grid(**kw_interp),\n label=f\"Cyclonic - OW {method}\",\n color=m[0].get_edgecolor(),\n ls=\"--\",\n **kw_hist,\n )\nax.axvline(threshold, color=\"r\")\nax.set_yscale(\"log\")\nax.grid()\nax.set_ylim(1, 1e4)\nax.set_xlim(-15e-10, 15e-10)\nax.legend()"
149+
]
150+
},
151+
{
152+
"cell_type": "markdown",
153+
"metadata": {},
154+
"source": [
155+
"Catch eddies with bad OW\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(\"Eddies with a min OW in speed contour over threshold\")\now_min = a.interp_grid(**kw_interp)\na_bad_ow = a.index(where(ow_min > threshold)[0])\na_bad_ow.display(ax, color=\"r\", label=\"Anticyclonic\")\now_min = c.interp_grid(**kw_interp)\nc_bad_ow = c.index(where(ow_min > threshold)[0])\nc_bad_ow.display(ax, color=\"b\", label=\"Cyclonic\")\nax.legend()"
167+
]
168+
},
169+
{
170+
"cell_type": "markdown",
171+
"metadata": {},
172+
"source": [
173+
"Display Radius and amplitude of 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+
"fig = plt.figure(figsize=(12, 5))\nfig.suptitle(\n \"Parameter distribution (solid line) and cumulative distribution (dashed line)\"\n)\nax_amp, ax_rad = fig.add_subplot(121), fig.add_subplot(122)\nax_amp_c, ax_rad_c = ax_amp.twinx(), ax_rad.twinx()\nax_amp_c.set_ylim(0, 1), ax_rad_c.set_ylim(0, 1)\nkw_a = dict(xname=\"amplitude\", bins=arange(0, 2, 0.002).astype(\"f4\"))\nkw_r = dict(xname=\"radius_s\", bins=arange(0, 500e6, 2e3).astype(\"f4\"))\nfor d, label, color in (\n (a, \"Anticyclonic all\", \"r\"),\n (a_bad_ow, \"Anticyclonic bad OW\", \"orange\"),\n (c, \"Cyclonic all\", \"blue\"),\n (c_bad_ow, \"Cyclonic bad OW\", \"lightblue\"),\n):\n x, y = d.bins_stat(**kw_a)\n ax_amp.plot(x * 100, y, label=label, color=color)\n ax_amp_c.plot(\n x * 100, y.cumsum() / y.sum(), label=label, color=color, ls=\"-.\", lw=0.5\n )\n x, y = d.bins_stat(**kw_r)\n ax_rad.plot(x * 1e-3, y, label=label, color=color)\n ax_rad_c.plot(\n x * 1e-3, y.cumsum() / y.sum(), label=label, color=color, ls=\"-.\", lw=0.5\n )\n\nax_amp.set_xlim(0, 12.5), ax_amp.grid(), ax_amp.set_ylim(0), ax_amp.legend()\nax_rad.set_xlim(0, 120), ax_rad.grid(), ax_rad.set_ylim(0)\nax_amp.set_xlabel(\"Amplitude (cm)\"), ax_amp.set_ylabel(\"Nb eddies\")\nax_rad.set_xlabel(\"Speed radius (km)\")"
185+
]
186+
}
187+
],
188+
"metadata": {
189+
"kernelspec": {
190+
"display_name": "Python 3",
191+
"language": "python",
192+
"name": "python3"
193+
},
194+
"language_info": {
195+
"codemirror_mode": {
196+
"name": "ipython",
197+
"version": 3
198+
},
199+
"file_extension": ".py",
200+
"mimetype": "text/x-python",
201+
"name": "python",
202+
"nbconvert_exporter": "python",
203+
"pygments_lexer": "ipython3",
204+
"version": "3.7.7"
205+
}
206+
},
207+
"nbformat": 4,
208+
"nbformat_minor": 0
209+
}

0 commit comments

Comments
 (0)