eof是什么意思,eof实际应用
原理
经验正交函数分析方法(empirical orthogonal function,缩写为EOF),也称特征向量分析(eigenvector analysis),或者主成分分析(principal component annlysis,缩写为PCA),是一种分析矩阵数据中的结构特征,提取主要数据特征量的一种方法。特征向量对应的是空间样本,也称空间特征向量或者空间模态(EOF),在一定程度上反映要素场的空间分布特点;主成分对应的是时间变化,也称时间系数(PC),反映相应空间模态随时间的权重变化。因此称EOF分析为时空分解,即X=EOFm×m×PCm×n,m为站点数,n为年数。算法如下:
(1)选定要分析的数据,进行数据预处理,通常处理为距平的形式,得到一个数据矩阵Xm×n
(2)计算X与其转置矩阵XT的交叉积,得到方阵
(3)计算方阵C的特征根(λ1,…,λm)和特征向量Vm×m,二者满足
其中E是m×m维对角阵,即
一般将特征根λ按从大到小顺序排列,即λ1>λ2>…>λm。因为数据X是真实的观测值,所以对应的λ大于或等于0。每个非0的特征根对应一列特征向量值,也称EOF。如λ1对应的特征向量值称第一个EOF模态,也就是V的第一列。
(4)计算主成分。将EOF投影到原始资料矩阵上,就得到所有空间特征向量对应的时间系数(即主成分),即
其中PC中每行数据就是对应每个特征向量的时间系数。
(5)计算贡献率。矩阵X的方差大小可以简单的用特征根的大小来表示。λ越高说明其对应的模态越重要,对总方差的贡献越大。
第k个模态对总的方差解释率为
(6)显著性检验。即使是随机数或者虚假数据,放在一起进行EOF分析,也可以将其分解成一系列的空间特征向量和主成分。因此,实际资料分析中得到的空间模态是否是随机的,需要记性统计检验。研究指出,在95%置信度水平下的特征根的误差:

