北大西洋涛动（North Atlantic oscillation）
北大西洋涛动（ NAO ) 是北大西洋上空冰岛低压和亚速尔群岛高压海平面气压差（SLP）波动的一种天气现象。它是由吉尔伯特·沃克爵士 (Sir Gilbert Walker) 在上世纪20年代发现的，他此前曾发现太平洋中的南方涛动(SO) ，而 厄尔尼诺-南方涛动 (ENSO) 现象与该现象有关。与 ENSO 一样，NAO 是区域范围内的决定性气候因素，通过冰岛低压和亚速尔高压强度的波动，它控制着西风的强度和方向以及横跨北大西洋的风暴路径(Storm track) 的位置。
北大西洋涛动与北极涛动(AO)（或北方环模 (NAM)） 密切相关，但不应与大西洋多年代际涛动(AMO) 混淆。
NAO定义：
NAO有多种定义。最容易理解的是基于测量站点之间季节性平均气压差，例如：
里斯本(Lisbon) 和 斯蒂基松鲁姆/雷克雅未克（Stykkishólmur/Reykjavík） 蓬塔德尔加达(Ponta Delgada)、亚速尔群岛和斯蒂基松鲁姆/雷克雅未克 亚速尔群岛 (1865–2002)、直布罗陀(1821–2007) 和雷克雅未克
共同点：冰岛的北点（该地区唯一有长期记录的站点）；以及各个南方点。所有这些都试图通过选择亚速尔群岛高压区和冰岛低压区（如上图所示）这两个稳定气压区中的监测站来捕捉相同的变化模式。
更复杂的定义基于海平面气压的经验正交函数（EOF）或者旋转经验正交函数（REOF）。
NAO相位:
负NAO-相位 ：海平面气压负异常从亚速尔群岛延伸至西地中海，正异常从格陵兰岛延伸至其北部海域。斜压区不太明显且更南。 正NAO+相位 ：副热带地区海平面气压正异常，高纬地区负异常。大西洋和西欧西风纬向环流加强，斜压区向北移动。
赫雷尔NAO指数(Hurrell NAO Index)
定义中基于里斯本和葡萄牙里斯本和冰岛斯蒂基斯霍尔米/雷克雅未克之间标准化海平面压力 (SLP) 的差异计算的NAO指数也被称为赫雷尔NAO指数。
该指数的优点：
基于站点的指数可以追溯到19世纪中叶或更早。 容易计算和理解。
缺点：
这些站点在空间上是固定的，因此可能无法跟踪NAO中心在年度周期中的移动情况。 由于与NAO无关的小规模和瞬态天气现象，导致各个观测站点的海平面气压可能会有噪声。
由于没有特定的方法来定义 NAO 的空间结构，因此也没有普遍接受的指标来描述该现象的时间演变。大多数现代 NAO 指数要么源自南北不同地点之间海平面气压异常的简单差异，要么源自SLP（通常是区域性）的EOF 的 PC 时间序列。前者的例子有很多，通常基于 NAO 中心附近各个站的仪器记录，但有时来自网格化的SLP分析。
采用PC时间序列方法的一个优点是，此类指数能更好地反映整个NAO空间模式；然而，由于这些指数基于网格化的 SLP 数据，因此只能计算 20 世纪至今的部分时间，具体取决于数据来源。
Python计算NAO以及比较
冬季（DJFM）赫雷尔NAO指数可以从以下网站下载：https://climatedataguide.ucar.edu/sites/default/files/2023-07/nao_station_djfm.txt。
ERA5的海平面气压数据已经使用CDO(https://code.mpimet.mpg.de/projects/cdo/)进行处理，以便仅保留北大西洋地区DJFM月份的平均值：
保留北大西洋地区：
cdo sellonlatbox,-100,40,20,80 era5.mslp.mon.mean.nc era5.mslp.mon.mean.atl.nc
保留DJFM月份：
cdo timselmean,4,11,8 era5.mslp.mon.mean.atl.nc era5.mslp.djfm.mean.atl.nc
绘制ERA5海平面气压气候态
import os
import numpy as np
import xarray as xr
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import matplotlib.path as mpath
from eofs.standard import Eof
Lisbon = [-9, 39]
Stykkisholmur = [-22.5, 65]
def plot_background(ax):
ax.coastlines()
ax.gridlines(crs=ccrs.PlateCarree(), linewidth=0.5, color='gray', alpha=0.5, linestyle='-')
boundary_path = make_boundary_path(lon, lat)
ax.set_boundary(boundary_path, transform=ccrs.PlateCarree())
return ax
def make_boundary_path(lon,lat):
lons,lats=np.meshgrid(lon,lat)
boundary_path = np.array([lons[-1,:],lats[-1,:]])
boundary_path = np.append(boundary_path,np.array([lons[::-1,-1],lats[::-1,-1]]),axis=1)
boundary_path = np.append(boundary_path,np.array([lons[1,::-1],lats[1,::-1]]),axis=1)
boundary_path = np.append(boundary_path,np.array([lons[:,1],lats[:,1]]),axis=1)
boundary_path = mpath.Path(np.swapaxes(boundary_path, 0, 1))
return boundary_path
def plot_town(ax):
ax.plot(Lisbon[0], Lisbon[1], marker='o', color='red', markersize=8,
alpha=0.7, transform=ccrs.PlateCarree())
ax.text(Lisbon[0], Lisbon[1], u'Lisbon',
verticalalignment='center', horizontalalignment='right',
transform=ccrs.PlateCarree()._as_mpl_transform(ax),
bbox=dict(facecolor='sandybrown', alpha=0.5, boxstyle='round'))
ax.plot(Stykkisholmur[0], Stykkisholmur[1], marker='o', color='red', markersize=8,
alpha=0.7, transform=ccrs.PlateCarree())
ax.text(Stykkisholmur[0], Stykkisholmur[1], u'Stykkisholmur',
verticalalignment='center', horizontalalignment='right',
transform=ccrs.PlateCarree()._as_mpl_transform(ax),
bbox=dict(facecolor='sandybrown', alpha=0.5, boxstyle='round'))
def plot_pc(ax):
plt.xlabel('Year')
plt.axhline(0, color='k')
plt.xlim(int(yr1), int(yr2)+1)
plt.xticks(xi, range(int(yr1), int(yr2)+1,5))
ax.axvline(1950, color='grey', linestyle='--', linewidth=0.5)
ax.axvline(1960, color='grey', linestyle='--', linewidth=0.5)
ax.axvline(1970, color='grey', linestyle='--', linewidth=0.5)
ax.axvline(1980, color='grey', linestyle='--', linewidth=0.5)
ax.axvline(1990, color='grey', linestyle='--', linewidth=0.5)
ax.axvline(2000, color='grey', linestyle='--', linewidth=0.5)
ax.axvline(2010, color='grey', linestyle='--', linewidth=0.5)
plt.ylim(-1500, 1500)
return ax
dir_data='./data/'
dir_figs='./figs/'
if not os.path.exists(dir_figs):
os.makedirs(dir_figs)
nc_file = dir_data+'era5.mslp.djfm.mean.natl.nc'
data = xr.open_dataset(nc_file)
latS=20
latN=80
lonW=-80
lonE=30
yr1 = '1981'
yr2 = '2020'
years = np.arange(int(yr1), int(yr2)+1)
xi = [i for i in range(int(yr1), int(yr2)+1,5)]
data = xr.open_dataset(nc_file).sel(latitude=slice(latN,latS)).sel(longitude=slice(lonW,lonE)).sel(time=slice(yr1,yr2))
data_clim = data.mean('time')
print(data)
print(data_clim)
lat = data.latitude.values
lon = data.longitude.values
time = data.time.values
slp, slp_clim, = xr.broadcast(data['msl']/100, data_clim['msl']/100)
#-- Manage dates
time_str=[x for x in range(len(time))]
date_str=[x for x in range(len(time))]
for i in range(len(time)):
time_str[i] = str(time[i])
date_str[i] = time_str[i][0:10]
slp_anom=(slp-slp_clim)
levels = np.arange(980,1042,2)
levels_a = np.arange(-16,18,2)
projection = ccrs.Orthographic(central_longitude=(lonW+lonE)/2, central_latitude=(latS+latN)/2)
fig = plt.figure(figsize=(15., 10.), dpi = 300)
ax = fig.add_subplot(1, 1, 1, projection=projection)
ax.set_title('Mean Sea Level Pressure (hPa) - DJFM '+ yr1+'-'+yr2, loc='center')
plot_background(ax)
cf = ax.contourf(lon, lat, slp_clim[0,:,:], levels, cmap='jet', transform=ccrs.PlateCarree(), extend='both')
c = ax.contour(lon, lat, slp_clim[0,:,:], levels, colors='black', linewidths=0.2, transform=ccrs.PlateCarree())
cb = fig.colorbar(cf, orientation='horizontal', aspect=65, shrink=0.5, pad=0.05, extendrect='True')
cb.set_label('hPa', size='large')
plot_town(ax)
plt.show()
计算EOF前四个模态及其贡献率
wgts = np.sqrt(np.cos(np.deg2rad(lat)))[:, np.newaxis]
solver = Eof(np.array(slp_anom), weights=wgts, center=True)
eigenvalues = solver.eigenvalues()
total_variance = solver.totalAnomalyVariance()
varfrac = solver.varianceFraction()
eofs = solver.eofs()
pcs = solver.pcs()
pcs_norm = solver.pcs(pcscaling=1)
clevs = np.linspace(-0.005, 0.005, 21)
绘制SLP EOF第一模态以及PC1时间序列
fig = plt.figure(figsize=(20, 8), dpi = 300)
fig.suptitle('EOF1 and PC1 : MSLP DJFM '+yr1+'-'+yr2, fontsize=16)
ax = fig.add_subplot(121, projection=projection)
plt.title('EOF1 ('+str(int(varfrac[0]*100))+'%)', fontsize=10, loc='center')
plot_background(ax)
plot_town(ax)
cf = ax.contourf(lon, lat, eofs[0], levels=clevs, transform=ccrs.PlateCarree(), cmap='RdBu_r', extend='both')
c = ax.contour(lon, lat, eofs[0], levels=clevs, colors='black', linewidths=1, transform=ccrs.PlateCarree())
cb = fig.colorbar(cf, orientation='horizontal', aspect=65, shrink=1, pad=0.05, extendrect='True')
ax = fig.add_subplot(122)
plt.ylabel('PC1')
plot_pc(ax)
colormat=np.where(pcs[:,0]>0, 'red','blue')
plt.bar(years, pcs[:,0], width=1, color=colormat, edgecolor = 'k')
plt.show()
比较赫雷尔NAO指数以及PC1时间序列
fig = plt.figure(figsize=(15., 8.), dpi = 300)
nao_file=dir_data+'nao_station_djfm_0.txt'
X, Y = [], []
for line in open(nao_file, 'r'):
values = [float(s) for s in line.split()]
X.append(values[0])
Y.append(values[1])
pc1=pcs_norm[:,0]
XX=np.asarray(X)
YY=np.asarray(Y)
id=np.where((XX >= int(yr1)) & (XX <= int(yr2)))
nao_year=XX[id]
nao_index=YY[id]
cor=np.corrcoef(pc1, nao_index)
plt.title('Normalized PC1 VS NAO station-based index : '+yr1+'-'+yr2, loc='center')
plt.title('Correlation : '+str(round(cor[0,1],2)), loc='right')
plt.xlabel('Year')
plt.ylabel('Normalized Units')
yrs = range(int(yr1), int(yr2)+1)
xi = [i for i in range(int(yr1), int(yr2)+1,2)]
plt.plot(yrs, pc1, color='b', linewidth=2, label='PC1 from EOF analysis')
plt.plot(yrs, nao_index, color='r', linewidth=2, label='NAO index')
plt.axhline(0, color='k')
plt.axhline(2, color='grey', linestyle='--', linewidth=0.5)
plt.axhline(-2, color='grey', linestyle='--', linewidth=0.5)
plt.xlim(int(yr1), int(yr2)+1)
plt.xticks(xi, range(int(yr1), int(yr2)+1,2))
plt.ylim(-6, 6)
locs, labels = plt.xticks()
plt.setp(labels, rotation=90)
plt.legend()
plt.show()
两者相关系数达到0.96。
Github 仓库
代码和数据已上传至Github仓库，可一键跑通：https://github.com/SQYQianYe/NAO-Index
完整代码
import os
import numpy as np
import xarray as xr
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import matplotlib.path as mpath
from eofs.standard import Eof
Lisbon = [-9, 39]
Stykkisholmur = [-22.5, 65]
def plot_background(ax):
ax.coastlines()
ax.gridlines(crs=ccrs.PlateCarree(), linewidth=0.5, color='gray', alpha=0.5, linestyle='-')
boundary_path = make_boundary_path(lon, lat)
ax.set_boundary(boundary_path, transform=ccrs.PlateCarree())
return ax
def make_boundary_path(lon,lat):
lons,lats=np.meshgrid(lon,lat)
boundary_path = np.array([lons[-1,:],lats[-1,:]])
boundary_path = np.append(boundary_path,np.array([lons[::-1,-1],lats[::-1,-1]]),axis=1)
boundary_path = np.append(boundary_path,np.array([lons[1,::-1],lats[1,::-1]]),axis=1)
boundary_path = np.append(boundary_path,np.array([lons[:,1],lats[:,1]]),axis=1)
boundary_path = mpath.Path(np.swapaxes(boundary_path, 0, 1))
return boundary_path
def plot_town(ax):
ax.plot(Lisbon[0], Lisbon[1], marker='o', color='red', markersize=8,
alpha=0.7, transform=ccrs.PlateCarree())
ax.text(Lisbon[0], Lisbon[1], u'Lisbon',
verticalalignment='center', horizontalalignment='right',
transform=ccrs.PlateCarree()._as_mpl_transform(ax),
bbox=dict(facecolor='sandybrown', alpha=0.5, boxstyle='round'))
ax.plot(Stykkisholmur[0], Stykkisholmur[1], marker='o', color='red', markersize=8,
alpha=0.7, transform=ccrs.PlateCarree())
ax.text(Stykkisholmur[0], Stykkisholmur[1], u'Stykkisholmur',
verticalalignment='center', horizontalalignment='right',
transform=ccrs.PlateCarree()._as_mpl_transform(ax),
bbox=dict(facecolor='sandybrown', alpha=0.5, boxstyle='round'))
def plot_pc(ax):
plt.xlabel('Year')
plt.axhline(0, color='k')
plt.xlim(int(yr1), int(yr2)+1)
plt.xticks(xi, range(int(yr1), int(yr2)+1,5))
ax.axvline(1950, color='grey', linestyle='--', linewidth=0.5)
ax.axvline(1960, color='grey', linestyle='--', linewidth=0.5)
ax.axvline(1970, color='grey', linestyle='--', linewidth=0.5)
ax.axvline(1980, color='grey', linestyle='--', linewidth=0.5)
ax.axvline(1990, color='grey', linestyle='--', linewidth=0.5)
ax.axvline(2000, color='grey', linestyle='--', linewidth=0.5)
ax.axvline(2010, color='grey', linestyle='--', linewidth=0.5)
plt.ylim(-1500, 1500)
return ax
dir_data='./data/'
dir_figs='./figs/'
if not os.path.exists(dir_figs):
os.makedirs(dir_figs)
nc_file = dir_data+'era5.mslp.djfm.mean.natl.nc'
data = xr.open_dataset(nc_file)
latS=20
latN=80
lonW=-80
lonE=30
yr1 = '1981'
yr2 = '2020'
years = np.arange(int(yr1), int(yr2)+1)
xi = [i for i in range(int(yr1), int(yr2)+1,5)]
data = xr.open_dataset(nc_file).sel(latitude=slice(latN,latS)).sel(longitude=slice(lonW,lonE)).sel(time=slice(yr1,yr2))
data_clim = data.mean('time')
print(data)
print(data_clim)
lat = data.latitude.values
lon = data.longitude.values
time = data.time.values
slp, slp_clim, = xr.broadcast(data['msl']/100, data_clim['msl']/100)
#-- Manage dates
time_str=[x for x in range(len(time))]
date_str=[x for x in range(len(time))]
for i in range(len(time)):
time_str[i] = str(time[i])
date_str[i] = time_str[i][0:10]
slp_anom=(slp-slp_clim)
levels = np.arange(980,1042,2)
levels_a = np.arange(-16,18,2)
projection = ccrs.Orthographic(central_longitude=(lonW+lonE)/2, central_latitude=(latS+latN)/2)
fig = plt.figure(figsize=(15., 10.), dpi = 300)
ax = fig.add_subplot(1, 1, 1, projection=projection)
ax.set_title('Mean Sea Level Pressure (hPa) - DJFM '+ yr1+'-'+yr2, loc='center')
plot_background(ax)
cf = ax.contourf(lon, lat, slp_clim[0,:,:], levels, cmap='jet', transform=ccrs.PlateCarree(), extend='both')
c = ax.contour(lon, lat, slp_clim[0,:,:], levels, colors='black', linewidths=0.2, transform=ccrs.PlateCarree())
cb = fig.colorbar(cf, orientation='horizontal', aspect=65, shrink=0.5, pad=0.05, extendrect='True')
cb.set_label('hPa', size='large')
plot_town(ax)
plt.show()
wgts = np.sqrt(np.cos(np.deg2rad(lat)))[:, np.newaxis]
solver = Eof(np.array(slp_anom), weights=wgts, center=True)
eigenvalues = solver.eigenvalues()
total_variance = solver.totalAnomalyVariance()
varfrac = solver.varianceFraction()
eofs = solver.eofs()
pcs = solver.pcs()
pcs_norm = solver.pcs(pcscaling=1)
clevs = np.linspace(-0.005, 0.005, 21)
fig = plt.figure(figsize=(20, 8), dpi = 300)
fig.suptitle('EOF1 and PC1 : MSLP DJFM '+yr1+'-'+yr2, fontsize=16)
ax = fig.add_subplot(121, projection=projection)
plt.title('EOF1 ('+str(int(varfrac[0]*100))+'%)', fontsize=10, loc='center')
plot_background(ax)
plot_town(ax)
cf = ax.contourf(lon, lat, eofs[0], levels=clevs, transform=ccrs.PlateCarree(), cmap='RdBu_r', extend='both')
c = ax.contour(lon, lat, eofs[0], levels=clevs, colors='black', linewidths=1, transform=ccrs.PlateCarree())
cb = fig.colorbar(cf, orientation='horizontal', aspect=65, shrink=1, pad=0.05, extendrect='True')
ax = fig.add_subplot(122)
plt.ylabel('PC1')
plot_pc(ax)
colormat=np.where(pcs[:,0]>0, 'red','blue')
plt.bar(years, pcs[:,0], width=1, color=colormat, edgecolor = 'k')
plt.show()
fig = plt.figure(figsize=(15., 8.), dpi = 300)
nao_file=dir_data+'nao_station_djfm_0.txt'
X, Y = [], []
for line in open(nao_file, 'r'):
values = [float(s) for s in line.split()]
X.append(values[0])
Y.append(values[1])
pc1=pcs_norm[:,0]
XX=np.asarray(X)
YY=np.asarray(Y)
id=np.where((XX >= int(yr1)) & (XX <= int(yr2)))
nao_year=XX[id]
nao_index=YY[id]
cor=np.corrcoef(pc1, nao_index)
plt.title('Normalized PC1 VS NAO station-based index : '+yr1+'-'+yr2, loc='center')
plt.title('Correlation : '+str(round(cor[0,1],2)), loc='right')
plt.xlabel('Year')
plt.ylabel('Normalized Units')
yrs = range(int(yr1), int(yr2)+1)
xi = [i for i in range(int(yr1), int(yr2)+1,2)]
plt.plot(yrs, pc1, color='b', linewidth=2, label='PC1 from EOF analysis')
plt.plot(yrs, nao_index, color='r', linewidth=2, label='NAO index')
plt.axhline(0, color='k')
plt.axhline(2, color='grey', linestyle='--', linewidth=0.5)
plt.axhline(-2, color='grey', linestyle='--', linewidth=0.5)
plt.xlim(int(yr1), int(yr2)+1)
plt.xticks(xi, range(int(yr1), int(yr2)+1,2))
plt.ylim(-6, 6)
locs, labels = plt.xticks()
plt.setp(labels, rotation=90)
plt.legend()
plt.show()
References
https://www.ncei.noaa.gov/access/monitoring/nao/
