更新完整版，修复了几个小问题，提供了测试数据下载

最新版见此

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 Eof

lat = 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 plt

import cartopy.crs as ccrs

import cartopy.feature as cfeature

from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

import cartopy.mpl.ticker as cticker

import cartopy.io.shapereader as shpreader

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

fig2 = plt.figure(figsize=(15,15))

proj = ccrs.PlateCarree(central_longitude=115)

#设置一个圆柱投影坐标，中心经度115°E

leftlon, 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即为折线图