λ是特征根,是数据的有效自由度。将λ按顺序依次检查,标上误差范围。如果前后两个λ之间误差范围有重叠,则没有通过显著性检验。
代码实现
import numpy as npimport cartopy.crs as ccrsimport cartopy.feature as cfeatureimport matplotlib.pyplot as plt from eofs.multivariate.standard import MultivariateEoffrom eofs.standard import Eofdef eof_analys(data, lat): # 计算权重:纬度cos值开方 coslat = np.cos(np.deg2rad(lat)) wgts = np.sqrt(coslat)[..., np.newaxis] # 做EOF分析 solver = Eof(data, weights=wgts) EOFs = solver.eofsAsCorrelation() # 空间模态 PCs = solver.pcs(npcs=3, pcscaling=1) # 时间序列,取前三个模态 # 方差贡献 eigen_Values = solver.eigenvalues() percentContrib = eigen_Values * 100. / np.sum(eigen_Values) # 返回空间模态,时间序列和方差贡献 return EOFs, PCs, percentContrib# 加入一个可视化程序,打包成一个计算和可视化一体的函数def mapart(ax): # 添加地图元素 projection = ccrs.PlateCarree() ax.coastlines(color='k', lw=0.5) ax.add_feature(cfeature.LAND, facecolor='white') # 设置地图范围 ax.set_extent([100, 150, 0, 50], crs=ccrs.PlateCarree()) # 设置经纬度标签 ax.set_xticks([100, 115, 130, 145], crs=projection) ax.set_yticks([0, 10, 20, 30, 40, 50], crs=projection) lon_formatter = LongitudeFormatter(zero_direction_label=True) lat_formatter = LatitudeFormatter() ax.xaxis.set_major_formatter(lon_formatter) ax.yaxis.set_major_formatter(lat_formatter)def eof_contourf(EOFs, PCs, pers, name): plt.close fig = plt.figure(figsize=(12, 12)) projection = ccrs.PlateCarree() year = range(1982, 2020) ax1 = fig.add_subplot(3, 2, 1, projection=projection) mapart(ax1) p = ax1.contourf(lon, lat, EOFs[0], np.linspace(-1, 1, 21), cmap=cmaps.BlueWhiteOrangeRed) ax1.set_title('mode1 (%s' % (round(pers[0], 2)) + "%)", font2, loc='left') ax2 = fig.add_subplot(3, 2, 2) ax2.plot(year, PCs[:, 0], color='k', linewidth=1.2, linestyle='--') y1 = np.poly1d(np.polyfit(year, PCs[:, 0], 1)) ax2.plot(year, y1(year), 'k--', linewidth=1.2) b = ax2.bar(year, PCs[:, 0], color='r') # 对y值大于0设置为蓝色 小于0的柱设置为绿色 for bar, height in zip(b, PCs[:, 0]): if height < 0: bar.set(color='blue') ax2.set_title('PC1' % (round(pers[0], 2)), font2, loc='left') ax3 = fig.add_subplot(3, 2, 3, projection=projection) mapart(ax3) pp = ax3.contourf(lon, lat, EOFs[1], np.linspace(-1, 1, 21), cmap=cmaps.BlueWhiteOrangeRed) ax3.set_title('mode2 (%s' % (round(per2, 2)) + "%)", font2, loc='left') ax4 = fig.add_subplot(3, 2, 4) ax4.plot(year, PCs[:, 1], color='k', linewidth=1.2, linestyle='--') ax4.set_title('PC2' % (round(per2, 2)), font2, loc='left') print(np.polyfit(year, PCs[:, 1], 1)) y2 = np.poly1d(np.polyfit(year, PCs[:, 1], 1)) # print(y2) ax4.plot(year, y2(year), 'k--', linewidth=1.2) bb = ax4.bar(year, PCs[:, 1], color='r') # 对y值大于0设置为蓝色 小于0的柱设置为绿色 for bar, height in zip(bb, PCs[:, 1]): if height < 0: bar.set(color='blue') ax5 = fig.add_subplot(3, 2, 5, projection=projection) mapart(ax5) ppp = ax5.contourf(lon, lat, EOFs[2], np.linspace(-1, 1, 21), cmap=cmaps.BlueWhiteOrangeRed) ax5.set_title('mode3 (%s' % (round(pers[2], 2)) + "%)", font2, loc='left') ax6 = fig.add_subplot(3, 2, 6) ax6.plot(year, PCs[:, 2], color='k', linewidth=1.2, linestyle='--') ax6.set_title('PC3' % (round(per3, 2)), font2, loc='left') y3 = np.poly1d(np.polyfit(year, PCs[:, 2], 1)) # print(y3) ax6.plot(year, y3(year), 'k--', linewidth=1.2) bbb = ax6.bar(year, pers[:, 2], color='r') # 对y值大于0设置为蓝色 小于0的柱设置为绿色 for bar, height in zip(bbb, PCs[:, 2]): if height < 0: bar.set(color='blue') # 添加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='-') # 在图下边留白边放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) c = plt.colorbar(pp, cax=cbar_ax, orientation='horizontal', aspect=20, pad=0.1) c.ax.tick_params(labelsize=14) plt.subplots_adjust(wspace=-0.12, hspace=0.3) plt.savefig('eof_%s.jpg' % name, dpi=300, format='jpg', bbox_inches='tight', transparent=True, pad_inches=0) plt.show()def eof_analyze(data, lat, name): EOFs, PCs, per = eof_anyls(data, lat) print('前三个模态的方差贡献分别是:%s,%s,%s' % (round(per1, 2), round(per2, 2), round(per3, 2))) eof_contourf(EOFs, PCs, percentContrib, name)# 调用上面写的eof分析和可视化函数# 输入一个三维数据,纬度,和要保存的图片名字eof_analyze(data, lat, 'fig_name')
结果呈现

小编说
在本篇文章中我们学习了如何在Python中利用EOF分析进行实际运用。在实际运行的时候发现,由于pip环境出了问题导致导入Cartopy库的时候一直导入不成功,最后重新安装了pip环境,又对代码进行了再一次运行。
文中分享的方法作者自定义了eof分析函数和画图函数,然后直接调用所写的这两个函数做分析,这样写代码可以简洁一些而且方便使用。把作者分享的所有代码都复制粘贴到一个文件里,然后最后一行eof_analyze(data,lat,“”)就是在调用前面自定义那个函数,大家只要在这里输入所需要的数据就可以了。当然还需要在mapart这个函数里更改你的地图范围。
本文地址:百科问答频道 https://www.neebe.cn/wenda/930920.html,易企推百科一个免费的知识分享平台,本站部分文章来网络分享,本着互联网分享的精神,如有涉及到您的权益,请联系我们删除,谢谢!