更新完整版，修復(fù)了幾個小問題，提供了測試數(shù)據(jù)下載
最新版見此

EOF(經(jīng)驗正交分解)是氣候研究中常用的研究變量時空變化特征的分析方法，短期氣候課中都學(xué)過中國東部夏季降水通過EOF分解可以分為三類雨型，在NCL中，EOF有直接的函數(shù)可以調(diào)用，而在python中，也有對應(yīng)的庫可以直接使用。
完整的說明見：
https://github.com/ajdawson/eofs
通過

conda install -c conda-forge eofs

可以直接安裝。
由于數(shù)據(jù)以及EOF所選范圍原因，可能結(jié)果存在一些差異，但是問題不大233333
先上圖吧：

圖的分析就不獻(xiàn)丑了。主要還是分享一下繪制方法。
首先是EOF的計算

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

這是用來分解的數(shù)據(jù)，55年夏季降水，緯度和經(jīng)度。接下來進(jìn)行EOF分解：

from eofs.standard import Eof

lat = np.array(pre_lat)
coslat = np.cos(np.deg2rad(lat))
wgts = np.sqrt(coslat)[..., np.newaxis]
#計算緯度權(quán)重
solver = Eof(summer_mean_tmp, weights=wgts)
#創(chuàng)建EOF函數(shù)
eof = solver.eofsAsCorrelation(neofs=3)
#獲取前三個模態(tài)
pc = solver.pcs(npcs=3, pcscaling=1)
var = solver.varianceFraction()
#獲取對應(yīng)的PC序列和解釋方差

這樣我們獲得了：

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

至于方差為啥是(55)，我也不記得了，但是var[n]對應(yīng)第n模態(tài)的方差，結(jié)果與NCL計算的結(jié)果基本一致。
得到了以上的結(jié)果，就可以用于繪圖了。
首先是模態(tài)的繪制，即帶地圖底圖的填色圖
我推薦使用matplotlib + cartopy實現(xiàn)

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

這次引入的東西有點多，分別用于繪圖，添加地圖投影，添加地圖特征(海岸線，湖泊等等)，設(shè)置地理坐標(biāo)(XX°E)，設(shè)置正確的中國國境線(默認(rèn)的國境線把一些沖突地區(qū)偷吃了)

fig2 = plt.figure(figsize=(15,15))
proj = ccrs.PlateCarree(central_longitude=115) 
#設(shè)置一個圓柱投影坐標(biāo)，中心經(jīng)度115°E
leftlon, rightlon, lowerlat, upperlat = (70,140,15,55)
#設(shè)置地圖邊界范圍

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)

通過函數(shù)的字面意思應(yīng)該可以理解大部分語句的意思，'bou2_4l.dbf'是中國國界線shp文件，氣象家園有下載。通過類似的方法，縮小子圖，覆蓋在大圖上便可實現(xiàn)南海的繪制。九段線同樣有對應(yīng)的shp文件。
EOF的另一部分便是PC序列的繪制，紅正藍(lán)負(fù)的實現(xiàn)通過如下循環(huán)換實現(xiàn)：

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即為折線圖