本文的 jupyter notebook 已上传至 github，欢送 fork 测试~
0. 导言
咱们都晓得，气象钻研的时间跨度个别都较长，根本都在 30 年以上，这就意味着对应的数据集非常宏大，既不能简略地对数据进行形容，也无奈轻易地从数据中提取特色。那么面对如此宏大的数据集，咱们如何能力从中提取出最能形容它的次要特色呢?
在这种场景下，EOF 合成就显示出它的弱小劣势了。它能够把随着工夫变动的气象要素场，合成为 空间函数局部 和 工夫函数局部，从而便于咱们发展剖析和钻研，让咱们可能从宏大的气象数据中抓住他们的次要特色。本文将介绍 EOF 的基本原理，以及在 python 中如何实现 EOF 合成。如有不妥之处，还请大家不吝赐教!
1. 什么是教训正交合成
教训正交合成 (Empirical orthogonal function，EOF，以下简称“EOF 合成”) 又被称为主成分剖析(Principal component analysis，PCA)，最早由 Pearson 提出，20 世纪 50 年代 Lorenz 将该办法引入大气迷信畛域，随后被广泛应用至今。
EOF 合成是一种剖析矩阵数据中的结构特征、提取次要数据特色量的办法。
EOF 合成的剖析对象，就是咱们在气象科研和业务工作中，常常要剖析各种气象要素场，如海表温度场、气压场、降水场等。这些气象要素场大多由不规则的网格点所组成。如果咱们抽取这些场的某一段历史期间的材料，就形成一组以网格点为空间点 (可看作多个变量) 的随工夫变动的样本。
EOF 合成所做的，就是把随工夫变动的气象要素场合成为 空间函数局部 和 工夫函数局部，且这两局部互相正交。
其中，空间函数局部可在肯定水平上概括因素场的空间散布特点，这部分是不随工夫变动的；而工夫函数局部则由空间点 (即多个变量) 的线性组合所形成，称为主重量。
因为这些主重量中前几个就可能占有原空间点 (即多个变量) 总方差的很大部分，即通过 EOF 合成能够很容易地将原始因素场的变动信息稀释在前几个主重量及其对应的空间函数上，所以，EOF 合成常被用于气象要素场时空变动特色法则的提取和钻研。
2. 基本原理
在气象科研和业务工作中，咱们常常会剖析某一气象要素的时空数据集，它们大多是 3 维的，包含 2 维的空间场以及 1 维的工夫序列。其中 2 维空间场对应着一个个的空间点，1 维的工夫序列对应着观测的工夫点。
当初假如咱们要对某一个气象要素的时空数据集进行钻研，它共有 m 个空间点 ($x_1,x_2,···,x_m$), n 个工夫点 ($t_1,t_2,···,t_n$)。那对于这样一个时空数据集，咱们就能够用一个 $m×n$ 的矩阵 $F$ 来示意。
其中 $F$ 的行对应着某个空间点在观测期内所有工夫点上的值，$F$ 的列对应着某个工夫点在观测区域内所有空间点上的值。
$$
F=
\begin{pmatrix}
f_{11} & f_{12} & ··· & f_{1n}\\
f_{21} & f_{22} & ··· & f_{2n}\\
··· & ··· & ··· & ···\\
f_{m1} & f_{m2} & ··· & f_{mn}\\
\end{pmatrix}
$$
而咱们的 EOF 合成办法，就是将时空数据集 $F$ 合成为空间函数 $V$ 和工夫函数 $Z$ 两局部，即
$$
F=VZ
$$
$$
\begin{pmatrix}
f_{11} & f_{12} & ··· & f_{1n}\\
f_{21} & f_{22} & ··· & f_{2n}\\
··· & ··· & ··· & ···\\
f_{m1} & f_{m2} & ··· & f_{mn}\\
\end{pmatrix}
=
\begin{pmatrix}
v_{11} & v_{12} & ··· & v_{1m}\\
v_{21} & v_{22} & ··· & v_{2m}\\
··· & ··· & ··· & ···\\
v_{m1} & v_{m2} & ··· & v_{mm}\\
\end{pmatrix}
\begin{pmatrix}
z_{11} & z_{12} & ··· & z_{1n}\\
z_{21} & z_{22} & ··· & z_{2n}\\
··· & ··· & ··· & ···\\
z_{m1} & z_{m2} & ··· & z_{mn}\\
\end{pmatrix}
$$
也可示意为
$$
F_{it}=\sum_{k=1}^{m}v_{ik}z_{kt}=v_{i1}z_{1t}+v_{i2}z_{2t}+···+v_{im}z_{mt}
$$
其中下标 $i=1,2,···,m$ 示意空间点，下标 $t=1,2,···,n$ 示意工夫点。
因而，该合成式的含意为：气象要素场中第 $i$ 个空间点上的第 $t$ 次观测值，能够看作是 m 个空间函数 $v_{i1},v_{i2},···,v_{im}$ 和 m 个工夫函数 $z_{1t},z_{2t},···,z_{mt}$ 的线性组合。
在空间函数 $V$ 中，它的每一列(如第 k 列，可示意为 $\mathbf{v_k}=(v_{1k},v_{2k},···,v_{mk})$ 都能够看作为一个典型场，称为一个空间模态(EOF)，它只跟空间点无关，而与工夫无关。
在工夫函数 $Z$ 中，它的每一行(如第 k 行，可示意为 $\mathbf{z_k}=(z_{k1},z_{k2},···,z_{kn})$ 都能够看作是一个工夫序列，称为工夫系数(PC)。
在此基础上，时空数据集 $F$ 中第 t 时刻的空间场(共有 m 个空间点)，就能够示意为
$$
\begin{pmatrix}
f_{1t} \\
f_{2t} \\
··· \\
f_{mt} \\
\end{pmatrix}
=
\begin{pmatrix}
v_{11} \\
v_{21} \\
··· \\
v_{m1} \\
\end{pmatrix}
z_{1t}
+
\begin{pmatrix}
v_{12} \\
v_{22} \\
··· \\
v_{m2} \\
\end{pmatrix}
z_{2t}
+
···
+
\begin{pmatrix}
v_{1m} \\
v_{2m} \\
··· \\
v_{mm} \\
\end{pmatrix}
z_{mt}
$$
该式表明，时空数据集中第 t 时刻的空间场能够示意为 m 个空间模态($\mathbf{v_1},\mathbf{v_2},···,\mathbf{v_m}$)，依照不同的权重线性 ($z_{1t},z_{2t},···,z_{mt}$) 叠加而成。
对于任意一个空间模态 $\mathbf{v_k}$ 而言，其在不同时刻对应的权重有所不同，若将其在全副观测工夫点所对应的权重 $(z_{k1},z_{k2},···,z_{kn})$ 列在一起，就是咱们下面提到的工夫系数 $\mathbf{z_k}$。
每一个空间模态都对应着一个本人的工夫系数 $\mathbf{z_k}$，如果 $\mathbf{z_k}$ 中的正值的绝对值越大，表明该工夫点所对应的数据集与该工夫系数所对应的空间模态类似度越高；如果 $\mathbf{z_k}$ 中的负值的绝对值越大，则表明该工夫点所对应的数据集越偏向于和该工夫系数所对应的空间模态呈相同的状态。
除此以外，咱们还能够依据每一个空间模态 (或工夫系数) 计算该空间模态 (或工夫系数) 对应的 方差贡献率 ，方差贡献率越高，表明该空间模态(或工夫系数) 所蕴含的对于原始数据的变动信息越多。
咱们通常依据方差贡献率的大小对空间模态进行排序，方差贡献率最大的称为第一空间模态，方差贡献率第二大的称为第二空间模态，以此类推。
因为所有空间模态 (或工夫系数) 所对应的方差贡献率加起来等于 1，所以，当方差贡献率最高的几个空间模态 (或工夫系数) 加起来失去的累计方差贡献率较大时(比方达到 80%)，咱们就会用它们几个来示意原始数据中所蕴含的次要特色。
但须要强调的是，EOF 合成只是一种统计办法，其失去的合成后果自身没有任何的物理意义，须要咱们结合实际及专业知识去对后果进行了解和解释的。
3. 在 Python 实现 EOF 合成
import struct
import numpy as np
import xarray as xr
from eofs.xarray import Eof
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import matplotlib.ticker as ticker
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import (LongitudeFormatter, LatitudeFormatter, LatitudeLocator, LongitudeLocator)
import cmaps
import warnings
warnings.filterwarnings('ignore')3.1 进行 EOF 剖析
# 读取数据
data = np.fromfile('ssta.dat', dtype='float32')
data = np.reshape(data,(56, 11, 336), order='F')
lon = np.arange(160, 272, 2) #经度 160E~90W（距离 2 度)
lat = np.arange(-10, 12, 2) #纬度 10S~10N（距离 2 度）time = xr.cftime_range(start='1979-01', end='2006-12', freq='MS') #工夫 1979.1~2006.12（逐月)
ds = xr.Dataset({'ssta': ([ 'lon', 'lat', 'time'], data)},
coords={'lon': lon, 'lat': lat, 'time': time})
#测验数据读取
fig = plt.figure(figsize=(18, 12), dpi=300)
projection = ccrs.PlateCarree(central_longitude=180) #指定投影为经纬度投影，核心经纬度为 180°
ax1 = fig.add_subplot(1, 1, 1, projection=projection) #绘制第一个子图
projection = ccrs.PlateCarree(central_longitude=180) #指定投影为经纬度投影，并指定核心经度为 180°
# 设置地图范畴，经度为（160, 270)，纬度为（-10, 10）ax1.set_extent([-20, 90, -10, 10], crs=ccrs.PlateCarree(central_longitude=180))
# 设置经纬度标签
ax1.set_xticks([-20, -10, 0, 10, 20, 30, 40, 50, 60, 70, 80, 90], crs=projection)
ax1.set_yticks([-10, -5, 0, 5, 10], crs=projection)
lon_formatter = LongitudeFormatter(zero_direction_label=True) #给标签增加对应的 N,S,E,W
lat_formatter = LatitudeFormatter()
ax1.xaxis.set_major_formatter(lon_formatter)
ax1.yaxis.set_major_formatter(lat_formatter)
ax1.contourf(ds.lon,ds.lat,ds.ssta[:,:,0].transpose('lat', 'lon'), np.linspace(-1,1,51),cmap=cmaps.BlueWhiteOrangeRed,transform=ccrs.PlateCarree()) #在第一个地图上绘制第一空间模态的填充色图，EOFs[0]示意从 EOF 中取出第一个空间模态
ax1.coastlines(color='k', lw=0.5) #增加海岸线
ax1.add_feature(cfeature.LAND, facecolor='white') #增加海洋# 转置数据
ds = ds.transpose('time', 'lat', 'lon')
#计算网格点的权重
coslat = np.cos(np.deg2rad(ds.coords['lat'].values))
wgts = np.sqrt(coslat)[..., np.newaxis]
# 计算 EOF 合成
solver = Eof(ds['ssta'], weights=wgts, center=False)
eofs = solver.eofs(neofs=10, eofscaling=2)
pcs = solver.pcs(npcs=10, pcscaling=1)
variance = solver.varianceFraction(neigs=10)
# 绘制方差解释率图
plt.figure(figsize=(8, 6),dpi=200)
plt.plot(np.arange(1, 11), variance[:10]*100, 'o-')
plt.xlabel('EOF modes')
plt.ylabel('Variance explained (%)')
plt.title('Variance explained by EOF modes')
plt.grid()
# 在每个点上增加文本标签
for i, v in enumerate(variance.values[:10]*100):
plt.text(i+1, v+0.9, f'{v:.2f}', fontsize=8, color='black')
plt.show()# 计算特征值和 NorthTest 的后果
eigenvalues = solver.eigenvalues(neigs=10)
errors = solver.northTest(neigs=10)
# 绘制 Error Bar 图
plt.figure(figsize=(8, 6),dpi=200)
plt.errorbar(np.arange(1, 11), eigenvalues, yerr=errors, fmt='o-', capsize=5)
plt.xlabel('EOF modes')
plt.ylabel('Eigenvalues')
plt.title('NorthTest(with error bars)')
plt.grid()
plt.show()3.2 给出前 4 个模态的空间散布和工夫序列（前 4 个 modes 的误差棒无相重叠局部）
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LongitudeFormatter, LatitudeFormatter
import cmaps
def mapart(ax):
'''增加地图元素'''
# 指定投影为经纬度投影，并指定核心经度为 180°
projection = ccrs.PlateCarree(central_longitude=180)
# 设置地图范畴，经度为（160, 270)，纬度为（-10, 10）ax.set_extent([-20, 90, -10, 10], crs=ccrs.PlateCarree(central_longitude=180))
# 设置经纬度标签
ax.set_xticks([-20, -10, 0, 10, 20, 30, 40, 50, 60, 70, 80, 90], crs=projection)
ax.set_yticks([-10, -5, 0, 5, 10], crs=projection)
# 给标签增加对应的 N,S,E,W
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
# 增加海岸线
#ax.coastlines(color='k', lw=0.5)
# 增加海洋
#ax.add_feature(cfeature.LAND, facecolor='white')
def eof_contourf(EOFs, PCs, pers):
'''绘制 EOF 填充图'''
plt.close
# 将方差转换为百分数的模式，如 0.55 变为 55%
pers=(pers*100).values
# 指定画布大小以及像素
fig = plt.figure(figsize=(18, 12), dpi=300)
# 指定投影为经纬度投影，核心经纬度为 180°
projection = ccrs.PlateCarree(central_longitude=180)
# 将横坐标转换为日期格局
dates_num = [mdates.date2num(t) for t in time]
# 绘制第一个子图
ax1 = fig.add_subplot(4, 2, 1, projection=projection)
# 为第一个子图增加地图底图
mapart(ax1)
# 在第一个地图上绘制第一空间模态的填充色图，EOFs[0]示意从 EOF 中取出第一个空间模态
p = ax1.contourf(ds.lon,
ds.lat,
EOFs[0],
np.linspace(-1, 1, 51),
cmap=cmaps.BlueWhiteOrangeRed,
transform=ccrs.PlateCarree(),
extend='both')
# 为第一个子图增加题目，其中，pers[0]示意从 pers 中取出第一个空间模态对应的方差贡献率
ax1.set_title('mode1 (%s' % (round(pers[0], 2))+"%)", loc='left')
ax1.set_aspect(1.6)
# 绘制第二个子图
ax2 = fig.add_subplot(4, 2, 2)
ax2.set_xticks(dates_num[::108]) # 每 9 年的 1 月作为横坐标标签
ax2.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))
ax2.xaxis.set_major_locator(ticker.FixedLocator(dates_num[::108])) # 每 9 年的 1 月作为横坐标刻度线
b = ax2.bar(dates_num, PCs[:, 0], width=50, color='r')
# 对工夫系数值小于 0 的柱子设置为蓝色
for bar, height in zip(b, PCs[:, 0]):
if height < 0:
bar.set(color='blue')
# 为第二个子图增加题目
ax2.set_title('PC1' % (round(pers[0], 2)), loc='left')
# 前面 ax3 和 ax5 与 ax1 绘制办法雷同，ax4 和 ax6 与 ax2 绘制办法雷同，不再正文
ax3 = fig.add_subplot(4, 2, 3, projection=projection)
mapart(ax3)
pp = ax3.contourf(ds.lon,
ds.lat,
EOFs[1],
np.linspace(-1, 1, 51),
cmap=cmaps.BlueWhiteOrangeRed,
transform=ccrs.PlateCarree(),
extend='both')
ax3.set_title('mode2 (%s' % (round(pers[1], 2))+"%)", loc='left')
ax3.set_aspect(1.6)
ax4 = fig.add_subplot(4, 2, 4)
ax4.set_xticks(dates_num[::108]) # 每 9 年的 1 月作为横坐标标签
ax4.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))
ax4.xaxis.set_major_locator(ticker.FixedLocator(dates_num[::108])) # 每 9 年的 1 月作为横坐标刻度线
ax4.set_title('PC2' % (round(pers[1], 2)), loc='left')
bb = ax4.bar(dates_num, PCs[:, 1], width=50, color='r')
for bar, height in zip(bb, PCs[:, 1]):
if height < 0:
bar.set(color='blue')
ax5 = fig.add_subplot(4, 2, 5, projection=projection)
mapart(ax5)
ppp = ax5.contourf(ds.lon,
ds.lat,
EOFs[2],
np.linspace(-1, 1, 51),
cmap=cmaps.BlueWhiteOrangeRed,
transform=ccrs.PlateCarree(),
extend='both')
ax5.set_title('mode3 (%s' % (round(pers[2], 2))+"%)", loc='left')
ax5.set_aspect(1.6)
ax6 = fig.add_subplot(4, 2, 6)
ax6.set_xticks(dates_num[::108]) # 每 9 年的 1 月作为横坐标标签
ax6.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))
ax6.xaxis.set_major_locator(ticker.FixedLocator(dates_num[::108])) # 每 9 年的 1 月作为横坐标刻度线
ax6.set_title('PC3' % (round(pers[2], 2)), loc='left')
bbb = ax6.bar(dates_num, PCs[:, 2], width=50, color='r')
for bar, height in zip(bbb, PCs[:, 2]):
if height < 0:
bar.set(color='blue')
ax7 = fig.add_subplot(4, 2, 7, projection=projection)
mapart(ax7)
pppp = ax7.contourf(ds.lon,
ds.lat,
EOFs[3],
np.linspace(-1, 1, 51),
cmap=cmaps.BlueWhiteOrangeRed,
transform=ccrs.PlateCarree(),
extend='both')
ax7.set_title('mode4 (%s' % (round(pers[3], 2))+"%)", loc='left')
ax7.set_aspect(1.6)
ax8 = fig.add_subplot(4, 2, 8)
ax8.set_xticks(dates_num[::108]) # 每 9 年的 1 月作为横坐标标签
ax8.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))
ax8.xaxis.set_major_locator(ticker.FixedLocator(dates_num[::108])) # 每 9 年的 1 月作为横坐标刻度线
#ax8.set_xlim(dates_num[0], dates_num[-1])
bbbb = ax8.bar(dates_num, PCs[:, 3], width=50, color='r')
for bar, height in zip(bbbb, PCs[:, 3]):
if height < 0:
bar.set(color='blue')
ax8.set_title('PC4' % (round(pers[3], 2)), loc='left')
# 为柱状图增加 0 标准线
ax2.axhline(y=0, linewidth=1, color='k', linestyle='-')
ax4.axhline(y=0, linewidth=1, color='k', linestyle='-')
ax6.axhline(y=0, linewidth=1, color='k', linestyle='-')
ax8.axhline(y=0, linewidth=1, color='k', linestyle='-')
# 在图下边留白边放 colorbar
fig.subplots_adjust(bottom=0.1)
# colorbar 地位：左 下 宽 高
l = 0.25
b = 0.04
w = 0.6
h = 0.015
# 对应 l,b,w,h；设置 colorbar 地位；rect = [l, b, w, h]
cbar_ax = fig.add_axes(rect)
# 绘制 colorbar
c = plt.colorbar(
p,
cax=cbar_ax,
orientation='horizontal',
aspect=20,
pad=0.1)
# 设置 colorbar 的标签大小
c.ax.tick_params(labelsize=14)
# 调整子图之间的间距
plt.subplots_adjust(wspace=0.1, hspace=0.3)
plt.show()
eof_contourf(eofs[:4]*(-1),pcs[:,:4]*(-1),variance[:4])参考：
教训正交方程 EOF- 原理，代码，气象学实例，文献解读
《气象统计分析与预报办法第 4 版 黄嘉佑》
USTC 气象统计办法课程 - 郑建秋
教训正交合成(EOF) 夏子涵 夏子涵 ”)
本文由 mdnice 多平台公布