WRFOUT 单层水汽通量散度与整层水汽通量散度实现2.0
镜像：气象分析3.7
数据：本人模拟2022年七月西南暴雨的wrfout文件
🎓 作者：酷炫用户名
📢 版权声明：公益性质转载需联系作者本人获取授权。转载本文时，请务必文字注明“来自：和鲸社区：酷炫用户名”，并附带本项目超链接。
前言
本项目旨在实现WRF模式中的单层水汽通量散度和整层水汽通量散度计算方法。WRF（Weather Research and Forecasting）模式是一种广泛应用于天气和气候预测研究的数值模式。水汽通量散度在天气和气候研究中具有重要作用。本项目将针对WRF模式的输出数据（WRFOUT）进行处理和分析，实现单层水汽通量散度和整层水汽通量散度的计算。
在实现该功能的过程中，下面将详细介绍所采用的公式原理，并给出相应的代码示例和使用说明。同时会对计算结果进行可视化展示，以便更好地理解和分析水汽通量散度的空间分布和变化规律。
PS : 基于评论区读者提醒，优化以下部分
更正混合比转化比湿
优化地图绘制
概念简介
水汽通量散度是衡量水汽输送量变化的一个指标。
水汽通量散度表示单位时间内和单位面积上的水汽通量变化率。它可以反映水汽是否聚集或分散,更能准确地判断是否有利于对流天气的形成。
水汽通量散度的负值表示水汽聚集,这有利于对流天气的形成;反之,正值表示水汽分散,不利于对流天气。
水汽通量散度公式为
本文计算部分参考了这里
单层与整层的概念可以阅读此处
单层水汽实现
1.导入库
!pip install frykit -i https://pypi.mirrors.ustc.edu.cn/simple/
Looking in indexes: https://pypi.mirrors.ustc.edu.cn/simple/
Requirement already satisfied: frykit in /opt/conda/lib/python3.9/site-packages (0.6.9)
Requirement already satisfied: cartopy>=0.20.0 in /opt/conda/lib/python3.9/site-packages (from frykit) (0.23.0)
Requirement already satisfied: numpy>=1.20.0 in /opt/conda/lib/python3.9/site-packages (from frykit) (1.26.4)
Requirement already satisfied: pandas>=1.2.0 in /opt/conda/lib/python3.9/site-packages (from frykit) (2.0.3)
Requirement already satisfied: shapely>=1.7 in /opt/conda/lib/python3.9/site-packages (from cartopy>=0.20.0->frykit) (1.8.5.post1)
Requirement already satisfied: pyshp>=2.3 in /opt/conda/lib/python3.9/site-packages (from cartopy>=0.20.0->frykit) (2.3.1)
Requirement already satisfied: matplotlib>=3.5 in /opt/conda/lib/python3.9/site-packages (from cartopy>=0.20.0->frykit) (3.8.3)
Requirement already satisfied: packaging>=20 in /opt/conda/lib/python3.9/site-packages (from cartopy>=0.20.0->frykit) (23.2)
Requirement already satisfied: pyproj>=3.3.1 in /opt/conda/lib/python3.9/site-packages (from cartopy>=0.20.0->frykit) (3.4.1)
Requirement already satisfied: python-dateutil>=2.8.2 in /opt/conda/lib/python3.9/site-packages (from pandas>=1.2.0->frykit) (2.8.2)
Requirement already satisfied: tzdata>=2022.1 in /opt/conda/lib/python3.9/site-packages (from pandas>=1.2.0->frykit) (2024.1)
Requirement already satisfied: pytz>=2020.1 in /opt/conda/lib/python3.9/site-packages (from pandas>=1.2.0->frykit) (2022.1)
Requirement already satisfied: cycler>=0.10 in /opt/conda/lib/python3.9/site-packages (from matplotlib>=3.5->cartopy>=0.20.0->frykit) (0.11.0)
Requirement already satisfied: contourpy>=1.0.1 in /opt/conda/lib/python3.9/site-packages (from matplotlib>=3.5->cartopy>=0.20.0->frykit) (1.2.0)
Requirement already satisfied: kiwisolver>=1.3.1 in /opt/conda/lib/python3.9/site-packages (from matplotlib>=3.5->cartopy>=0.20.0->frykit) (1.4.2)
Requirement already satisfied: pyparsing>=2.3.1 in /opt/conda/lib/python3.9/site-packages (from matplotlib>=3.5->cartopy>=0.20.0->frykit) (3.0.9)
Requirement already satisfied: importlib-resources>=3.2.0 in /opt/conda/lib/python3.9/site-packages (from matplotlib>=3.5->cartopy>=0.20.0->frykit) (5.7.1)
Requirement already satisfied: pillow>=8 in /opt/conda/lib/python3.9/site-packages (from matplotlib>=3.5->cartopy>=0.20.0->frykit) (9.4.0)
Requirement already satisfied: fonttools>=4.22.0 in /opt/conda/lib/python3.9/site-packages (from matplotlib>=3.5->cartopy>=0.20.0->frykit) (4.33.3)
Requirement already satisfied: certifi in /opt/conda/lib/python3.9/site-packages (from pyproj>=3.3.1->cartopy>=0.20.0->frykit) (2024.2.2)
Requirement already satisfied: six>=1.5 in /opt/conda/lib/python3.9/site-packages (from python-dateutil>=2.8.2->pandas>=1.2.0->frykit) (1.16.0)
Requirement already satisfied: zipp>=3.1.0 in /opt/conda/lib/python3.9/site-packages (from importlib-resources>=3.2.0->matplotlib>=3.5->cartopy>=0.20.0->frykit) (3.8.0)
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
2. 读取数据
# 用 netCDF4 包读取 WRF 模拟数据
wrf_file = Dataset('/home/mw/input/typhoon9537/wrfout_d01_2019-08-08_23_00_00')
# 提取要素
lon = getvar(wrf_file, 'lon')
lat = getvar(wrf_file, 'lat')
u = getvar(wrf_file, 'ua')
v = getvar(wrf_file, 'va')
q = getvar(wrf_file, 'QVAPOR')
z = getvar(wrf_file, 'z')
# 提取WRF模拟的经纬度数组
lats, lons = latlon_coords(u)
# 提取WRF模拟的地图投影
wrf_proj = get_cartopy(u)
3. 散度计算过程
p = getvar(wrf_file, 'pressure')
u850 = interplevel(u, p, 850)
v850 = interplevel(v, p, 850)
q850 = interplevel(q, p, 850)
z850 = interplevel(z, p, 850)
## 混合比转为比湿
q850 = q850/(1+q850)
dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat)
# 计算水汽通量散度
q_flux_divergence = mpcalc.divergence(
to_np(
u850 *
q850/constants.g),
to_np(
v850 *
q850/constants.g),
dx=to_np(dx),
dy=to_np(dy))
q_flux_divergence
| Magnitude | [[5.416015133860801e-09 2.5922963398730183e-08 3.294845301376941e-08 ... |
|---|---|
| Units | 1/meter |
4. 绘图
import frykit.plot as fplt
# 创建画布
fig = plt.figure(figsize=(10, 8))
# 设置地图投影
ax = plt.axes(projection=wrf_proj)
# 设置地图范围
ax.set_xlim(cartopy_xlim(u))
ax.set_ylim(cartopy_ylim(u))
# 绘制分布（用contourf方法实现等值线填充）
plt.contourf(to_np(lons),
to_np(lats),
to_np(q_flux_divergence),
levels=np.arange(-50, 50, 5),
extend='both',
transform=ccrs.PlateCarree(),
cmap=cmaps.MPL_BuPu_r)
# 添加colorbar
cbar = plt.colorbar(ax=ax, extend='both', shrink=1)
cbar.set_label('flux', fontdict={'size': 20})
cbar.ax.tick_params(labelsize=20)
fplt.add_cn_province(ax)
# 设置刻度标签的字体大小
plt.tick_params(labelsize=15)
# 添加标题
plt.title(str(u.Time.values)[0:16]+'(UTC)', loc='left', fontsize=20)
plt.show()
整层水汽
测试一次性插值多个高度层变量
levels = [1000,925,850,700,600,500,400,300]
zinterp = interplevel(z, p, levels)
zinterp
u_interp = interplevel(u, p, levels)
v_interp = interplevel(v, p, levels)
q_interp = interplevel(q, p, levels)
## 混合比转换比湿
q_interp = q_interp/(1+q_interp)
lev = u_interp.level
# 计算水汽通量散度
qu = u_interp *q_interp/constants.g
qv = v_interp *q_interp/constants.g
q_flux_divergence_all = np.zeros((lev.shape[0],lat.shape[0],lon.shape[1]))
for i in range(lev.shape[0]):
q_flux_divergence_all[i] = mpcalc.divergence(u = to_np(qu[i]),v = to_np(qv[i]),dx = to_np(dx) ,dy = to_np(dy),x_dim=-1, y_dim=-2)
# 将 q_flux_divergence_all 中的 NaN 值替换为 0
q_flux_divergence_all = np.nan_to_num(q_flux_divergence_all, nan=0)
total_div_qv = np.trapz(q_flux_divergence_all,lev,axis=0)
total_div_qv[2]
array([-1.23045415e-05, -1.92274235e-06, 3.30437385e-06, 2.36010008e-06,
-1.92641520e-06, -1.62218247e-06, -4.20900298e-06, -8.60500777e-06,
-1.35983031e-05, -1.27191566e-05, -9.48728556e-06, -7.97110291e-06,
-7.43931893e-06, -7.10214819e-06, -7.65573923e-06, -8.71555968e-06,
-1.02940113e-05, -1.15697656e-05, -1.23363635e-05, -1.28024959e-05,
-1.27916642e-05, -1.26266829e-05, -1.27376472e-05, -1.32539894e-05,
-1.43531046e-05, -1.59188653e-05, -1.79339798e-05, -1.96122505e-05,
-2.08179349e-05, -2.14130566e-05, -2.17538690e-05, -2.16238217e-05,
-2.11385941e-05, -2.01178170e-05, -1.87489909e-05, -1.68306736e-05,
-1.46711831e-05, -1.25545910e-05, -1.05596372e-05, -8.63732272e-06,
-6.80231058e-06, -5.00241246e-06, -3.54421184e-06, -2.21593321e-06,
-8.51034011e-07, 1.44932865e-07, 9.42846013e-07, 1.89380313e-06,
2.54864418e-06, 2.97459698e-06, 3.06910836e-06, 2.62002480e-06,
1.50500995e-06, -8.79092281e-07, -3.88421081e-06, -6.40408921e-06,
-7.90258950e-06, -8.79398371e-06, -9.35923824e-06, -9.06693035e-06,
# 创建画布
fig = plt.figure(figsize=(10, 8))
# 设置地图投影
ax = plt.axes(projection=wrf_proj)
# 设置地图范围
ax.set_xlim(cartopy_xlim(u))
ax.set_ylim(cartopy_ylim(u))
# 用contourf方法实现等值线填充）
plt.contourf(to_np(lons),
to_np(lats),
to_np(total_div_qv*10e5),
levels=np.arange(-500, 500, 5),
extend='both',
transform=ccrs.PlateCarree(),
cmap=cmaps.MPL_BuPu_r)
# 添加colorbar
cbar = plt.colorbar(ax=ax, extend='both', shrink=1)
cbar.set_label('flux * 10e-5', fontdict={'size': 20})
cbar.ax.tick_params(labelsize=20)
fplt.add_cn_province(ax)
# 设置刻度标签的字体大小
plt.tick_params(labelsize=15)
# 添加标题
plt.title(str(u.Time.values)[0:16]+'(UTC)', loc='left', fontsize=20)
plt.show()
最后
metpy要注意的点挺多的，什么单位，维度，格式都要对齐才能好好运行，大家只能多测试。
注意这里的整层的水汽通量散度只是简单插值几层高度层进行积分。
