前言
最近有粉丝在后台私信问我："空间场，海温场这些如何滤波？","一直以为时间序列才需要滤波"。很久之前做过相关的工作，今天就分享给大家。
高通滤波
今天我们以高通滤波举例，高通滤波（High-pass filtering） 是一种信号处理技术，用于移除信号中频率较低的成分（即长时间尺度的波动），保留频率较高的成分（即短时间尺度的波动）。
对于海温等变量，存在陆地掩膜或者缺测值填充的情况，我们今天就以sst变量举例，对三维变量场的高通滤波的代码如下：
def high_pass_filter(data: xr.DataArray | np.array,
cutoff: np.array
) -> np.array:
"""
Apply a high-pass filter to remove frequencies lower than the cutoff.
Parameters:
- data: Input data array (time, lat, lon)
- cutoff: Cutoff period in years
Returns:
- Filtered data array
"""
if isinstance(data, xr.DataArray):
data = data.values
fft_data = fft(data, axis=0)
freqs = np.fft.fftfreq(data.shape[0])
high_pass_mask = np.abs(freqs) > 1 / cutoff
filtered_data = ifft(fft_data * high_pass_mask[:, None, None], axis=0)
return np.real(filtered_data)
对于以上代码的检验，我们今天以论文复现的方式来检验，请大家继续往下看。
论文复现
今天复现的论文是《秋季北大西洋马蹄型海温异常与初冬我国气温年际变化的联系》。
该论文主要关注秋季北大西洋马蹄型海温异常对 冬季NAO 和 对流层亚洲极涡等大气环流异常的影响。
论文中的图(1)如下：
该图是我们今天需要复现的图，接下来我们来看数据在论文中的如何处理的：
英国气象局哈德莱中心（The UK Meteorological Office Hadley Center）提供的月平均海表温度（sea surface temperature，SST）资料（Rayner et al., 2003），水平分辨率为 1°×1°，时间段为1951～2020 年，数据下载地址： https://www.metoffice.gov.uk/hadobs/hadisst/.文中秋季指的是 9 月、10 月和 11 月三个月平均，初冬指的是 12 月。 研究主要关注年际变率，因此所有资料均去除 9 年以上的年代际变率和线性趋势。
有了以上处理信息，那么就很简单了。从文中的图片风格来看，作者使用了NCL来绘图，下面我们使用Python来复现这张图。
Python代码复现
前期准备
import numpy as np
import xarray as xr
from eofs.standard import Eof
from matplotlib import pyplot as plt
import cartopy.mpl.ticker as cticker
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import os
from scipy.fftpack import fft, ifft
def high_pass_filter(data: xr.DataArray | xr.Dataset,
cutoff: np.array
) -> np.array:
"""
Apply a high-pass filter to remove frequencies lower than the cutoff.
Parameters:
- data: Input data array (time, lat, lon)
- cutoff: Cutoff period in years
Returns:
- Filtered data array
"""
if isinstance(data, xr.DataArray):
data = data.values
fft_data = fft(data, axis=0)
freqs = np.fft.fftfreq(data.shape[0])
high_pass_mask = np.abs(freqs) > 1 / cutoff
filtered_data = ifft(fft_data * high_pass_mask[:, None, None], axis=0)
return np.real(filtered_data)
读取数据
Extend = [-80,0,10,70]
Lat_Range = slice(70, 10)
Lon_Range = slice(-80, 0)
d = xr.open_dataset('HadISST_sst.nc').sel(latitude = Lat_Range, longitude = Lon_Range, time = slice("1951", "2020"))
lat = d.latitude
lon = d.longitude
sst = d.sst
9年高通滤波
cutoff_period = 9 * 12
sst = high_pass_filter(sst, cutoff_period)
去除季节循环和去除全球趋势
去除季节循环是为了：
去除季节性影响，突出长期变化 避免季节性变化主导 EOF 模式 降低噪声，提高鲁棒性
sst = xr.DataArray(sst,dims=['time','latitude','longitude'], coords=[d.time, lat, lon])
Meansst= sst.groupby("time.month").mean()
ssta = (sst.groupby("time.month") - Meansst).where(lambda m: m.time.dt.month.isin([9, 10, 11])).groupby('time.year').mean("time")
ssta = ssta - ssta.weighted(np.cos(np.deg2rad(lat))).mean(('latitude','longitude'))
EOF
lat = np.array(lat)
coslat = np.cos(np.deg2rad(lat))
wgts = np.sqrt(coslat)[..., np.newaxis]
ssta = np.array(ssta)
solver = Eof(ssta, weights=wgts,center = True)
eof1 = solver.eofsAsCorrelation().squeeze()
pc1 = solver.pcs(pcscaling=1).squeeze()
var = solver.varianceFraction()[0]
pc1 = pc1[:,0]
绘图
plt.figure(figsize=(10,4),facecolor = 'w',dpi = 300)
plt.subplots_adjust(left=0.05, right=0.98, top=0.9, bottom=0.1, wspace=0.1)
proj = ccrs.PlateCarree()
leftlon, rightlon, lowerlat, upperlat = (Extend[0],Extend[1],Extend[2],Extend[3])
lon_formatter = cticker.LongitudeFormatter()
lat_formatter = cticker.LatitudeFormatter()
ax = plt.subplot(121, projection=ccrs.PlateCarree())
ax.set_extent([leftlon, rightlon, lowerlat, upperlat], crs=ccrs.PlateCarree())
ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax.add_feature(cfeature.LAKES, alpha=0.5)
ax.set_xticks(np.arange(leftlon,rightlon+10,20), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(lowerlat,upperlat+10,20), crs=ccrs.PlateCarree())
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
C = plt.contourf(lon,lat, eof1[0,:,:], levels=np.arange(-1,1.01,0.05), transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)
Cf = plt.contour(lon,lat, eof1[0,:,:], levels=np.arange(-1,1.01,0.2),
colors = 'k',transform=ccrs.PlateCarree())
ax.coastlines(color="#888888",lw = 2)
plt.title('(a) EOF1',loc='left',fontsize =15)
plt.title( '%.2f%%' % (var*100),loc='right',fontsize =15)
cax = plt.axes([0.06, -0.02, 0.42, 0.03])
plt.colorbar(C, cax =cax, orientation="horizontal", ticks = [-1, -0.5, 0, 0.5, 1])
ax2 = plt.subplot(122)
plt.title('(b) PC1',loc='left',fontsize =15)
rects2 = plt.bar(np.arange(1951,2021),np.where(pc1>0,pc1,np.nan),color='red',
width = 1,edgecolor='k'
)
rects3 = plt.bar(np.arange(1951,2021),np.where(pc1<0,pc1,np.nan),color='blue',
width = 1,edgecolor='k'
)
plt.xlim(1950, 2020)
plt.show()
Output
第一模态贡献率为20.72% ,论文中则是 20.5% ，可以说是相当接近了。
完整代码
import numpy as np
import xarray as xr
from eofs.standard import Eof
from matplotlib import pyplot as plt
import cartopy.mpl.ticker as cticker
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from scipy.fftpack import fft, ifft
def high_pass_filter(data: xr.DataArray | np.array,
cutoff: np.array
) -> np.array:
"""
Apply a high-pass filter to remove frequencies lower than the cutoff.
Parameters:
- data: Input data array (time, lat, lon)
- cutoff: Cutoff period in years
Returns:
- Filtered data array
"""
if isinstance(data, xr.DataArray):
data = data.values
fft_data = fft(data, axis=0)
freqs = np.fft.fftfreq(data.shape[0])
high_pass_mask = np.abs(freqs) > 1 / cutoff
filtered_data = ifft(fft_data * high_pass_mask[:, None, None], axis=0)
return np.real(filtered_data)
Extend = [-80,0,10,70]
Lat_Range = slice(70, 10)
Lon_Range = slice(-80, 0)
d = xr.open_dataset('HadISST_sst.nc').sel(latitude = Lat_Range, longitude = Lon_Range, time = slice("1951", "2020"))
lat = d.latitude
lon = d.longitude
sst = d.sst
# 9年高通滤波
cutoff_period = 9 * 12
sst = high_pass_filter(sst, cutoff_period)
# 去季节循环
sst = xr.DataArray(sst,dims=['time','latitude','longitude'], coords=[d.time, lat, lon])
Meansst= sst.groupby("time.month").mean()
ssta = (sst.groupby("time.month") - Meansst).where(lambda m: m.time.dt.month.isin([9, 10, 11])).groupby('time.year').mean("time")
ssta = ssta - ssta.weighted(np.cos(np.deg2rad(lat))).mean(('latitude','longitude'))
# EOF
lat = np.array(lat)
coslat = np.cos(np.deg2rad(lat))
wgts = np.sqrt(coslat)[..., np.newaxis]
ssta = np.array(ssta)
solver = Eof(ssta, weights=wgts,center = True)
eof1 = solver.eofsAsCorrelation().squeeze()
pc1 = solver.pcs(pcscaling=1).squeeze()
var = solver.varianceFraction()[0]
pc1 = pc1[:,0]
plt.figure(figsize=(10,4),facecolor = 'w',dpi = 300)
plt.subplots_adjust(left=0.05, right=0.98, top=0.9, bottom=0.1, wspace=0.1)
proj = ccrs.PlateCarree()
leftlon, rightlon, lowerlat, upperlat = (Extend[0],Extend[1],Extend[2],Extend[3])
lon_formatter = cticker.LongitudeFormatter()
lat_formatter = cticker.LatitudeFormatter()
ax = plt.subplot(121, projection=ccrs.PlateCarree())
ax.set_extent([leftlon, rightlon, lowerlat, upperlat], crs=ccrs.PlateCarree())
ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax.add_feature(cfeature.LAKES, alpha=0.5)
ax.set_xticks(np.arange(leftlon,rightlon+10,20), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(lowerlat,upperlat+10,20), crs=ccrs.PlateCarree())
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
C = plt.contourf(lon,lat, eof1[0,:,:], levels=np.arange(-1,1.01,0.05), transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)
Cf = plt.contour(lon,lat, eof1[0,:,:], levels=np.arange(-1,1.01,0.2),
colors = 'k',transform=ccrs.PlateCarree())
ax.coastlines(color="#888888",lw = 2)
plt.title('(a) EOF1',loc='left',fontsize =15)
plt.title( '%.2f%%' % (var*100),loc='right',fontsize =15)
cax = plt.axes([0.06, -0.02, 0.42, 0.03])
plt.colorbar(C, cax =cax, orientation="horizontal", ticks = [-1, -0.5, 0, 0.5, 1])
ax2 = plt.subplot(122)
plt.title('(b) PC1',loc='left',fontsize =15)
rects2 = plt.bar(np.arange(1951,2021),np.where(pc1>0,pc1,np.nan),color='red',
width = 1,edgecolor='k'
)
rects3 = plt.bar(np.arange(1951,2021),np.where(pc1<0,pc1,np.nan),color='blue',
width = 1,edgecolor='k'
)
plt.xlim(1950, 2020)
plt.show()
References
李忠贤, 王庭轩, 曾刚, 等. 2024. 秋季北大西洋马蹄型海温异常与初冬我国气温年际变化的联系 [J]. 大气科学, 48(3): 1131−1143. LI Zhongxian, WANG Tingxuan, ZENG Gang, et al. 2024. Characteristics of North Atlantic Horseshoe Sea Surface Temperature Anomaly in Autumn and Relationship with Interannual Variation in Early Winter Temperature in China [J]. Chinese Journal of Atmospheric Sciences (in Chinese), 48(3): 1131−1143. doi:10.3878/j.issn.1006-9895.2209.22106
往期回顾
Python | 解决cmaps库报错的问题 Python | 北大西洋涛动 | NAO指数 | EOF Python | HadSLP2 ASCII转NetCDF Python | 半球间经度转换 Python | 大气科学 | 偏相关 37个CMIP6模型对流参数化方案
