对 MRG 与 Kelvin 波时空滤波

fft_time_lon.png

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
"""使用FFT 时空滤波，提取 MRG 波与 Kelvin 波的信号"""

import numpy as np
import xarray as xr
import moisten_ew as mew
import matplotlib.pyplot as plt


# 加载数据
data = xr.open_dataset("~/Documents/equator_waves/data/olr.day.mean.nc").olr.sel(
    lat=slice(10, -10), time="2000-11"
)

# ==== 时空滤波 ====

# 设置 MRGW 滤波区域
mrg_area = mew.MRGFilterArea() # 使用默认参数
# 滤波
mrg_data = mew.fft_time_lon_area_filter(
    data, mrg_area,
    lon_name='lon', time_name='time'
)

# 设置 Kelvin 波滤波区域
kelvin_area = mew.KelvinFilterArea() # 使用默认参数
# 滤波
kelvin_data = mew.fft_time_lon_area_filter(
    data, kelvin_area,
    lon_name='lon', time_name='time'
)

# ==== 画图 ====
fig = plt.figure(figsize=(10, 4), constrained_layout=True)
gs = fig.add_gridspec(1, 3, hspace=0.3)

# 画滤波区域
ax = fig.add_subplot(gs[0, 0])

# 区域为 shapely Polygon 对象，画出其外边界
ax.plot(mrg_area.shape.exterior.xy[0], mrg_area.shape.exterior.xy[1],
        label='MRG Wave Area')
ax.plot(kelvin_area.shape.exterior.xy[0], kelvin_area.shape.exterior.xy[1],
        label='Kelvin Wave Area')

ax.set_xlabel('Zonal Wavenumber')
ax.set_ylabel('Frequency (cpd)')
ax.set_xlim(-20, 20)
ax.set_ylim(0, 0.5)
ax.axvline(0, color='k', lw=0.5, ls='--')
ax.set_title('Filter Areas', loc='left')
ax.legend()

# 画滤波后的数据
for i, wave_data in enumerate([mrg_data, kelvin_data]):
    ax = fig.add_subplot(gs[0, i+1])
    im = ax.pcolormesh(
        wave_data.lon, wave_data.time, wave_data.mean(dim='lat'),
        cmap='jet'
    )
    ax.set_xlabel('Longitude')
    ax.set_title(['MGR Wave', 'Kelvin Wave'][i], loc='left')

fig.savefig("docs/images/fft_time_lon.png", dpi=300)