Skip to content

Commit 2780f46

Browse files
AntSimiAntSimi
committed
Add notebook about filtering
1 parent 2aa6191 commit 2780f46

File tree

4 files changed

+354
-2
lines changed

4 files changed

+354
-2
lines changed

‎examples/06_grid_manipulation/pet_filter.py‎

Lines changed: 144 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,144 @@
1+
"""
2+
Grid filtering in PET
3+
=====================
4+
5+
How filter work in py eddy tracker. This implementation maybe doesn't respect state art, but ...
6+
7+
We code a specific filter in order to filter grid with same wavelength at each pixel.
8+
"""
9+
10+
from py_eddy_tracker.dataset.grid import RegularGridDataset
11+
from py_eddy_tracker import data
12+
from matplotlib import pyplot as plt
13+
14+
15+
def start_axes(title):
16+
fig = plt.figure(figsize=(13, 5))
17+
ax = fig.add_axes([0.03, 0.03, 0.90, 0.94])
18+
ax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)
19+
ax.set_aspect("equal")
20+
ax.set_title(title)
21+
return ax
22+
23+
24+
def update_axes(ax, mappable=None):
25+
ax.grid()
26+
if mappable:
27+
plt.colorbar(m, cax=ax.figure.add_axes([0.95, 0.05, 0.01, 0.9]))
28+
29+
30+
# %%
31+
# All information will be for regular grid
32+
g = RegularGridDataset(
33+
data.get_path("dt_med_allsat_phy_l4_20160515_20190101.nc"), "longitude", "latitude",
34+
)
35+
# %%
36+
# Kernel
37+
# ------
38+
# Shape of kernel will increase in x when latitude increase
39+
fig = plt.figure(figsize=(12, 8))
40+
for i, latitude in enumerate((15, 35, 55, 75)):
41+
k = g.kernel_bessel(latitude, 500, order=3).T
42+
ax0 = plt.subplot(
43+
2,
44+
2,
45+
i + 1,
46+
title=f"Kernel at {latitude}° of latitude\nfor 1/8° grid, shape : {k.shape}",
47+
aspect="equal",
48+
)
49+
m = ax0.pcolormesh(k, vmin=-0.5, vmax=2, cmap="viridis_r")
50+
plt.colorbar(m, cax=fig.add_axes((0.92, 0.05, 0.01, 0.9)))
51+
52+
# %%
53+
# Kernel applying
54+
# ---------------
55+
# Original grid
56+
ax = start_axes("ADT")
57+
m = g.display(ax, "adt", vmin=-0.15, vmax=0.15)
58+
update_axes(ax, m)
59+
60+
# %%
61+
# We will select wavelength of 300 km
62+
#
63+
# Low frequency
64+
ax = start_axes("ADT low frequency")
65+
g.copy("adt", "adt_low")
66+
g.bessel_low_filter("adt_low", 300, order=3)
67+
m = g.display(ax, "adt_low", vmin=-0.15, vmax=0.15)
68+
update_axes(ax, m)
69+
70+
# %%
71+
# High frequency
72+
ax = start_axes("ADT high frequency")
73+
g.copy("adt", "adt_high")
74+
g.bessel_high_filter("adt_high", 300, order=3)
75+
m = g.display(ax, "adt_high", vmin=-0.15, vmax=0.15)
76+
update_axes(ax, m)
77+
78+
# %%
79+
# Clues
80+
# -----
81+
# wavelength : 80km
82+
83+
g.copy("adt", "adt_high_80")
84+
g.bessel_high_filter("adt_high_80", 80, order=3)
85+
g.copy("adt", "adt_low_80")
86+
g.bessel_low_filter("adt_low_80", 80, order=3)
87+
88+
area = dict(llcrnrlon=11.75, urcrnrlon=21, llcrnrlat=33, urcrnrlat=36.75)
89+
90+
# %%
91+
# Spectrum
92+
fig = plt.figure(figsize=(10, 6))
93+
ax = fig.add_subplot(111)
94+
ax.set_title("Spectrum")
95+
ax.set_xlabel("km")
96+
97+
lon_spec, lat_spec = g.spectrum_lonlat("adt", area=area)
98+
mappable = ax.loglog(*lat_spec, label="lat raw")[0]
99+
ax.loglog(*lon_spec, label="lon raw", color=mappable.get_color(), linestyle="--")
100+
101+
lon_spec, lat_spec = g.spectrum_lonlat("adt_high_80", area=area)
102+
mappable = ax.loglog(*lat_spec, label="lat high")[0]
103+
ax.loglog(*lon_spec, label="lon high", color=mappable.get_color(), linestyle="--")
104+
105+
lon_spec, lat_spec = g.spectrum_lonlat("adt_low_80", area=area)
106+
mappable = ax.loglog(*lat_spec, label="lat low")[0]
107+
ax.loglog(*lon_spec, label="lon low", color=mappable.get_color(), linestyle="--")
108+
109+
ax.set_xlim(10, 1000)
110+
ax.set_ylim(1e-6, 1)
111+
ax.set_xscale("log")
112+
ax.legend()
113+
ax.grid()
114+
115+
# %%
116+
# Spectrum ratio
117+
fig = plt.figure(figsize=(10, 6))
118+
ax = fig.add_subplot(111)
119+
ax.set_title("Spectrum ratio")
120+
ax.set_xlabel("km")
121+
122+
lon_spec, lat_spec = g.spectrum_lonlat(
123+
"adt_high_80", area=area, ref=g, ref_grid_name="adt"
124+
)
125+
mappable = ax.plot(*lat_spec, label="lat high")[0]
126+
ax.plot(*lon_spec, label="lon high", color=mappable.get_color(), linestyle="--")
127+
128+
lon_spec, lat_spec = g.spectrum_lonlat(
129+
"adt_low_80", area=area, ref=g, ref_grid_name="adt"
130+
)
131+
mappable = ax.plot(*lat_spec, label="lat low")[0]
132+
ax.plot(*lon_spec, label="lon low", color=mappable.get_color(), linestyle="--")
133+
134+
ax.set_xlim(10, 1000)
135+
ax.set_ylim(0, 1)
136+
ax.set_xscale("log")
137+
ax.legend()
138+
ax.grid()
139+
140+
# %%
141+
# Old filter
142+
# ----------
143+
# To do ...
144+

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

