基于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()

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

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