前言
没想到食堂又出现小龙虾的尾巴，经理惦记上了捏
有读者留言想要知道怎么处理wrf的垂直速度，故写一个
首先关于上升的有两个变量，一个是wa，官网的描述是W-component of Wind on Mass Points 单位是m/s
这应该是读者关心的变量
另一个则是omega（dp/dt）,单位是Pa/s，具体内容翻开天气学原理和方法p120，小编天气学很菜就不多说了
气象家园的帖子有说，链接是https://bbs.06climate.com/forum.php?mod=viewthread&tid=57957&highlight=omega
使用omega是p坐标下的铅直速度速度，单位是hpa/s，omega=dp/dt,负数表示上升，正数表示下沉运动，
由于omega和v值数量级差太多，故而乘以-100，
w是z坐标下的垂直速度，单位是m/s，w=dz/dt,omega=-ρgw，天气动力学书中有此公式
在wrfPython中变量直接用getvar获取即可
直接上代码
台风利奇马风速剖面分布
In [12]:
import numpy as npimport matplotlib.pyplot as pltfrom matplotlib.cm import get_cmapimport cartopy.crs as crsfrom cartopy.feature import NaturalEarthFeaturefrom netCDF4 import Datasetfrom wrf import to_np, getvar, CoordPair, vertcross# Open the NetCDF filefilename = "/home/mw/input/typhoon9537/wrfout_d01_2019-08-09_03_00_00"ncfile = Dataset(filename)# Extract the model height and wind speedz = getvar(ncfile, "z")wspd = getvar(ncfile, "uvmet_wspd_wdir")[0,:]# Create the start point and end point for the cross sectionstart_point = CoordPair(lat=24, lon=120)end_point = CoordPair(lat=32, lon=128)# Compute the vertical cross-section interpolation. Also, include the# lat/lon points along the cross-section.wspd_cross = vertcross(wspd, z, wrfin=ncfile, start_point=start_point,end_point=end_point, latlon=True, meta=True)# Create the figurefig = plt.figure(figsize=(20,12))ax = plt.axes()# Make the contour plotwspd_contours = ax.contourf(to_np(wspd_cross), cmap=get_cmap("jet"))# Add the color barplt.colorbar(wspd_contours, ax=ax)# Set the x-ticks to use latitude and longitude labels.coord_pairs = to_np(wspd_cross.coords["xy_loc"])x_ticks = np.arange(coord_pairs.shape[0])x_labels = [pair.latlon_str(fmt="{:.2f}, {:.2f}")for pair in to_np(coord_pairs)]ax.set_xticks(x_ticks[::20])ax.set_xticklabels(x_labels[::20], rotation=45, fontsize=8)# Set the y-ticks to be height.vert_vals = to_np(wspd_cross.coords["vertical"])v_ticks = np.arange(vert_vals.shape[0])ax.set_yticks(v_ticks[::20])ax.set_yticklabels(vert_vals[::20], fontsize=8)# Set the x-axis and y-axis labelsax.set_xlabel("Latitude, Longitude", fontsize=12)ax.set_ylabel("Height (m)", fontsize=12)plt.title("Vertical Cross Section of Wind Speed (m/s)")plt.show()
/tmp/ipykernel_53/2475736227.py:32: MatplotlibDeprecationWarning: The get_cmap function was deprecated in Matplotlib 3.7 and will be removed two minor releases later. Use ``matplotlib.colormaps[name]`` or ``matplotlib.colormaps.get_cmap(obj)`` instead.
wspd_contours = ax.contourf(to_np(wspd_cross), cmap=get_cmap("jet"))
台风利奇马WA剖面
In [13]:
import numpy as npimport matplotlib.pyplot as pltfrom matplotlib.cm import get_cmapimport cartopy.crs as crsfrom cartopy.feature import NaturalEarthFeaturefrom netCDF4 import Datasetfrom wrf import to_np, getvar, CoordPair, vertcross# Open the NetCDF filefilename = "/home/mw/input/typhoon9537/wrfout_d01_2019-08-09_03_00_00"ncfile = Dataset(filename)# Extract the model height and wind speedz = getvar(ncfile, "z")wa = getvar(ncfile, "wa")# Create the start point and end point for the cross sectionstart_point = CoordPair(lat=24, lon=120)end_point = CoordPair(lat=32, lon=128)# Compute the vertical cross-section interpolation. Also, include the# lat/lon points along the cross-section.wa_cross = vertcross(wa, z, wrfin=ncfile, start_point=start_point,end_point=end_point, latlon=True, meta=True)# Create the figurefig = plt.figure(figsize=(20,12))ax = plt.axes()# Make the contour plotwa_contours = ax.contourf(to_np(wa_cross), cmap=get_cmap("jet"))# Add the color barplt.colorbar(wa_contours, ax=ax)# Set the x-ticks to use latitude and longitude labels.coord_pairs = to_np(wspd_cross.coords["xy_loc"])x_ticks = np.arange(coord_pairs.shape[0])x_labels = [pair.latlon_str(fmt="{:.2f}, {:.2f}")for pair in to_np(coord_pairs)]ax.set_xticks(x_ticks[::20])ax.set_xticklabels(x_labels[::20], rotation=45, fontsize=8)# Set the y-ticks to be height.vert_vals = to_np(wa_cross.coords["vertical"])v_ticks = np.arange(vert_vals.shape[0])ax.set_yticks(v_ticks[::20])ax.set_yticklabels(vert_vals[::20], fontsize=8)# Set the x-axis and y-axis labelsax.set_xlabel("Latitude, Longitude", fontsize=12)ax.set_ylabel("Height (m)", fontsize=12)plt.title("Vertical Cross Section of wa (m/s)")plt.show()
/tmp/ipykernel_53/3001559905.py:32: MatplotlibDeprecationWarning: The get_cmap function was deprecated in Matplotlib 3.7 and will be removed two minor releases later. Use ``matplotlib.colormaps[name]`` or ``matplotlib.colormaps.get_cmap(obj)`` instead.
wa_contours = ax.contourf(to_np(wa_cross), cmap=get_cmap("jet"))
台风利奇马OMEGA剖面
In [14]:
import numpy as npimport matplotlib.pyplot as pltfrom matplotlib.cm import get_cmapimport cartopy.crs as crsfrom cartopy.feature import NaturalEarthFeaturefrom netCDF4 import Datasetfrom wrf import to_np, getvar, CoordPair, vertcross# Open the NetCDF filefilename = "/home/mw/input/typhoon9537/wrfout_d01_2019-08-09_03_00_00"ncfile = Dataset(filename)# Extract the model height and wind speedz = getvar(ncfile, "z")omega = getvar(ncfile, "omega")# Create the start point and end point for the cross sectionstart_point = CoordPair(lat=24, lon=120)end_point = CoordPair(lat=32, lon=128)# Compute the vertical cross-section interpolation. Also, include the# lat/lon points along the cross-section.omega_cross = vertcross(omega, z, wrfin=ncfile, start_point=start_point,end_point=end_point, latlon=True, meta=True)# Create the figurefig = plt.figure(figsize=(20,12))ax = plt.axes()# Make the contour plotomega_contours = ax.contourf(to_np(omega_cross), cmap=get_cmap("jet"))# Add the color barplt.colorbar(omega_contours, ax=ax)# Set the x-ticks to use latitude and longitude labels.coord_pairs = to_np(omega_cross.coords["xy_loc"])x_ticks = np.arange(coord_pairs.shape[0])x_labels = [pair.latlon_str(fmt="{:.2f}, {:.2f}")for pair in to_np(coord_pairs)]ax.set_xticks(x_ticks[::20])ax.set_xticklabels(x_labels[::20], rotation=45, fontsize=8)# Set the y-ticks to be height.vert_vals = to_np(omega_cross.coords["vertical"])v_ticks = np.arange(vert_vals.shape[0])ax.set_yticks(v_ticks[::20])ax.set_yticklabels(vert_vals[::20], fontsize=8)# Set the x-axis and y-axis labelsax.set_xlabel("Latitude, Longitude", fontsize=12)ax.set_ylabel("Height (m)", fontsize=12)plt.title("Vertical Cross Section of omega (pa/s)")plt.show()
/tmp/ipykernel_53/4265223123.py:32: MatplotlibDeprecationWarning: The get_cmap function was deprecated in Matplotlib 3.7 and will be removed two minor releases later. Use ``matplotlib.colormaps[name]`` or ``matplotlib.colormaps.get_cmap(obj)`` instead.
omega_contours = ax.contourf(to_np(omega_cross), cmap=get_cmap("jet"))
小结
1. 当然大家使用时注意一下wa和omega数值上是反的
omega>0的时候是下降，反之是上升
2. 还有就是wa在普通过程中数值是非常小的，能有0.1m/s算是十分大了。
通常会乘个100。台风就不乘了
3. 想在剖面加风矢量箭头可以参考教程 https://www.heywhale.com/mw/project/618e8ca486227300178d33df3.
