WRFOUT 涡度平流和温度平流计算与可视化
本帖最后由 气ython风雨 于 2024-3-9 19:02 编辑前言
涡度平流和温度平流是两种常见的气象诊断量，可以帮助我们更好地理解大气运动和热力学过程。
以下代码将计算上述气象诊断量并可视化。
WRFOUT 相对涡度平流
单位：s－2；量级为10^-10
表征由水平风引起的涡度输送，其中相对涡度平流的作用是使槽脊移动。高空槽前的正涡度平流可引起辐散，槽后的负涡度平流可引起辐合。
导入库
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
相对涡度数据获取与处理
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")
# 提前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
Out:
Magnitude 62143e-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相对涡度分布图
# 准备映射投影
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(, crs=ccrs.PlateCarree())
ax.set_title("500 MB Height (dm), vorticity advection")
plt.show()
ERA5再分析数据绘制涡度平流
相对涡度数据获取与处理
# 加载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, level=500, longitude=slice(lon_min, lon_max), latitude=slice(lat_min, lat_max))
v500re = ds1.v.sel(time=time, level=500, longitude=slice(lon_min, lon_max), latitude=slice(lat_min, lat_max))
vo = ds.vo.sel(time=time,level=500,longitude=slice(lon_min, lon_max), latitude=slice(lat_min, lat_max))
vo
绘制500hPa相对涡度分布图
# 准备映射投影
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(, crs=ccrs.PlateCarree())
ax.set_title("500 MB Height (dm), vorticity advection")
plt.show()
WRFOUT 位涡平流
位涡数据获取与处理
pvo = getvar(ncfile, 'pvo')
pvo = interplevel(pvo, p, 500)
注意，此处wrfout的位涡单位是PVU，这与下方era5数据的单位不同pvo = ndimage.gaussian_filter(pvo, sigma=3, order=0)
pvoadv = mpcalc.advection(pvo, u_500, v_500, dx=dx, dy=dy)
pvoadv
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位涡平流分布图
# 准备映射投影
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(, crs=ccrs.PlateCarree())
ax.set_title("500 MB Height (dm),potential vorticity advection")
plt.show()
ERA5 再分析数据绘制位涡平流
位涡数据获取与处理
pv = ds.pv.sel(time=time,level=500,longitude=slice(lon_min, lon_max), latitude=slice(lat_min, lat_max))
pv
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
绘制500hPa位涡平流分布图
# 准备映射投影
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(, crs=ccrs.PlateCarree())
ax.set_title("500 MB Height (dm),potential vorticity advection")
plt.show()
温度平流
单位K•s-1；量级为10^-5
温度的冷暖平流是表明大气斜压性的一种度量，大尺度天气系统的发生发展均与之有关。
温度数据获取与平流化处理
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)
绘制500hPa温度平流分布图
# 创建图形并设置大小
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(, crs=ccrs.PlateCarree())
ax.set_title("500 MB Height (dm), temperature advection *10^5")
plt.show()
完整代码与文件在这里
文章来源于微信公众号：气ython风雨
页:
[1]