使用python对1979-2004年的SSTA海温异常现象进行EOF分析并绘图

学习目标：

用python进行EOF分析

例子

EOF分析是气象分析中常见的一种分析方法，也被称为经验正交函数。经过EOF分析，可以将几十年的海温数据变成几个空间模态和时间序列，这样就可以通过空间模态大致分析一些变化趋势。

选取的是1979—2004年的海温数据，下载的网站是
https://www.metoffice.gov.uk/hadobs/hadisst/index.html
具体参考这篇文章https://blog.csdn.net/shayuxing/article/details/122291093
输出结果

相应的代码

#导入模块
import numpy as np #numpy一种开源计算库，提供和处理N维数组对象
import cartopy.crs as ccrs #导入地图绘图包cartopy
import cartopy.feature as cfeature
import matplotlib.pyplot as plt #导入2D绘图库matplotlib
from eofs.standard import Eof #导入正交分解函数EOF
import xarray as xr #导入处理多维数组数据集库xarray
from cartopy.mpl import ticker
 
#读取数据
path="C:\\Users\\huain\\Desktop\\HadISST_sst.nc"
SST=xr.open_dataset(path).sel(latitude=slice(30,-30),longitude=slice(-180,180),time=slice("1979","2004")) #选取数据范围和时间序列
sst1=SST.sst[:]
sst2=np.array(sst1)
lat=SST.latitude[:]
lon=SST.longitude[:] 
sst=np.array(sst1)
ano=sst1.groupby('time.month')-sst1.groupby('time.month').mean('time', skipna=True)
ano1=np.array(ano)
 
#计算纬度权重
lat=np.array(lat)
coslat=np.cos(np.deg2rad(lat))
wgts = np.sqrt(coslat)[..., np.newaxis]

#创建EOF分解器
solver=Eof(ano1,weights=wgts)
eof=solver.eofsAsCorrelation(neofs=4)
pc=solver.pcs(npcs=4,pcscaling=1)
var=solver.varianceFraction(neigs=4)

#绘图 
fig=plt.figure(figsize=(15,15))#创建画布
proj=ccrs.PlateCarree(central_longitude=180)
leftlon,rightlon,lowerlat,upperlat=(-180,180,-30,30)#绘制区域
lon_formatter=ticker.LongitudeFormatter()
lat_formatter=ticker.LatitudeFormatter()

#绘制第一模态
fig_ax1=fig.add_axes([0.05,0.75,0.4,0.3],projection=proj)#调整位置大小
fig_ax1.set_extent([leftlon,rightlon,lowerlat,upperlat],crs=ccrs.PlateCarree())
fig_ax1.add_feature(cfeature.OCEAN,edgecolor='black')
fig_ax1.add_feature(cfeature.LAKES,alpha=0.5)
fig_ax1.add_feature(cfeature.COASTLINE,lw=1)
fig_ax1.set_xticks(np.arange(leftlon,rightlon,40),crs=ccrs.PlateCarree())
fig_ax1.set_yticks(np.arange(lowerlat,upperlat+5,10),crs=ccrs.PlateCarree())
fig_ax1.xaxis.set_major_formatter(lon_formatter)
fig_ax1.yaxis.set_major_formatter(lat_formatter)
rivers_110m = cfeature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', '110m')
fig_ax1.set_title('(a) EOF1(HadISSTA from 1979-2004)',loc='left',fontsize =10)
fig_ax1.set_title( '%.2f%%' % (var[0]*100),loc='right',fontsize =10)
c1=fig_ax1.contourf(lon,lat, eof[0,:,:], levels=np.arange(-0.9,1.0,0.1), zorder=0, extend = 'both',transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)

#绘制第二模态 
fig_ax2=fig.add_axes([0.05,0.55,0.4,0.3],projection=proj)
fig_ax2.set_extent([leftlon,rightlon,lowerlat,upperlat],crs=ccrs.PlateCarree())
fig_ax2.add_feature(cfeature.OCEAN,edgecolor='black')
fig_ax2.add_feature(cfeature.LAKES,alpha=0.5)
fig_ax2.add_feature(cfeature.COASTLINE,lw=1)
fig_ax2.set_xticks(np.arange(leftlon,rightlon,40),crs=ccrs.PlateCarree())
fig_ax2.set_yticks(np.arange(lowerlat,upperlat+5,10),crs=ccrs.PlateCarree())
fig_ax2.xaxis.set_major_formatter(lon_formatter)
fig_ax2.yaxis.set_major_formatter(lat_formatter)
rivers_110m = cfeature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', '110m')
fig_ax2.set_title('(c) EOF2(HadISSTA from 1979-2004)',loc='left',fontsize =10)
fig_ax2.set_title( '%.2f%%' % (var[1]*100),loc='right',fontsize =10)
c2=fig_ax2.contourf(lon,lat, eof[1,:,:], levels=np.arange(-0.9,1.0,0.1), zorder=0, extend = 'both',transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)

