WRFOUT 涡度平流和温度平流计算与可视化

发布于

前言

涡度平流和温度平流是两种常见的气象诊断量，可以帮助我们更好地理解大气运动和热力学过程。 以下代码将计算上述气象诊断量并可视化。

WRFOUT 相对涡度平流

单位：s－2；量级为10^-10

表征由水平风引起的涡度输送，其中相对涡度平流的作用是使槽脊移动。高空槽前的正涡度平流可引起辐散，槽后的负涡度平流可引起辐合。

导入库

代码语言：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
import metpy.constants as constants
import scipy.ndimage as ndimage
代码语言：javascript
复制

相对涡度数据获取与处理

代码语言：javascript
复制
代码语言：javascript
复制
ncfile = Dataset('/home/mw/input/typhoon9537/wrfout_d01_2019-08-08_23_00_00')

p = getvar(ncfile, "pressure")

z = getvar(ncfile, "z", units="dm")

ua = getvar(ncfile, "ua", units="kt")

va = getvar(ncfile, "va", units="kt")

lon = getvar(ncfile, 'lon')

lat = getvar(ncfile, 'lat')

wspd = getvar(ncfile, "wspd_wdir", units="kt")[0,:]


提前interpolate减少loop绘图时间


ht_500 = interplevel(z, p, 500)

u_500 = interplevel(ua, p, 500)* units('m/s')

v_500 = interplevel(va, p, 500)* units('m/s')

wspd_500 = interplevel(wspd, p, 500)


avor = getvar(ncfile, 'avo', meta=False)

avor = interplevel(avor, p, 500)* units('1/s')

#高斯滤波处理

avor = ndimage.gaussian_filter(avor, sigma=3, order=0) * units('1/s')

#获得行星涡度

f = mpcalc.coriolis_parameter(np.deg2rad(to_np(lat))).to(units('1/sec'))

#计算相对涡度平流

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

adv = mpcalc.advection(avor*1e-5-f, u_500, v_500, dx=dx, dy=dy)

adv
代码语言：javascript
复制

Out[52]:

Magnitude

[[-1.2625431230062143e-09 -2.218037666702117e-09 -2.8729166506779293e-09 ... -2.1899866236972013e-08 -1.9723836456431697e-08 -1.5885956880594937e-08] [-1.5420121005685506e-09 -2.273706928282503e-09 -2.5385447798787676e-09 ... -4.100087479991979e-08 -4.0626488474781005e-08 -3.698737711448127e-08] [-1.8081066592293629e-09 -2.285932193097694e-09 -2.101475473748282e-09 ... -6.082706582541257e-08 -6.085030359714656e-08 -5.5982852395177244e-08] ... [3.810899411929663e-09 1.1687452609468359e-08 1.8845908895573483e-08 ... -3.24593680999357e-08 -4.511414890635478e-08 -4.5838578525088425e-08] [3.6971627837177484e-09 1.0667419877081673e-08 1.6601254825007813e-08 ... -2.191717542269422e-08 -3.528844328269863e-08 -3.744606815253251e-08] [3.845407496543239e-09 9.992944568984933e-09 1.504098584876701e-08 ... -4.490564614940227e-09 -1.334922737216065e-08 -1.9087568578885183e-08]]

Units

1/second2

绘制500hPa相对涡度分布图

代码语言：javascript
复制
代码语言：javascript
复制
# 准备映射投影

cart_proj = ccrs.PlateCarree()


创建图形并设置大小


fig = plt.figure(figsize=(12,9))

ax = plt.axes(projection=cart_proj)


绘制填色,预先定义层级间隔


levels = np.arange(-1000., 1000., 20.)

cmap = 'NCV_jaisnd'

abc = ax.pcolormesh(to_np(lon), to_np(lat), adv,

cmap=cmaps.BlueDarkRed18,# vmin=-30, vmax=30,

transform=ccrs.PlateCarree(),alpha=0.5)


添加颜色条


cb = plt.colorbar(abc, orientation="horizontal", pad=0.05)


#限定范围

ax.set_xlim(cartopy_xlim(ht_500))


ax.set_ylim(cartopy_ylim(ht_500))


