基于WRFOUT计算相对涡度，绝对涡度，位涡并可视化
版本：python3.7
数据：wrfout模拟数据
核心代码：metpy.calc.vorticity
前言
涡度是流体力学中的一个重要概念，用于描述流体运动中的旋转性质。它在研究流体动力学、天气现象等领域有着广泛的应用。
下面将展示如何从WRFOUT数据中计算相对涡度，绝对涡度，位涡及其可视化
相对涡度
实际上我们天气学所用的相对涡度应该称之为：相对涡度的垂直分量
导入计算与可视化库
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
提取所需变量
计算相对涡度所用的metpy.calc.vorticity的格式是
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是插值需要的参数
# 用 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)
计算850hPa相对涡度
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
绝对涡度获取
avo = getvar(wrf_file, 'avo', timeidx=0)
avo850 = interplevel(avo, p, 850)
avo850
由wrfpython内部计算得到的绝对涡度，其单位是10-5 s-1，所以在绘图时不需要乘1e5
绘制850hPa绝对涡度分布图
位涡
罗斯贝提出的一个类似位温的用于垂直涡度的概念 ,其公式为
位涡数据获取
pvo = getvar(wrf_file, 'pvo', timeidx=0)
pvo850 = interplevel(pvo, p, 850)
pvo850
绘制850hPa位涡分布图
验证相对涡度计算结果:使用avo减去利用metpy计算的行星涡度的垂直分量f
计算相对涡度
f = mpcalc.coriolis_parameter(np.deg2rad(to_np(lats))).to(units('1/sec'))
vor1= avo850*1e-5*units('1/sec') -f
绘制850hPa相对涡度分布图
大致上相对涡度分布一致，用它们的差值作图还是有点区别
diff=vor-vor1
diff.plot()
<matplotlib.collections.QuadMesh at 0x7f11fbdc3c10>
可见差别较小，使用metpy计算的结果可信
完整代码与文件在这里，文件在注册社区账号点击左侧文件标识可下载，代码需要右上角在线运行
文章来源于微信公众号：气ython风雨
好久没看到这样的好贴了，支持一下！
页:
[1]