对任意信号时空滤波

fft_filter_area.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
64
65
66
67
68
69
"""
使用时空滤波提取感兴趣范围的信号。
本示例演示了如何提取 OLR 频谱中波数为 14，频率为 0.1 cpd 附近的极强信号
（实则为 OLR 插值卫星扫描数据所产生的假信号）。
"""
import matplotlib.pyplot as plt
import xarray as xr
import numpy as np
import moisten_ew as mew

# 加载数据，并选择南北纬20度以内
data = xr.open_dataset("~/Documents/equator_waves/data/olr.day.mean.nc").olr.sel(
    lat=slice(20, -20),
    time=slice("2000", "2020")
)

# 计算WK99频谱，用于画图参考滤波区域
res = mew.wk99_spectrum(data, 'lon', 'lat', 'time')
res = res.sel(wavenumber=slice(-20, 20), freq=slice(0, 1))


# 设置一个椭圆滤波区域
area = mew.EllipseFilterArea(
    center_wavenumber=14.15 * mew.unit('zonal_wavenumber'),
    center_frequency=0.097 * mew.unit('cpd'),
    wavenumber_radius=1 * mew.unit('zonal_wavenumber'),
    freq_radius=0.02 * mew.unit('cpd')
)

# 计算滤波结果
filtered = mew.fft_time_lon_area_filter(
    data.sel(time='2000-11'), area,
    lon_name='lon', time_name='time'
)


# 画图
fig = plt.figure(figsize=(12, 4))
gs = fig.add_gridspec(1, 3, wspace=0.4, right=0.96, left=0.08)
for i in range(2):
    ax = fig.add_subplot(gs[0, i])

    # 画与背景的比值
    d = res.spec_ratio.isel(type=i)

    c = ax.contourf(d.wavenumber, d.freq, np.log(d),
                cmap='RdBu_r', levels=np.linspace(-1, 1, 21),
                extend='both')

    # 画滤波区域
    ax.plot(area.shape.exterior.xy[0], area.shape.exterior.xy[1],
            color='r', lw=1, label='Filter Area')

    ax.set_ylim(0, 0.4)
    ax.axvline(0, color='k', lw=0.5, ls='--')
    ax.set_xlabel('Zonal Wavenumber')
    ax.set_ylabel('Frequency (cpd)')
    ax.set_title(f'{d.type.values}', loc='left')
    ax.legend()

ax = fig.add_subplot(gs[0, 2])
im = ax.pcolormesh(
    filtered.lon, filtered.time, filtered.mean(dim='lat'),
    cmap='jet'
)
ax.set_xlabel('Longitude')
ax.set_title('OLR Filter result', loc='left')

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