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

学习目标：

用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() 








1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161

1