如何绘制wrfout文件的垂直速度变量

用户11172986

发布于 2024-06-20 19:27:46

760

发布于 2024-06-20 19:27:46

文章被收录于专栏：气python风雨气python风雨

前言

没想到食堂又出现小龙虾的尾巴，经理惦记上了捏

有读者留言想要知道怎么处理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 np
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
import cartopy.crs as crs
from cartopy.feature import NaturalEarthFeature
from netCDF4 import Dataset

from wrf import to_np, getvar, CoordPair, vertcross

# Open the NetCDF file
filename = "/home/mw/input/typhoon9537/wrfout_d01_2019-08-09_03_00_00"
ncfile = Dataset(filename)

# Extract the model height and wind speed
z = getvar(ncfile, "z")
wspd =  getvar(ncfile, "uvmet_wspd_wdir")[0,:]

# Create the start point and end point for the cross section
start_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 figure
fig = plt.figure(figsize=(20,12))
ax = plt.axes()

# Make the contour plot
wspd_contours = ax.contourf(to_np(wspd_cross), cmap=get_cmap("jet"))

# Add the color bar
plt.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 labels
ax.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 np
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
import cartopy.crs as crs
from cartopy.feature import NaturalEarthFeature
from netCDF4 import Dataset

from wrf import to_np, getvar, CoordPair, vertcross

# Open the NetCDF file
filename = "/home/mw/input/typhoon9537/wrfout_d01_2019-08-09_03_00_00"
ncfile = Dataset(filename)

# Extract the model height and wind speed
z = getvar(ncfile, "z")
wa =  getvar(ncfile, "wa")

# Create the start point and end point for the cross section
start_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 figure
fig = plt.figure(figsize=(20,12))
ax = plt.axes()

# Make the contour plot
wa_contours = ax.contourf(to_np(wa_cross), cmap=get_cmap("jet"))

# Add the color bar
plt.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 labels
ax.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 np
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
import cartopy.crs as crs
from cartopy.feature import NaturalEarthFeature
from netCDF4 import Dataset

from wrf import to_np, getvar, CoordPair, vertcross

# Open the NetCDF file
filename = "/home/mw/input/typhoon9537/wrfout_d01_2019-08-09_03_00_00"
ncfile = Dataset(filename)

# Extract the model height and wind speed
z = getvar(ncfile, "z")
omega =  getvar(ncfile, "omega")

# Create the start point and end point for the cross section
start_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 figure
fig = plt.figure(figsize=(20,12))
ax = plt.axes()

# Make the contour plot
omega_contours = ax.contourf(to_np(omega_cross), cmap=get_cmap("jet"))

# Add the color bar
plt.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 labels
ax.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.

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2024-06-08，如有侵权请联系 cloudcommunity@tencent.com 删除

import