#绘制第三模态 
fig_ax3=fig.add_axes([0.05,0.35,0.4,0.3],projection=proj)
fig_ax3.set_extent([leftlon,rightlon,lowerlat,upperlat],crs=ccrs.PlateCarree())
fig_ax3.add_feature(cfeature.OCEAN,edgecolor='black')
fig_ax3.add_feature(cfeature.LAKES,alpha=0.5)
fig_ax3.add_feature(cfeature.COASTLINE,lw=1)
fig_ax3.set_xticks(np.arange(leftlon,rightlon,40),crs=ccrs.PlateCarree())
fig_ax3.set_yticks(np.arange(lowerlat,upperlat+5,10),crs=ccrs.PlateCarree())
fig_ax3.xaxis.set_major_formatter(lon_formatter)
fig_ax3.yaxis.set_major_formatter(lat_formatter)
rivers_110m = cfeature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', '110m')
fig_ax3.set_title('(e) EOF3(HadISSTA from 1979-2004)',loc='left',fontsize =10)
fig_ax3.set_title( '%.2f%%' % (var[2]*100),loc='right',fontsize =10)
c3=fig_ax3.contourf(lon,lat, eof[2,:,:], levels=np.arange(-0.9,1.0,0.1), zorder=0, extend = 'both', transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)

#绘制第四模态 
fig_ax4=fig.add_axes([0.05,0.15,0.4,0.3],projection=proj)
fig_ax4.set_extent([leftlon,rightlon,lowerlat,upperlat],crs=ccrs.PlateCarree())
fig_ax4.add_feature(cfeature.OCEAN,edgecolor='black')
fig_ax4.add_feature(cfeature.LAKES,alpha=0.5)
fig_ax4.add_feature(cfeature.COASTLINE,lw=1)
fig_ax4.set_xticks(np.arange(leftlon,rightlon,40),crs=ccrs.PlateCarree())
fig_ax4.set_yticks(np.arange(lowerlat,upperlat+5,10),crs=ccrs.PlateCarree())
fig_ax4.xaxis.set_major_formatter(lon_formatter)
fig_ax4.yaxis.set_major_formatter(lat_formatter)
rivers_110m = cfeature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', '110m')
fig_ax4.set_title('(g) EOF4(HadISSTA from 1979-2004)',loc='left',fontsize =10)
fig_ax4.set_title( '%.2f%%' % (var[3]*100),loc='right',fontsize =10)
c4=fig_ax4.contourf(lon,lat, eof[3,:,:], levels=np.arange(-0.9,1.0,0.1), zorder=0, transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)


#绘制colorbar 
cbposition=fig.add_axes([0.25, 0.13, 0.5, 0.015])#位置
fig.colorbar(c1,cax=cbposition,orientation='horizontal',format='%.1f')


#绘制第一模态时间序列
ax=fig.add_axes([0.5,0.82,0.35,0.15])
bar=ax.bar(np.arange(1979,2005,1/12),height=pc[:,0],color='blue',align="center",width=0.1,linewidth=0.1,bottom=None,edgecolor='black') #PC1条形
for bar, height in zip(bar, pc[:, 0]):
        if height > 0:
            bar.set(color='red')
ax.set_title('(b) PC1',loc='left',fontsize = 10)
           
#绘制第二模态时间序列 
ax=fig.add_axes([0.5, 0.61, 0.35, 0.15])
bar=ax.bar(np.arange(1979,2005,1/12),height=pc[:,1],color='blue',align="center",width=0.1,linewidth=0.1,bottom=None,edgecolor='black') #PC2条形
for bar, height in zip(bar, pc[:, 1]):
        if height > 0:
            bar.set(color='red')
ax.set_title('(d) PC2',loc='left',fontsize = 10)

#绘制第三模态时间序列 
ax=fig.add_axes([0.5, 0.40, 0.35, 0.15])
bar=ax.bar(np.arange(1979,2005,1/12),height=pc[:,2],color='blue',align="center",width=0.1,linewidth=0.1,bottom=None,edgecolor='black') #PC3条形
for bar, height in zip(bar, pc[:, 2]):
        if height > 0:
            bar.set(color='red')
ax.set_title('(f) PC3',loc='left',fontsize = 10)

#绘制第四模态时间序列 
ax=fig.add_axes([0.5, 0.19, 0.35, 0.15])
bar=ax.bar(np.arange(1979,2005,1/12),height=pc[:,3],color='blue',align="center",width=0.1,linewidth=0.1,bottom=None,edgecolor='black') #PC4条形
for bar, height in zip(bar, pc[:, 3]):
        if height > 0:
            bar.set(color='red')
ax.set_title('(h) PC4',loc='left',fontsize = 10)
plt.show()
 
fig=plt.figure(figsize=(10,6))
ax=fig.add_axes([0.2,0.2,0.5,0.5])

#把PC1、PC2、PC3时间序列绘制在一张图上 
ax.plot(np.arange(1979,2005,1/12),pc[:,0],linewidth=1,linestyle='-',color='r',label='PC1') #PC1线

bar=ax.bar(np.arange(1979,2005,1/12),height=pc[:,1],color='blue',align="center",width=0.1,linewidth=0.1,bottom=None,edgecolor='black',label='PC2') #PC2条形

ax.plot(np.arange(1979,2005,1/12),pc[:,2],linestyle='--',linewidth=1,color='black',label='PC3')#PC3线
ax.set_ylim(-4,4)
ax.set_title("PC")
ax.set_xlabel("Time")
ax.set_ylabel("y")
plt.legend()
plt.grid() 
plt.show()