WRF | 如何使用wrfout计算锋生

用户11172986

发布于 2026-04-24 19:07:19

640

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

hello小伙伴们，这次带来的是wrf数据的锋生函数计算，依然是使用metpy函数计算

前期写过使用分析数据计算锋生，这次使用wrf数据进行计算与可视化

有什么有趣的提议可以评论区提出来

导入必要的库


from wrf import getvar, interplevel, latlon_coords, to_np
from netCDF4 import Dataset
import metpy.calc as mpcalc
from metpy.units import units
import numpy as np

数据处理

提取目标区域数据


# 打开 WRF 输出
ncfile = Dataset("/home/mw/input/typhoon9537/wrfout_d01_2019-08-08_19_00_00")  # 替换为你的文件名

# 获取位温（theta，单位 K）
theta = getvar(ncfile, "theta", timeidx=0)  # 单位：K

# 获取气压（P，单位 hPa）
p = getvar(ncfile, "pressure", timeidx=0)  # 单位：hPa

# 获取 U、V 风场（单位 m/s）
u = getvar(ncfile, "ua", timeidx=0)  # m/s
v = getvar(ncfile, "va", timeidx=0)  # m/s

# 插值到 1000 hPa
theta_800 = interplevel(theta, p, 800)  # 单位：K
u_800 = interplevel(u, p, 800)
v_800 = interplevel(v, p, 800)

# 获取经纬度
lat, lon = latlon_coords(theta_800)

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

# 计算锋生函数（单位：K / 100km / 3h）
frontogenesis = mpcalc.frontogenesis(
    theta_800 * units.kelvin,
    u_800 * units('m/s'),
    v_800 * units('m/s'),
    dx=dx, dy=dy
)

# 转换为常用单位：K / 100km / 3h
fronto_100km_3h = frontogenesis * 1e9 * (100 / 1000) * (3 * 3600)

/tmp/ipykernel_56/1626882428.py:26: UserWarning: More than one time coordinate present for variable  "theta_interp".
  frontogenesis = mpcalc.frontogenesis(
/tmp/ipykernel_56/1626882428.py:26: UserWarning: Horizontal dimension numbers not found. Defaulting to (..., Y, X) order.
  frontogenesis = mpcalc.frontogenesis(

数据可视化

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter


# 创建图形和坐标轴
fig = plt.figure(figsize=(12, 8))
ax = plt.axes(projection=ccrs.PlateCarree())

# 添加海岸线
ax.coastlines(resolution='10m', linewidth=0.8)

# 设置经纬度刻度和标签
ax.set_xticks(np.arange(-180, 181, 30), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(-90, 91, 30), crs=ccrs.PlateCarree())
ax.xaxis.set_major_formatter(LongitudeFormatter())
ax.yaxis.set_major_formatter(LatitudeFormatter())

# 绘制填色图
# 使用更好的颜色映射，避免使用'jet'
cs = ax.contourf(to_np(lon), to_np(lat), to_np(frontogenesis), 
                 levels=50, cmap='RdBu_r',  # 使用更科学的颜色映射
                 transform=ccrs.PlateCarree(),
                 extend='both')  # 扩展颜色条到数据范围之外

# # 添加等值线
# contour = ax.contour(to_np(lon), to_np(lat), to_np(frontogenesis),
#                      levels=10, colors='black', linewidths=0.5,
#                      transform=ccrs.PlateCarree())
# ax.clabel(contour, inline=True, fontsize=8, fmt='%.1f')

# 设置标题
plt.title("800 hPa 锋生函数 (K / 100km / 3h)", fontsize=14, pad=20)

# 添加颜色条
cbar = plt.colorbar(cs, ax=ax, orientation='horizontal', 
                   shrink=0.8, pad=0.1, aspect=40)
cbar.set_label('锋生强度 (K / 100km / 3h)', fontsize=12)

# 添加网格线
ax.gridlines(draw_labels=False, linewidth=0.5, color='gray', alpha=0.5)

# 设置图形布局
plt.tight_layout()

# 显示图形
plt.show()