Python气象数据处理与绘图(3)：以EOF为例画柱状图(折线

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

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