Lines changed: 205 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,205 @@
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+
"\nGrid filtering in PET\n=====================\n\nHow filter work in py eddy tracker. This implementation maybe doesn't respect state art, but ...\n\nWe code a specific filter in order to filter grid with same wavelength at each pixel.\n"
19+
]
20+
},
21+
{
22+
"cell_type": "code",
23+
"execution_count": null,
24+
"metadata": {
25+
"collapsed": false
26+
},
27+
"outputs": [],
28+
"source": [
29+
"from py_eddy_tracker.dataset.grid import RegularGridDataset\nfrom py_eddy_tracker import data\nfrom matplotlib import pyplot as plt\n\n\ndef 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(m, cax=ax.figure.add_axes([0.95, 0.05, 0.01, 0.9]))"
30+
]
31+
},
32+
{
33+
"cell_type": "markdown",
34+
"metadata": {},
35+
"source": [
36+
"All information will be for regular grid\n\n"
37+
]
38+
},
39+
{
40+
"cell_type": "code",
41+
"execution_count": null,
42+
"metadata": {
43+
"collapsed": false
44+
},
45+
"outputs": [],
46+
"source": [
47+
"g = RegularGridDataset(\n data.get_path(\"dt_med_allsat_phy_l4_20160515_20190101.nc\"), \"longitude\", \"latitude\",\n)"
48+
]
49+
},
50+
{
51+
"cell_type": "markdown",
52+
"metadata": {},
53+
"source": [
54+
"Kernel\n------\nShape of kernel will increase in x when latitude increase\n\n"
55+
]
56+
},
57+
{
58+
"cell_type": "code",
59+
"execution_count": null,
60+
"metadata": {
61+
"collapsed": false
62+
},
63+
"outputs": [],
64+
"source": [
65+
"fig = plt.figure(figsize=(12, 8))\nfor i, latitude in enumerate((15, 35, 55, 75)):\n k = g.kernel_bessel(latitude, 500, order=3).T\n ax0 = plt.subplot(\n 2,\n 2,\n i + 1,\n title=f\"Kernel at {latitude}\u00b0 of latitude\\nfor 1/8\u00b0 grid, shape : {k.shape}\",\n aspect=\"equal\",\n )\n m = ax0.pcolormesh(k, vmin=-0.5, vmax=2, cmap=\"viridis_r\")\nplt.colorbar(m, cax=fig.add_axes((0.92, 0.05, 0.01, 0.9)))"
66+
]
67+
},
68+
{
69+
"cell_type": "markdown",
70+
"metadata": {},
71+
"source": [
72+
"Kernel applying\n---------------\nOriginal grid\n\n"
73+
]
74+
},
75+
{
76+
"cell_type": "code",
77+
"execution_count": null,
78+
"metadata": {
79+
"collapsed": false
80+
},
81+
"outputs": [],
82+
"source": [
83+
"ax = start_axes(\"ADT\")\nm = g.display(ax, \"adt\", vmin=-0.15, vmax=0.15)\nupdate_axes(ax, m)"
84+
]
85+
},
86+
{
87+
"cell_type": "markdown",
88+
"metadata": {},
89+
"source": [
90+
"We will select wavelength of 300 km\n\nLow frequency\n\n"
91+
]
92+
},
93+
{
94+
"cell_type": "code",
95+
"execution_count": null,
96+
"metadata": {
97+
"collapsed": false
98+
},
99+
"outputs": [],
100+
"source": [
101+
"ax = start_axes(\"ADT low frequency\")\ng.copy(\"adt\", \"adt_low\")\ng.bessel_low_filter(\"adt_low\", 300, order=3)\nm = g.display(ax, \"adt_low\", vmin=-0.15, vmax=0.15)\nupdate_axes(ax, m)"
102+
]
103+
},
104+
{
105+
"cell_type": "markdown",
106+
"metadata": {},
107+
"source": [
108+
"High frequency\n\n"
109+
]
110+
},
111+
{
112+
"cell_type": "code",
113+
"execution_count": null,
114+
"metadata": {
115+
"collapsed": false
116+
},
117+
"outputs": [],
118+
"source": [
119+
"ax = start_axes(\"ADT high frequency\")\ng.copy(\"adt\", \"adt_high\")\ng.bessel_high_filter(\"adt_high\", 300, order=3)\nm = g.display(ax, \"adt_high\", vmin=-0.15, vmax=0.15)\nupdate_axes(ax, m)"
120+
]
121+
},
122+
{
123+
"cell_type": "markdown",
124+
"metadata": {},
125+
"source": [
126+
"Clues\n-----\nwavelength : 80km\n\n"
127+
]
128+
},
129+
{
130+
"cell_type": "code",
131+
"execution_count": null,
132+
"metadata": {
133+
"collapsed": false
134+
},
135+
"outputs": [],
136+
"source": [
137+
"g.copy(\"adt\", \"adt_high_80\")\ng.bessel_high_filter(\"adt_high_80\", 80, order=3)\ng.copy(\"adt\", \"adt_low_80\")\ng.bessel_low_filter(\"adt_low_80\", 80, order=3)\n\narea = dict(llcrnrlon=11.75, urcrnrlon=21, llcrnrlat=33, urcrnrlat=36.75)"
138+
]
139+
},
140+
{
141+
"cell_type": "markdown",
142+
"metadata": {},
143+
"source": [
144+
"Spectrum\n\n"
145+
]
146+
},
147+
{
148+
"cell_type": "code",
149+
"execution_count": null,
150+
"metadata": {
151+
"collapsed": false
152+
},
153+
"outputs": [],
154+
"source": [
155+
"fig = plt.figure(figsize=(10, 6))\nax = fig.add_subplot(111)\nax.set_title(\"Spectrum\")\nax.set_xlabel(\"km\")\n\nlon_spec, lat_spec = g.spectrum_lonlat(\"adt\", area=area)\nmappable = ax.loglog(*lat_spec, label=\"lat raw\")[0]\nax.loglog(*lon_spec, label=\"lon raw\", color=mappable.get_color(), linestyle=\"--\")\n\nlon_spec, lat_spec = g.spectrum_lonlat(\"adt_high_80\", area=area)\nmappable = ax.loglog(*lat_spec, label=\"lat high\")[0]\nax.loglog(*lon_spec, label=\"lon high\", color=mappable.get_color(), linestyle=\"--\")\n\nlon_spec, lat_spec = g.spectrum_lonlat(\"adt_low_80\", area=area)\nmappable = ax.loglog(*lat_spec, label=\"lat low\")[0]\nax.loglog(*lon_spec, label=\"lon low\", color=mappable.get_color(), linestyle=\"--\")\n\nax.set_xlim(10, 1000)\nax.set_ylim(1e-6, 1)\nax.set_xscale(\"log\")\nax.legend()\nax.grid()"
156+
]
157+
},
158+
{
159+
"cell_type": "markdown",
160+
"metadata": {},
161+
"source": [
162+
"Spectrum ratio\n\n"
163+
]
164+
},
165+
{
166+
"cell_type": "code",
167+
"execution_count": null,
168+
"metadata": {
169+
"collapsed": false
170+
},
171+
"outputs": [],
172+
"source": [
173+
"fig = plt.figure(figsize=(10, 6))\nax = fig.add_subplot(111)\nax.set_title(\"Spectrum ratio\")\nax.set_xlabel(\"km\")\n\nlon_spec, lat_spec = g.spectrum_lonlat(\n \"adt_high_80\", area=area, ref=g, ref_grid_name=\"adt\"\n)\nmappable = ax.plot(*lat_spec, label=\"lat high\")[0]\nax.plot(*lon_spec, label=\"lon high\", color=mappable.get_color(), linestyle=\"--\")\n\nlon_spec, lat_spec = g.spectrum_lonlat(\n \"adt_low_80\", area=area, ref=g, ref_grid_name=\"adt\"\n)\nmappable = ax.plot(*lat_spec, label=\"lat low\")[0]\nax.plot(*lon_spec, label=\"lon low\", color=mappable.get_color(), linestyle=\"--\")\n\nax.set_xlim(10, 1000)\nax.set_ylim(0, 1)\nax.set_xscale(\"log\")\nax.legend()\nax.grid()"
174+
]
175+
},
176+
{
177+
"cell_type": "markdown",
178+
"metadata": {},
179+
"source": [
180+
"Old filter\n----------\nTo do ...\n\n"
181+
]
182+
}
183+
],
184+
"metadata": {
185+
"kernelspec": {
186+
"display_name": "Python 3",
187+
"language": "python",
188+
"name": "python3"
189+
},
190+
"language_info": {
191+
"codemirror_mode": {
192+
"name": "ipython",
193+
"version": 3
194+
},
195+
"file_extension": ".py",
196+
"mimetype": "text/x-python",
197+
"name": "python",
198+
"nbconvert_exporter": "python",
199+
"pygments_lexer": "ipython3",
200+
"version": "3.7.7"
201+
}
202+
},
203+
"nbformat": 4,
204+
"nbformat_minor": 0
205+
}

