基于WRFOUT计算相对涡度，绝对涡度，位涡并可视化

发布于

版本：python3.7 数据：wrfout模拟数据 核心代码：metpy.calc.vorticity

前言

涡度是流体力学中的一个重要概念，用于描述流体运动中的旋转性质。它在研究流体动力学、天气现象等领域有着广泛的应用。

下面将展示如何从WRFOUT数据中计算相对涡度，绝对涡度，位涡及其可视化

相对涡度

实际上我们天气学所用的相对涡度应该称之为：相对涡度的垂直分量

导入计算与可视化库

代码语言：javascript
复制
代码语言：javascript
复制
from wrf import uvmet, to_np, getvar, interplevel, smooth2d, get_cartopy, cartopy_xlim, cartopy_ylim, latlon_coords
import numpy as np
from netCDF4 import Dataset
import xarray as xr
from metpy.units import units
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
from matplotlib.colors import from_levels_and_colors
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.io.shapereader as shpreader
from cartopy.feature import NaturalEarthFeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from cartopy.mpl.gridliner import LATITUDE_FORMATTER, LONGITUDE_FORMATTER
import cmaps
from glob import glob
import metpy.calc as mpcalc
from metpy.units import units
import metpy.constants as constants
代码语言：javascript
复制

提取所需变量

计算相对涡度所用的metpy.calc.vorticity的格式是

代码语言：javascript
复制
metpy.calc.vorticity(u, v, *, dx=None, dy=None, x_dim=-1, y_dim=-2, parallel_scale=None, meridional_scale=None, latitude=None, longitude=None, crs=None)

因此我们需要从wrfout文件中提取 u,v ,而p是插值需要的参数

代码语言：javascript
复制
代码语言：javascript
复制
# 用 netCDF4 包读取 WRF 模拟数据
wrf_file = Dataset('/home/mw/input/typhoon9537/wrfout_d01_2019-08-08_23_00_00')
# 提取要素
lon = getvar(wrf_file, 'lon', timeidx=0)
lat = getvar(wrf_file, 'lat', timeidx=0)
u = getvar(wrf_file, 'ua', timeidx=0)
v = getvar(wrf_file, 'va', timeidx=0)
p = getvar(wrf_file, 'pressure', timeidx=0)
# 提取WRF模拟的经纬度数组
lats, lons = latlon_coords(u)
# 提取WRF模拟的地图投影
wrf_proj = get_cartopy(u)
代码语言：javascript
复制

计算850hPa相对涡度

代码语言：javascript
复制
代码语言：javascript
复制
u850 = interplevel(u, p, 850)
v850 = interplevel(v, p, 850)

z = getvar(wrf_file, 'z', units="dm", timeidx=0)

z850 = interplevel(z, p, 850)


dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat)

vor = mpcalc.vorticity(u850,v850,dx=dx,dy=dy)

计算得到的相对涡度的单位是1/second，通常在绘图时会乘个1e5

绘制850hPa相对涡度分布图

绝对涡度

绝对涡度等于相对涡度加行星涡度f（也是垂直分量）
wrfpython可以直接使用getvar函数提取，变量名是avo

绝对涡度获取

代码语言：javascript
复制
代码语言：javascript
复制
avo = getvar(wrf_file, 'avo', timeidx=0)

avo850 = interplevel(avo, p, 850)

avo850

由wrfpython内部计算得到的绝对涡度，其单位是10-5 s-1，所以在绘图时不需要乘1e5

绘制850hPa绝对涡度分布图

位涡

罗斯贝提出的一个类似位温的用于垂直涡度的概念 ,其公式为

位涡数据获取

代码语言：javascript
复制
代码语言：javascript
复制
pvo = getvar(wrf_file, 'pvo', timeidx=0)
pvo850 = interplevel(pvo, p, 850)
pvo850

绘制850hPa位涡分布图

验证相对涡度计算结果:使用avo减去利用metpy计算的行星涡度的垂直分量f

计算相对涡度

代码语言：javascript
复制
代码语言：javascript
复制
f = mpcalc.coriolis_parameter(np.deg2rad(to_np(lats))).to(units('1/sec'))
vor1= avo850*1e-5*units('1/sec') -f
代码语言：javascript
复制

绘制850hPa相对涡度分布图

大致上相对涡度分布一致，用它们的差值作图还是有点区别

代码语言：javascript
复制
代码语言：javascript
复制
diff=vor-vor1
diff.plot()

<matplotlib.collections.QuadMesh at 0x7f11fbdc3c10>

可见差别较小，使用metpy计算的结果可信

完整代码与文件在这里，文件在注册社区账号点击左侧文件标识可下载，代码需要右上角在线运行

#file#import#可视化#数据#数组