python绘制风压和风向叠加图

学习目标：

掌握python绘制风压和风向叠加图

由下面两篇文章修改而来

1. https://blog.csdn.net/qiuylulu/article/details/115304020
2. https://huaweicloud.csdn.net/638079ebdacf622b8df88733.html?spm=1001.2101.3001.6650.4&utm_medium=distribute.pc_relevant.none-task-blog-2%7Edefault%7ECTRLIST%7Eactivity-4-119879764-blog-115304020.235%5Ev36%5Epc_relevant_default_base3&depth_1-utm_source=distribute.pc_relevant.none-task-blog-2%7Edefault%7ECTRLIST%7Eactivity-4-119879764-blog-115304020.235%5Ev36%5Epc_relevant_default_base3&utm_relevant_index=7

绘图结果

—

相应的绘图代码

###引入库包
#matplotlib用来做图，包含了画图的pyplot,ticker等函数
#numpy 简单计算
#cartopy 绘地图
#netCDF4 读取NC格式数据
import matplotlib.pyplot as plt
import numpy as np
import matplotlib as mpl
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.ticker as mticker
import netCDF4 as nc
import matplotlib.pyplot as plt
#此外，为了控制图内所有字体大小、轴线粗细、字体等信息一致，使图片更加美化，统一进行如下设置：
mpl.rcParams["font.family"] = 'Arial'  #默认字体类型
mpl.rcParams["mathtext.fontset"] = 'cm' #数学文字字体
mpl.rcParams["font.size"] = 12
mpl.rcParams["axes.linewidth"] = 1
#读取文件，这里用了ERA5的pressure文件来选取850hpa上的uv风场，single文件里则代表了地表的参量。
f1=nc.Dataset('C:\\Users\\huain\\Desktop\\msl.nc')#在ERA5下载的选择数据集ERA5 hourly data on single levels from 1940 to present，选择Mean sea level pressure
f2=nc.Dataset('C:\\Users\\huain\\Desktop\\pha.nc')#在ERA5下载的选择数据集ERA5 hourly data on pressure levels from 1940 to present，选择850hpa的u,v风速
#读取变量，需要用到msl（海平面气压），lat，lon，time，u，v等信息。
#注意：如果不知道变量名，可以print(f1)来看下具体有哪些变量。
#ERA5中lat是从50->30这样从北往南排列的。
#lat,lon 都是一维函数，msl是（时间，lat,lon）三维函数，u850和v850都是（时间，层数，lat，lon）的四位函数，绘图时只需要二维结果，所以后续要对数据进行挑选和加工。
msl=f1.variables["msl"][:]/100   #pa->hpa
lat=f1.variables['latitude'][:]
lon=f1.variables['longitude'][:]
time=f1.variables['time'][:]
u850=f2.variables["u"][:]
v850=f2.variables["v"][:]
#首先，要建立二维的lon2d,和lat2d，后面画填色图的时候经纬度信息要是二维的。
lon2d, lat2d = np.meshgrid(lon, lat)
msl_all=np.mean(msl[:,:,:],axis=0)
u850_aim=np.mean(u850[:,:,:],axis=0)
v850_aim=np.mean(v850[:,:,:],axis=0)
#使用的是最简单的lat,lon的投影，图片大小是10*8，这样就建立了画布，得到了图的轴ax。
proj = ccrs.PlateCarree(central_longitude=130)  #中国为左
fig = plt.figure(figsize=(10,8),dpi=550)  # 创建画布
ax = fig.subplots(1, 1, subplot_kw={'projection': proj})  # 创建子图
u_all=u850_aim
v_all=v850_aim
#-----------绘制地图-------------------------------------------
ax.add_feature(cfeature.LAND.with_scale('50m'))####添加陆地######
ax.add_feature(cfeature.COASTLINE.with_scale('50m'),lw=2)#####添加海岸线#########
ax.add_feature(cfeature.OCEAN.with_scale('50m'))######添加海洋########
#-----------添加经纬度---------------------------------------
extent=[90,180,-30,50]##经纬度范围
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=0., color='k', alpha=0.5, linestyle='--')
gl.top_labels = False ##关闭上侧坐标显示
gl.right_labels = False ##关闭右侧坐标显示
gl.xformatter = LONGITUDE_FORMATTER ##坐标刻度转换为经纬度样式
gl.yformatter = LATITUDE_FORMATTER 
gl.xlocator = mticker.FixedLocator(np.arange(extent[0], extent[1]+1, 10))
gl.ylocator = mticker.FixedLocator(np.arange(extent[2], extent[3]+1, 10))
gl.xlabel_style={'size':10}
gl.ylabel_style={'size':10}
ax.set_extent(extent)     #显示所选择的区域
#首先levels定义了填色图色标的范围，1004-1032之间，间隔1个给一个颜色。
#ax.contourf绘制了海平面气压的填色图，参数lon2d, lat2d ,msl_all都是二维的数组，colormap选择了’Spectral_r’，这个可以看个人的需求，网上有很多资源。
#ax.quiver是绘制风场的，也需要输入lon2d, lat2d，u,v的信息（此时uv为三维数组，需要挑选自己所需要的层数，变成二维的），scale用来调整画出来的风场的大小，值越大，风的线条越小。width是风线条的粗细，headwidth和headlength是箭头粗细和长度。uv和lat,lon都是每隔5个显示，是为了使风场图不那么密集，这个可以根据自己的需求进行调节。
#最重要的是，一定要加上transform=ccrs.PlateCarree()，图片中地图需要转换下，使地图和数据直接做出来的二维图对应上。
#绘制每个子图的标题，pad代表标题离轴的远近。
#返回了填色图和风场图，这样函数就完成了，后面只需要传递参数调用就行。
levels = np.arange(1004, 1032 + 1, 1)
cb = ax.contourf(lon2d,lat2d,msl_all, levels=levels, cmap='bwr',transform=ccrs.PlateCarree())
cq =ax.quiver(lon2d[::10,::10],lat2d[::10,::10],u_all[::10,::10],v_all[::10,::10],color='k',scale=200,zorder=10,width=0.003,headwidth=3,headlength=4.5,transform=ccrs.PlateCarree()) 

#箭头（q，位置左右，位置上下，要显示的风速（10m/s）,label，label在箭头右侧，以axes为坐标系，字体大小）
ax.quiverkey(cq, 0.85, 1.02, 10, '10 m/s', labelpos='E',coordinates='axes',fontproperties={'size':12})

#因为是填色图，所以一定要条件colorbar. ax=(ax[0][0],ax[0][1],ax[1][0],ax[1][1])，方向为垂直放置在cb3旁边，tick显示的范围为1004-1032，每隔4个显示。aspect是colorbar的长宽高比例，shrink是代表要显示多长，是百分比。pad是离图片的距离。
#ax.tick_params里的pad的意思是色标上的数字离色标的距离，length和width是色标tick的长断粗细。

#绘制colorbar 
cb = fig.colorbar(cb,ax=(ax),orientation='vertical',ticks=np.arange(1004,1032+4,4),aspect=30,shrink=0.9,pad=0.02)
cb.set_label('Pressure',fontsize=18)
#--------------------add text----------------------------
ax.text(0.85,0.87,'L',color='red',fontsize=15,weight='bold',va='bottom',ha='center',transform=ax.transAxes)
ax.text(0.92,0.12,'H',color='red',fontsize=15,weight='bold',va='bottom',ha='center',transform=ax.transAxes)
#标题
plt.suptitle('2023_05_12_23:00 850hPa Winds and Pressure',fontsize=14,y=0.93)
plt.savefig("Figure2.png")

---