‎src/py_eddy_tracker/data/__init__.py‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@
44
import tarfile
55
import lzma
66

7+
78
def get_path(name):
89
return path.join(path.dirname(__file__), name)
910

‎src/py_eddy_tracker/dataset/grid.py‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1471,7 +1471,8 @@ def spectrum_lonlat(self, grid_name, area=None, ref=None, **kwargs):
14711471
pws.append(pw)
14721472
if nb_invalid:
14731473
logger.warning("%d/%d columns invalid", nb_invalid, i + 1)
1474-
lat_content = 1 / f, array(pws).mean(axis=0)
1474+
with errstate(divide="ignore"):
1475+
lat_content = 1 / f, array(pws).mean(axis=0)
14751476

14761477
# Lon spectrum
14771478
fs, pws = list(), list()
@@ -1500,7 +1501,8 @@ def spectrum_lonlat(self, grid_name, area=None, ref=None, **kwargs):
15001501
for f, pw in zip(fs, pws)
15011502
]
15021503
).mean(axis=0)
1503-
lon_content = 1 / f_interp, pw_m
1504+
with errstate(divide="ignore"):
1505+
lon_content = 1 / f_interp, pw_m
15041506
if ref is None:
15051507
return lon_content, lat_content
15061508
else:

0 commit comments

Comments
 (0)