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

admin管理员组

文章数量:1438520

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

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

项目概述

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

环境设置

依赖库

• wrf-python
• metpy
• numpy
• xarray
• matplotlib
• cartopy

数据处理

数据读取

代码语言：javascript代码运行次数：0运行复制

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")

计算方法比较

方法一：传统循环计算

代码语言：javascript代码运行次数：0运行复制

%%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'}
)

代码语言：javascript代码运行次数：0运行复制

<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维度

代码语言：javascript代码运行次数：0运行复制

%%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)

代码语言：javascript代码运行次数：0运行复制

<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函数

代码语言：javascript代码运行次数：0运行复制

%%timeit

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

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

代码语言：javascript代码运行次数：0运行复制

<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)

方法比较与优化建议

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

可视化示例

代码语言：javascript代码运行次数：0运行复制

# 选择某一层进行可视化
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()

细化绘图

总结

1. WRF散度计算可以通过多种方法实现，推荐使用metpy的内置函数
2. 避免循环可以显著提高计算效率，特别是处理多层数据时
3. 计算结果可以通过xarray和matplotlib方便地进行可视化

参考资料

1. 气象家园论坛讨论
2. Stack Overflow上的相关讨论
3. MetPy官方文档

作者信息

• 公众号: 气python风雨

本文参与 腾讯云自媒体同步曝光计划，分享自微信公众号。原始发表：2025-04-21，如有侵权请联系 cloudcommunity@tencent 删除数据数组python函数可视化

本文标签： WRF散度计算步骤及Python可视化实现20