EOF(经验正交分解)是气候研究中常用的研究变量时空变化特征的分析方法，短期气候课中都学过中国东部夏季降水通过EOF分解可以分为三类雨型，在NCL中，EOF有直接的函数可以调用，而在python中，也有对应的库可以直接使用。
完整的说明见：
https://github.com/ajdawson/eofs
通过

conda install -c conda-forge eofs

可以直接安装。
由于数据以及EOF所选范围原因，可能结果存在一些差异，但是问题不大233333
先上图吧：

图的分析就不献丑了。主要还是分享一下绘制方法。
首先是EOF的计算

print(summer_mean_tmp.shape)
#(55, 82, 142)
print(pre_lat.shape)
#(82,)
print(pre_lon.shape)
#(142,)

这是用来分解的数据，55年夏季降水，纬度和经度。接下来进行EOF分解：

from eofs.standard import Eoflat = np.array(pre_lat)coslat = np.cos(np.deg2rad(lat))wgts = np.sqrt(coslat)[..., np.newaxis]#计算纬度权重solver = Eof(summer_mean_tmp, weights=wgts)#创建EOF函数eof = solver.eofsAsCorrelation(neofs=3)#获取前三个模态pc = solver.pcs(npcs=3, pcscaling=1)var = solver.varianceFraction()#获取对应的PC序列和解释方差

这样我们获得了：

print(eof.shape)#(3, 82, 142)print(pc.shape)#(55, 3)print(var.shape)#(55, )

至于方差为啥是(55)，我也不记得了，但是var[n]对应第n模态的方差，结果与NCL计算的结果基本一致。
得到了以上的结果，就可以用于绘图了。
首先是模态的绘制，即带地图底图的填色图
我推荐使用matplotlib + cartopy实现

conda install -c conda-forge cartopy
conda install matplotlib

完成安装后

import matplotlib as pltimport cartopy.crs as ccrsimport cartopy.feature as cfeaturefrom cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTERimport cartopy.mpl.ticker as ctickerimport cartopy.io.shapereader as shpreader

这次引入的东西有点多，分别用于绘图，添加地图投影，添加地图特征(海岸线，湖泊等等)，设置地理坐标(XX°E)，设置正确的中国国境线(默认的国境线把一些冲突地区偷吃了)

fig2 = plt.figure(figsize=(15,15))proj = ccrs.PlateCarree(central_longitude=115) #设置一个圆柱投影坐标，中心经度115°Eleftlon, rightlon, lowerlat, upperlat = (70,140,15,55)#设置地图边界范围f2_ax1 = fig2.add_axes([0.1, 0.8, 0.5, 0.3],projection = proj)f2_ax1.set_extent([leftlon, rightlon, lowerlat, upperlat], crs=ccrs.PlateCarree())f2_ax1.add_feature(cfeature.COASTLINE.with_scale('50m'))f2_ax1.add_feature(cfeature.LAKES, alpha=0.5)f2_ax1.set_xticks(np.arange(leftlon,rightlon+10,10), crs=ccrs.PlateCarree())f2_ax1.set_yticks(np.arange(lowerlat,upperlat+10,10), crs=ccrs.PlateCarree())lon_formatter = cticker.LongitudeFormatter()lat_formatter = cticker.LatitudeFormatter()f2_ax1.xaxis.set_major_formatter(lon_formatter)f2_ax1.yaxis.set_major_formatter(lat_formatter)f2_ax1.set_title('(a)',loc='left',fontsize =15)f2_ax1.set_title( '%.2f%%' % (var[0]*100),loc='right',fontsize =15)f2_ax1.contourf(pre_lon,pre_lat, eof[0,:,:], levels=np.arange(-0.9,1.0,0.1), zorder=0, extend = 'both',transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)china = shpreader.Reader('bou2_4l.dbf').geometries()f2_ax1.add_geometries(china, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)

通过函数的字面意思应该可以理解大部分语句的意思，'bou2_4l.dbf'是中国国界线shp文件，气象家园有下载。通过类似的方法，缩小子图，覆盖在大图上便可实现南海的绘制。九段线同样有对应的shp文件。
EOF的另一部分便是PC序列的绘制，红正蓝负的实现通过如下循环换实现：

color1=[]for i in range(1961,2016):
    if pc[i-1961,0] >=0:
        color1.append('red')
    elif pc[i-1961,0] <0:
        color1.append('blue')

或许有更简便的方法，目前我还没想到2333

f2_ax4 = fig2.add_axes([0.65, 0.8, 0.5, 0.3])
f2_ax4.set_title('(d)',loc='left')
f2_ax4.set_ylim(-2.5,2.5)
f2_ax4.axhline(0,linestyle="--")
f2_ax4.bar(years,pc[:,0],color=color1)
#f2_ax4.plot(years,pc[:,0])
#只需将bar换为plot即为折线图

作者：摸鱼咯
链接：https://www.jianshu.com/p/d062e779f2dd
来源：简书
著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。