Python气象绘图 Lambert扇形投影,以冬季北大西洋涛动(NAO)经典模态空间分布为例

nao_slp = ds.sel(time=slice('1979-12-01', '2022-02-01'), lat= slice(80, 20), lon= slice(-90, 40)).slp
DJF_sel = (nao_slp['time'].dt.month == 12) | (nao_slp['time'].dt.month == 1) | (nao_slp['time'].dt.month == 2)
nao_DJF_slp = nao_slp.sel(time=DJF_sel)
# 数据化为冬季
data_group = nao_DJF_slp.groupby('time.month')
mon1 = list(data_group)[0][1] #1月份
mon2 = list(data_group)[1][1] #2月份
mon12 = list(data_group)[2][1] #12月份
time = mon12['time']
new_time = pd.date_range(start= '1980-12-01',periods=len(time), freq='A').values
mon12['time'] = new_time
data_merge = xr.concat([mon1, mon2, mon12], dim='time')
DJF_data = data_merge.groupby('time.year').mean().rename({'year':'time'}) #冬季平均数据
lon = nao_slp.lon
lat = nao_slp.lat

2.EOF分析

n = 3 #修改n值即可得到需要的模态数
coslat = np.cos(np.deg2rad(lat.values))
wgts = np.sqrt(coslat)[..., np.newaxis]
solver = Eof(DJF_data.values, weights=wgts)
eof1 = solver.eofsAsCovariance(neofs=n)[0]
pc = solver.pcs(npcs=n, pcscaling= 1) #方差
var = solver.varianceFraction(neigs= n)

3.Lambert扇形投影绘图

首先是普通的绘图代码，注意投影proj = ccrs.LambertConforma（）

leftlon, rightlon, lowerlat, upperlat = [-90, 40, 20,80]
proj=ccrs.LambertConformal(central_longitude=(leftlon+rightlon)*0.5,
central_latitude=(lowerlat+upperlat)*0.5)
fig= plt.figure(figsize=(8,8),dpi=300)#设置画布大小
ax= fig.add_axes([0.2,0.3,0.5,0.5],projection =proj)
ax.add_feature(cfeature.COASTLINE.with_scale('110m'))
c1 = ax.contourf(lon, lat, eof1.squeeze(), zorder=0, extend='both', levels= np.arange(-3, 6.5,0.5), norm = mcolors.TwoSlopeNorm(vmin=-3, vmax = 6, vcenter=0.5), cmap= plt.cm.RdYlBu_r,transform=ccrs.PlateCarree())
position=fig.add_axes([0.2, 0.35, 0.5, 0.015])
fig.colorbar(c1,cax=position,orientation='horizontal')#,format='%.2f',)

下面是扇形绘图显示的关键部分，ax.set_boundary（）很重要（使用不多，不是很了解）

proj_to_data = ccrs.PlateCarree()._as_mpl_transform(ax) - ax.transData
rect = mpath.Path([[leftlon, lowerlat], [rightlon, lowerlat],
[rightlon, upperlat], [leftlon, upperlat], [leftlon, lowerlat]]).interpolated(50)
rect_in_target = proj_to_data.transform_path(rect)
ax.set_boundary(rect_in_target)
ax.set_xlim(rect_in_target.vertices[:,0].min(), rect_in_target.vertices[:,0].max())
ax.set_ylim(rect_in_target.vertices[:,1].min(), rect_in_target.vertices[:,1].max())

最后设置刻度的显示（这里了解不多，不能随心设置显示的刻度）

gl=ax.gridlines(draw_labels=True, x_inline=False, y_inline=False, linestyle='dashed')
gl.top_labels=False
gl.right_labels= False
gl.rotate_labels=False
gl.xlocator=ctk.LongitudeLocator(6)
gl.ylocator=ctk.LatitudeLocator(5)
gl.xformatter=ctk.LongitudeFormatter(zero_direction_label=True)
gl.yformatter=ctk.LatitudeFormatter()

最终结果rt

注：360°经度转180°的代码源于CSDN，时间久远不太记得；扇形Lambert投影设置应该是来自摸鱼。以上，原作者看到可联系。

————————————————

版权声明：本文为知乎博主「momo」的原创文章，转载请附上原文出处链接及本声明。

原文链接：https://zhuanlan.zhihu.com/p/699465671