ax.coastlines()

ax.set_extent([118, 128, 22, 32], crs=ccrs.PlateCarree())

ax.set_title("500 MB Height (dm), vorticity advection")

plt.show()
代码语言：javascript
复制

ERA5再分析数据绘制涡度平流

相对涡度数据获取与处理

代码语言：javascript
复制
代码语言：javascript
复制
# 加载ERA5数据集

ds = xr.open_dataset('/home/mw/input/vor3873/vor_20190809.nc')

ds1 = xr.open_dataset('/home/mw/input/ERA5_Lekima4742/ERA5_Lekima.nc')

time = ds.time


选择数据并进行裁剪


lon_min, lon_max = 70, 140  # 经度范围

lat_min, lat_max = 60, 15  # 纬度范围


u500re = ds1.u.sel(time=time[23], level=500, longitude=slice(lon_min, lon_max), latitude=slice(lat_min, lat_max))

v500re = ds1.v.sel(time=time[23], level=500, longitude=slice(lon_min, lon_max), latitude=slice(lat_min, lat_max))

vo = ds.vo.sel(time=time[23],level=500,longitude=slice(lon_min, lon_max), latitude=slice(lat_min, lat_max))

vo
代码语言：javascript
复制

绘制500hPa相对涡度分布图

代码语言：javascript
复制
代码语言：javascript
复制
# 准备映射投影

cart_proj = ccrs.PlateCarree()


创建图形并设置大小


fig = plt.figure(figsize=(12,9))

ax = plt.axes(projection=cart_proj)


绘制填色,预先定义层级间隔


#levels = np.arange(-1000., 1000., 20.)

#cmap = 'NCV_jaisnd'

bdc = ax.pcolormesh(to_np(lon1), to_np(lat1), adv1,

cmap=cmaps.BlueDarkRed18, #vmin=-30, vmax=30,

transform=ccrs.PlateCarree(),alpha=0.5)


添加颜色条


cb = plt.colorbar(bdc, orientation="horizontal", pad=0.05)

ax.coastlines()

#限定范围

ax.set_extent([118, 128, 22, 32], crs=ccrs.PlateCarree())


ax.set_title("500 MB Height (dm), vorticity advection")

plt.show()
代码语言：javascript
复制
代码语言：javascript
复制

WRFOUT 位涡平流

位涡数据获取与处理

代码语言：javascript
复制
代码语言：javascript
复制
pvo = getvar(ncfile, 'pvo')

pvo = interplevel(pvo, p, 500)

注意，此处wrfout的位涡单位是PVU，这与下方era5数据的单位不同

代码语言：javascript
复制
代码语言：javascript
复制
pvo = ndimage.gaussian_filter(pvo, sigma=3, order=0)

pvoadv = mpcalc.advection(pvo, u_500, v_500, dx=dx, dy=dy)

pvoadv
代码语言：javascript
复制

Magnitude

[[-9.08094373792328e-06 -1.667895758996559e-05 -2.2075712824100678e-05 ... -0.00019570278548612473 -0.00018182904426596296 -0.000151876006159788] [-1.3596305101121857e-05 -1.971777285013539e-05 -2.257579169818516e-05 ... -0.000372265612215262 -0.000371158390938355 -0.00034213397237166387] [-1.7466678992253744e-05 -2.1820846610943972e-05 -2.1539021684347855e-05 ... -0.000548509724308836 -0.0005492788795706235 -0.0005102732673498781] ... [3.503358621261398e-05 0.00010836654620009474 0.00017538392407312573 ... -0.0002009477336337225 -0.0002791951364593981 -0.00028337437471440607] [3.444130022446092e-05 9.974574082071261e-05 0.00015552461786049978 ... -0.00013468424512472142 -0.00021723250766898444 -0.00023007488313979492] [3.6150355900239555e-05 9.401744065303945e-05 0.0001415819941154479 ... -2.5650905689909083e-05 -8.05448927333894e-05 -0.00011643640401079955]]

Units

1/second

绘制500hPa位涡平流分布图

代码语言：javascript
复制
代码语言：javascript
复制
# 准备映射投影

cart_proj = ccrs.PlateCarree()


