WRFOUT 绘制台站探空图与简单分析

#

#库
from wrf import uvmet, to_np, getvar, interplevel, smooth2d, get_cartopy, cartopy_xlim, cartopy_ylim, latlon_coords,ll_to_xy
import numpy as np
from netCDF4 import Dataset
import matplotlib.pyplot as plt
import matplotlib as m
import metpy.calc as mpcalc
from metpy.plots import Hodograph, SkewT
from metpy.units import units
import metpy.calc as mpcalc
from metpy.cbook import get_test_data
from metpy.plots import add_metpy_logo, SkewT
from metpy.units import units
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1.inset_locator import inset_axes

#读取WRF输出文件
wrfin = Dataset('/home/mw/input/wrfout3385/wrfout_d02_2022-07-14_0700.nc')

#指定要提取的经纬度坐标点
lat_lon = [35, 104]

#将经纬度坐标转换为模型坐标系（x, y）
x_y = ll_to_xy(wrfin, lat_lon[0], lat_lon[1])

#提取所需变量数据
p = getvar(wrfin, "pressure", timeidx=0)[:, x_y[0], x_y[1]] * units.hPa
T = getvar(wrfin, "tc", timeidx=0)[:, x_y[0], x_y[1]] * units.degC
Td = getvar(wrfin, "td", timeidx=0)[:, x_y[0], x_y[1]] * units.degC
u = getvar(wrfin, "ua", timeidx=0)[:, x_y[0], x_y[1]] * units('m/s')
v = getvar(wrfin, "va", timeidx=0)[:, x_y[0], x_y[1]] * units('m/s')
h = getvar(wrfin, "height", timeidx=0)[:, x_y[0], x_y[1]] 

# 计算风向和风速
wind_dir = mpcalc.wind_direction(u, v)
wind_speed = mpcalc.wind_speed(u, v)

#创建9*9英寸，100dpi的画布
fig = plt.figure(figsize=(9, 9), dpi=100)

#在画布上添加SkewT对象并设置旋转角度为45度
skew = SkewT(fig, rotation=45)

#绘制温度和露点温度线
skew.plot(p, T, 'r')
skew.plot(p, Td, 'g')

#绘制风羽标志
skew.plot_barbs(p, u, v)

#设置y轴范围为1050hPa到100hPa，x轴范围为0到40度
skew.ax.set_ylim(1050, 100)
skew.ax.set_xlim(-30, 40)

#添加特殊线，如干绝热线、湿绝热线和混合比线
skew.plot_dry_adiabats()
skew.plot_moist_adiabats()
skew.plot_mixing_lines()

# 创建高空风速分析图
ax_hod = inset_axes(skew.ax, '30%', '30%', loc=1)
h = Hodograph(ax_hod, component_range=70.)
h.add_grid(increment=20)
h.plot_colormapped(u, v, wind_speed)  # Plot a line colored by wind speed

#设置x轴和y轴标签
skew.ax.set_xlabel('温度 (°C)')
skew.ax.set_ylabel('压强 (hPa)')

Text(0, 0.5, '压强 (hPa)')

cape, cin = mpcalc.cape_cin(p, T, Td, prof, which_lfc='bottom', which_el='top')
print(cape, cin)

1360.9249456932523 joule / kilogram -104.38382141033901 joule / kilogram