风场剖面叠加气温等值线图

本文作者：第八星系-李智

联系邮箱：lizhi258147369@163.com

分别以流线和风向杆两种形式展示风场剖面，数据来源于ERA5再分析资料

import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
import cartopy.crs as ccrs
import cartopy.feature as cf
import cartopy.io.shapereader as shpreader
from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.ticker as mticker
from matplotlib import rcParams


filename=r'E:\北大兼职\20230310-11臭氧污染事件\亚运会同期历史数据\纵向风场\ERA520220501.nc'

f=xr.open_dataset(filename)
lon=f["longitude"][:]#读取经度，一维的
lat=f["latitude"][:]#读取纬度，一维的

#print(lon)
#print(lat)

T0=f['t'][:][:][:][:]-273.15#读取气温，4维的
T3=f['t'][:][:][:][:]-273.15#读取气温，4维的
RH_P0_L100_GLL0=f['r'][:][:][:][:]#读取相对湿度，4维的
RH_P0_L100_GLL3=f['r'][:][:][:][:]#读取相对湿度，4维的
UGRD_P0_L100_GLL0=f['u'][:][:][:][:]#读取纬向风，4维的
VGRD_P0_L100_GLL0=f['v'][:][:][:][:]#读取经向风，4维的
UGRD_P0_L100_GLL3=f['u'][:][:][:][:]#读取纬向风，4维的
VGRD_P0_L100_GLL3=f['v'][:][:][:][:]#读取经向风，4维的
print(UGRD_P0_L100_GLL0)
pres= f['level'][:]#读取气压层数，一维的

config = {
    "font.family": 'serif',
    "font.size": 7,  #
    "mathtext.fontset": 'stix',  # matplotlib渲染数学字体时使用的字体，和Times New Roman差别不大
    "font.serif": ['SimSun'],  # 宋体
    'axes.unicode_minus': False  # 处理负号，即-号
}
rcParams.update(config)

plt.rcParams['axes.unicode_minus']=False  #有的情况下，绘图负号显示不正常，此处更正

proj = ccrs.PlateCarree()
fig, axes = plt.subplots(1,2,figsize=(6,4),constrained_layout=True,dpi=600)
ax1 = plt.subplot(121)
ax2 = plt.subplot(122)

ax1.invert_yaxis()#反转纵轴，使1000hPa作为起点
ax1.set_yticks([1000,925,850,700,500])#突出我们感兴趣的气压层
ax1.set_xticks(np.arange(119,123,1))
ax1.xaxis.set_major_formatter(LongitudeFormatter())
ax1.set(ylim=(1000,500))#设置气压层绘图范围，此处我们只显示到500hPa
ax1.set_xlabel('longitude',fontsize=7)
ax1.set_ylabel('pressure(hPa)',fontsize=7)
ax1.tick_params(axis='both',which='both',labelsize=9)
##################################################################################
ac1=ax1.contourf(lon[24:40],pres[:],T0[0,:,19,24:40],   #温度等值线做底图
                 cmap='seismic',levels=np.arange(-20,20,2),extend='both',alpha=0.75)

q1 = ax1.quiver(lon[24:40],pres[:],UGRD_P0_L100_GLL0[0,:,19,24:40],VGRD_P0_L100_GLL0[0,:,19,24:40],  #创建流线图
        color="k", width=0.005, scale=50, zorder=3)
ax1.quiverkey(q1, 0.48, 0.85, U=5, angle=-90, label="5 m/s",
        labelpos="E", color="k", labelcolor="k",  coordinates='figure')
plt.colorbar(ac1,extend='both',shrink=1,label='Air temperature（℃）',pad=0.08)
ax1.axvline(121.14,ls='--',c='m')#对感兴趣的经度做出标记

ax1.set_title('2022年5月1日06时气温与风场剖面图',fontsize=8)


ax2.invert_yaxis()#反转纵轴，使1000hPa作为起点
ax2.set_yticks([1000,925,850,700,500])#突出我们感兴趣的气压层
ax2.set_xticks(np.arange(119,123,1))
ax2.xaxis.set_major_formatter(LongitudeFormatter())
ax2.set(ylim=(1000,500))#设置气压层绘图范围，此处我们只显示到500hPa
ax2.set_xlabel('longitude',fontsize=7)
ax2.set_ylabel('pressure(hPa)',fontsize=7)
ax2.tick_params(axis='both',which='both',labelsize=9)
##################################################################################
ac2=ax2.contourf(lon[24:40],pres[:],T3[3,:,19,24:40],      #温度等值线做底图
               cmap='seismic',levels=np.arange(-20,20,2),extend='both',alpha=0.75)
ax2.barbs(lon[24:40],pres[:],UGRD_P0_L100_GLL3[3,:,19,24:40],VGRD_P0_L100_GLL3[3,:,19,24:40],  #创建风向杆
         barb_increments={'half':2,'full':4,'flag':20},length=5)  #使风向杆符合中国标准

ax2.axvline(121.14,ls='--',c='m')#对感兴趣的经度做出标记
ax2.set_title('2022年5月1日23时气温与风场剖面图',fontsize=8)
plt.colorbar(ac2,extend='both',shrink=1,label='Air temperature（℃）',pad=0.08)

plt.savefig('E:/北大兼职/20230310-11臭氧污染事件/亚运会同期历史数据/纵向风场/test1.png', dpi=600, bbox_inches='tight')
plt.show()