创建图形并设置大小


fig = plt.figure(figsize=(12,9))

ax = plt.axes(projection=cart_proj)


绘制填色,预先定义层级间隔


#levels = np.arange(-1000., 1000., 20.)

#cmap = 'NCV_jaisnd'

cbd = ax.pcolormesh(to_np(lon), to_np(lat), pvoadv,

cmap=cmaps.BlueDarkRed18, #vmin=-30, vmax=30,

transform=ccrs.PlateCarree())


添加颜色条


cb = plt.colorbar(cbd, orientation="horizontal", pad=0.05)


#限定范围

ax.set_xlim(cartopy_xlim(ht_500))


ax.set_ylim(cartopy_ylim(ht_500))


读取地图边界数据


ax.coastlines()

ax.set_extent([118, 128, 22, 32], crs=ccrs.PlateCarree())

ax.set_title("500 MB Height (dm),potential vorticity advection")

plt.show()
代码语言：javascript
复制

ERA5 再分析数据绘制位涡平流

位涡数据获取与处理

代码语言：javascript
复制
代码语言：javascript
复制
pv = ds.pv.sel(time=time[23],level=500,longitude=slice(lon_min, lon_max), latitude=slice(lat_min, lat_max))

pv
代码语言：javascript
复制
lon2= pv.longitude

lat2=pv.latitude

d2x,d2y =  mpcalc.lat_lon_grid_deltas(lon2, lat2)

pvadv = mpcalc.advection(pv, u500re, v500re, dx=d2x, dy=d2y)

pvadv
代码语言：javascript
复制

绘制500hPa位涡平流分布图

代码语言：javascript
复制
代码语言：javascript
复制
# 准备映射投影
cart_proj = ccrs.PlateCarree()

创建图形并设置大小


fig = plt.figure(figsize=(12,9))

ax = plt.axes(projection=cart_proj)


绘制填色,预先定义层级间隔


#levels = np.arange(-1000., 1000., 20.)

#cmap = 'NCV_jaisnd'

bdc = ax.pcolormesh(to_np(lon1), to_np(lat1), pvadv,

cmap=cmaps.BlueDarkRed18, #vmin=-30, vmax=30,

transform=ccrs.PlateCarree(),alpha=0.5)


添加颜色条


cb = plt.colorbar(bdc, orientation="horizontal", pad=0.05)


ax.coastlines()

ax.set_extent([118, 128, 22, 32], crs=ccrs.PlateCarree())

ax.set_title("500 MB Height (dm),potential vorticity advection")

plt.show()
代码语言：javascript
复制

温度平流

单位K•s-1；量级为10^-5
　　温度的冷暖平流是表明大气斜压性的一种度量，大尺度天气系统的发生发展均与之有关。

温度数据获取与平流化处理

代码语言：javascript
复制
代码语言：javascript
复制
temp = getvar(ncfile, 'temp')

temp_500 = interplevel(temp, p, 500)-273.15

tempadv = mpcalc.advection(temp_500, u_500, v_500, dx=dx, dy=dy)
代码语言：javascript
复制

绘制500hPa温度平流分布图

代码语言：javascript
复制
代码语言：javascript
复制
# 创建图形并设置大小

fig = plt.figure(figsize=(12,9))

ax = plt.axes(projection=cart_proj)


绘制填色,预先定义层级间隔


levels = np.arange(-1000., 1000., 20.)

adc = ax.pcolormesh(to_np(lon), to_np(lat), tempadv,

cmap=cmaps.BlueDarkRed18, #vmin=-30, vmax=30,

transform=ccrs.PlateCarree(),alpha=0.5)


添加颜色条


cb = plt.colorbar(adc, orientation="horizontal", pad=0.05)


#限定范围

ax.set_xlim(cartopy_xlim(ht_500))


ax.set_ylim(cartopy_ylim(ht_500))


读取地图边界数据


ax.coastlines()

ax.set_extent([118, 128, 22, 32], crs=ccrs.PlateCarree())

ax.set_title("500 MB Height (dm), temperature advection *10^5")

plt.show()
代码语言：javascript
复制

#import#地图#可视化#流计算#数据