WRF散度计算步骤及Python可视化实现2.0

用户11172986

发布于 2025-04-22 04:59:13

11100

代码可运行

文章被收录于专栏：气python风雨气python风雨

运行总次数：0

代码可运行

WRF散度计算步骤及Python可视化实现2.0

项目概述

在气象家园论坛中，有用户对WRF模式输出的散度计算感到困惑。本项目旨在介绍WRF散度的计算方法，并提供Python可视化实现。在后续研究中发现不使用循环也可以完成计算，因此进行了测试比较。

环境设置

依赖库

wrf-python
metpy
numpy
xarray
matplotlib
cartopy

数据处理

数据读取

import numpy as np
import xarray as xr
import metpy.calc as mpcalc
from metpy.interpolate import cross_section
from wrf import getvar, interplevel, to_np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from netCDF4 import Dataset
from metpy.units import units
import metpy
# 读取WRF输出文件
file_path = "/home/mw/input/typhoon9537/wrfout_d01_2019-08-08_19_00_00"
ncfile = Dataset(file_path)

# 获取需要的变量
u = getvar(ncfile, "ua")  # 纬向风 (m/s)
v = getvar(ncfile, "va")  # 经向风 (m/s)
p = getvar(ncfile, "pressure")  # 气压 (hPa)
z = getvar(ncfile, "z")
lat = getvar(ncfile, "lat")
lon = getvar(ncfile, "lon")

计算方法比较

方法一：传统循环计算

%%timeit
# 计算网格间距
dx, dy = mpcalc.lat_lon_grid_deltas(lon.values, lat.values)

# 初始化散度数组
div = np.zeros_like(u.values)

# 逐层计算散度
for k in range(u.shape[0]):
    div[k, :, :] = mpcalc.divergence(
        u[k, :, :] * units('m/s'), 
        v[k, :, :] * units('m/s'),
        dx=dx, 
        dy=dy
    )

# 将散度数组转为xarray.DataArray
div_da = xr.DataArray(
    div,
    dims=u.dims,
    coords=u.coords,
    attrs={'units': 's^-1', 'description': 'divergence'}
)

<magic-timeit>:9: UserWarning: More than one time coordinate present for variable  "ua".
<magic-timeit>:9: UserWarning: Horizontal dimension numbers not found. Defaulting to (..., Y, X) order.


737 ms ± 15.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

方法二：扩展dx和dy维度

%%timeit
# 计算网格间距
dx, dy = mpcalc.lat_lon_grid_deltas(lon.values, lat.values)
# 扩展dx和dy的维度
dx = dx[None, :]
dy = dy[None, :]

# 直接计算散度
div = mpcalc.divergence(u, v, dx=dx, dy=dy)

<magic-timeit>:8: UserWarning: More than one time coordinate present for variable  "ua".
<magic-timeit>:8: UserWarning: Horizontal dimension numbers not found. Defaulting to (..., Y, X) order.


454 ms ± 7.29 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

方法三：使用grid_deltas_from_dataarray函数

%%timeit

# 直接从数据数组计算网格间距
dx, dy = metpy.xarray.grid_deltas_from_dataarray(u)

# 计算散度
div = mpcalc.divergence(u, v, dx=dx, dy=dy)

<magic-timeit>:2: UserWarning: More than one time coordinate present for variable  "XLAT".
<magic-timeit>:2: UserWarning: y and x dimensions unable to be identified. Assuming [..., y, x] dimension order.
<magic-timeit>:5: UserWarning: More than one time coordinate present for variable  "ua".
<magic-timeit>:5: UserWarning: Horizontal dimension numbers not found. Defaulting to (..., Y, X) order.


513 ms ± 4.88 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

方法比较与优化建议

计算效率：
- 传统循环方法计算速度较慢，特别是对于多层数据
- 方法二和方法三避免了循环，计算效率更高，方法二速度最快
代码简洁性：
- 方法三最为简洁，直接利用metpy的内置函数
- 方法二需要手动扩展维度
推荐方法：
- 对于新版本metpy，推荐使用方法三
- 对于旧版本，可以使用方法二
注意事项：
- 使用时会收到维度警告，可以安全忽略或通过指定维度消除
- 计算结果可能存在微小差异，主要源于浮点精度和算法实现

可视化示例

# 选择某一层进行可视化
level = 10# 选择第10层
div_slice = div.isel(bottom_top=level)

# 创建地图投影
proj = ccrs.PlateCarree()

# 创建图形
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection=proj)

# 绘制散度
contour = ax.contourf(
    lon.values, lat.values, div_slice.values*1e5,
    levels=np.linspace(-50, 50, 21),
    cmap='RdBu_r',
    transform=proj
)

# 添加地图要素
ax.coastlines()
ax.gridlines()
plt.colorbar(contour, ax=ax, label='Divergence (s$^{-1}$)')

# 添加标题
plt.title(f'WRF Divergence at Level {level}')

plt.show()