xwrf | 如何基于xwrf绘制wrfout变量的剖面图

用户11172986

发布于 2026-04-24 18:55:36

800

xwrf | 如何基于xwrf绘制wrfout变量的剖面图

项目概述

在气象研究和数值模拟中，WRF（Weather Research and Forecasting）模式输出的数据（wrfout）通常包含多维气象变量（如温度、湿度、风场等）。为了分析这些变量在垂直方向上的分布特征，垂直剖面图（Cross-Section Plot）是一种常用的可视化方法。它能够直观地展示气象要素沿某一水平路径的垂直变化，有助于研究锋面、急流、对流系统等天气现象的结构特征。

Python中的 xwrf 库与metpy提供了便捷的WRF数据分析和可视化功能，可以高效地计算和绘制垂直剖面。本文将介绍如何使用 xwrf 提取 wrfout 数据，计算沿自定义路径的垂直剖面，并绘制符合学术出版标准的高质量图形。

环境设置

安装依赖

!pip install xwrf -i https://pypi.mirrors.ustc.edu.cn/simple/

Looking in indexes: https://pypi.mirrors.ustc.edu.cn/simple/
Collecting xwrf
  Downloading https://mirrors.ustc.edu.cn/pypi/packages/80/ff/304123d20207e272ccf05431352db8d828c10706809f1d64b764c77c4525/xwrf-0.0.4-py3-none-any.whl (19 kB)
Requirement already satisfied: xarray!=0.20.0,!=0.20.1,>=0.18 in /opt/conda/lib/python3.9/site-packages (from xwrf) (2024.2.0)
Requirement already satisfied: pyproj>=2.4.1 in /opt/conda/lib/python3.9/site-packages (from xwrf) (3.4.1)
Requirement already satisfied: donfig>=0.6.0 in /opt/conda/lib/python3.9/site-packages (from xwrf) (0.8.1.post0)
Requirement already satisfied: pyyaml in /opt/conda/lib/python3.9/site-packages (from donfig>=0.6.0->xwrf) (5.4.1)
Requirement already satisfied: certifi in /opt/conda/lib/python3.9/site-packages (from pyproj>=2.4.1->xwrf) (2024.2.2)
Requirement already satisfied: packaging>=22 in /opt/conda/lib/python3.9/site-packages (from xarray!=0.20.0,!=0.20.1,>=0.18->xwrf) (23.2)
Requirement already satisfied: numpy>=1.23 in /opt/conda/lib/python3.9/site-packages (from xarray!=0.20.0,!=0.20.1,>=0.18->xwrf) (1.26.4)
Requirement already satisfied: pandas>=1.5 in /opt/conda/lib/python3.9/site-packages (from xarray!=0.20.0,!=0.20.1,>=0.18->xwrf) (2.0.3)
Requirement already satisfied: python-dateutil>=2.8.2 in /opt/conda/lib/python3.9/site-packages (from pandas>=1.5->xarray!=0.20.0,!=0.20.1,>=0.18->xwrf) (2.8.2)
Requirement already satisfied: tzdata>=2022.1 in /opt/conda/lib/python3.9/site-packages (from pandas>=1.5->xarray!=0.20.0,!=0.20.1,>=0.18->xwrf) (2024.1)
Requirement already satisfied: pytz>=2020.1 in /opt/conda/lib/python3.9/site-packages (from pandas>=1.5->xarray!=0.20.0,!=0.20.1,>=0.18->xwrf) (2022.1)
Requirement already satisfied: six>=1.5 in /opt/conda/lib/python3.9/site-packages (from python-dateutil>=2.8.2->pandas>=1.5->xarray!=0.20.0,!=0.20.1,>=0.18->xwrf) (1.16.0)
Installing collected packages: xwrf
Successfully installed xwrf-0.0.4

导入库与读取

import xarray as xr
import xwrf
import matplotlib.pyplot as plt
import numpy as np
from metpy.interpolate import cross_section
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from matplotlib.colors import BoundaryNorm
# 1. 读取WRF输出文件并预处理
file_path = '/home/mw/input/typhoon9537/wrfout_d01_2019-08-08_23_00_00'  # 替换为你的文件路径
ds = xr.open_dataset(file_path).xwrf.postprocess()
print(ds)

<xarray.Dataset> Size: 617MB
Dimensions:                    (Time: 1, y: 437, x: 447, soil_layers_stag: 4,
                                z: 34, x_stag: 448, y_stag: 438, z_stag: 35,
                                seed_dim_stag: 2)
Coordinates: (12/15)
    XLAT                       (y, x) float32 781kB 20.89 20.89 ... 32.78 32.78
    XLONG                      (y, x) float32 781kB 116.5 116.6 ... 130.2 130.3
    XTIME                      (Time) datetime64[ns] 8B ...
    XLAT_U                     (y, x_stag) float32 783kB 20.89 20.89 ... 32.78
    XLONG_U                    (y, x_stag) float32 783kB 116.5 116.6 ... 130.3
    XLAT_V                     (y_stag, x) float32 783kB 20.88 20.88 ... 32.79
    ...                         ...
  * z_stag                     (z_stag) float32 140B 1.0 0.993 ... 0.013 0.0
  * Time                       (Time) datetime64[ns] 8B 2019-08-08T23:00:00
  * x_stag                     (x_stag) float64 4kB -6.705e+05 ... 6.705e+05
  * y                          (y) float64 3kB -6.54e+05 -6.51e+05 ... 6.54e+05
  * y_stag                     (y_stag) float64 4kB -6.555e+05 ... 6.555e+05
  * x                          (x) float64 4kB -6.69e+05 -6.66e+05 ... 6.69e+05
Dimensions without coordinates: soil_layers_stag, seed_dim_stag
Data variables: (12/160)
    Times                      (Time) |S19 19B b'2019-08-08_23:00:00'
    LU_INDEX                   (Time, y, x) float32 781kB ...
    ZS                         (Time, soil_layers_stag) float32 16B ...
    DZS                        (Time, soil_layers_stag) float32 16B ...
    VAR_SSO                    (Time, y, x) float32 781kB ...
    BATHYMETRY_FLAG            (Time) int32 4B ...
    ...                         ...
    geopotential_height        (Time, z_stag, y, x) float32 27MB 0.0 ... 2.09...
    wind_east                  (Time, z, y, x) float32 27MB 8.973 ... -8.632
    wind_north                 (Time, z, y, x) float32 27MB 0.8171 ... -1.929
    wind_east_10               (Time, y, x) float32 781kB 8.247 8.262 ... -1.364
    wind_north_10              (Time, y, x) float32 781kB 0.75 0.7638 ... 1.235
    wrf_projection             object 8B +proj=lcc +x_0=0 +y_0=0 +a=6370000 +...
Attributes: (12/134)
    TITLE:                            OUTPUT FROM WRF V4.4.2 MODEL
    START_DATE:                      2019-08-08_18:00:00
    SIMULATION_START_DATE:           2019-08-08_18:00:00
    WEST-EAST_GRID_DIMENSION:        448
    SOUTH-NORTH_GRID_DIMENSION:      438
    BOTTOM-TOP_GRID_DIMENSION:       35
    ...                              ...
    ISLAKE:                          21
    ISICE:                           15
    ISURBAN:                         13
    ISOILWATER:                      14
    HYBRID_OPT:                      2
    ETAC:                            0.2

剖面插值

# 2. 添加必要的metpy单位信息
ds = ds.metpy.parse_cf()
# 3. 定义剖面起点和终点 (经纬度坐标)
start_point = (25.0, 120.0)  # (纬度, 经度)
end_point = (30.0, 125.0)

# 4. 创建垂直剖面
cross = cross_section(
    ds['air_potential_temperature'],  # 温度场
    start_point,
    end_point,
    steps=100,  # 剖面点数
    interp_type='linear'
)
print(cross)

<xarray.DataArray 'air_potential_temperature' (Time: 1, z: 34, index: 100)> Size: 27kB
<Quantity([[[302.7520626  302.79776204 302.84359141 ... 301.93200097 302.06725171
   302.18286203]
  [302.88180945 302.93319434 302.96634203 ... 301.93839125 302.07870443
   302.2031821 ]
  [303.00872577 303.03460497 303.06402959 ... 301.99097732 302.12458208
   302.24540202]
  ...
  [401.08947976 399.95359748 398.63346022 ... 400.49328285 400.08702435
   399.62845682]
  [431.16560796 432.79708979 433.78063201 ... 424.39098679 423.81428777
   423.71475879]
  [470.19364105 470.54106535 470.01827391 ... 466.17964465 465.91242329
   465.96753058]]], 'kelvin')>
Coordinates:
    XTIME      (Time) datetime64[ns] 8B 2019-08-08T23:00:00
  * Time       (Time) datetime64[ns] 8B 2019-08-08T23:00:00
    metpy_crs  object 8B Projection: lambert_conformal_conic
    XLAT       (index) float64 800B 25.0 25.05 25.1 25.15 ... 29.9 29.95 30.0
    XLONG      (index) float64 800B 120.0 120.0 120.1 ... 124.9 124.9 125.0
    CLAT       (index) float64 800B 25.0 25.05 25.1 25.15 ... 29.9 29.95 30.0
  * z          (z) float32 136B 0.9965 0.988 0.9765 ... 0.033 0.0195 0.0065
    x          (index) float64 800B -2.988e+05 -2.939e+05 ... 1.896e+05
    y          (index) float64 800B -2.156e+05 -2.101e+05 ... 3.303e+05
  * index      (index) int64 800B 0 1 2 3 4 5 6 7 8 ... 92 93 94 95 96 97 98 99
Attributes:
    standard_name:  air_potential_temperature
    grid_mapping:   wrf_projection

简单绘图

import cmaps
# 5. 准备绘图数据
# 获取气压坐标并转换为hPa
pressure = cross['z']*1000
cross_section_dist = np.linspace(0, 1, len(cross['index']))  # 标准化距离

# 创建等值线层级
levels =np.arange(300,380,5)

# 6. 创建图形
fig, ax = plt.subplots(figsize=(12, 8))

# 绘制填色等值线
cf = ax.contourf(
    cross_section_dist,
    pressure,
    cross[0],
    levels=levels,
    cmap='YlGnBu'
)

# 添加等高线
cs = ax.contour(
    cross_section_dist,
    pressure,
    cross[0],
    levels=levels,
    colors='k',
    linewidths=0.5
)
ax.clabel(cs, fontsize=8, fmt='%d')

# 设置坐标轴
ax.set_ylabel('Pressure (hPa)')
ax.set_xlabel('Cross-Section Distance')
ax.set_ylim(1000, 100)  # 地面到100hPa
# ax.invert_yaxis()  # 气压反向

# 添加色标
cbar = fig.colorbar(cf, ax=ax, orientation='vertical')
cbar.set_label('Temperature (K)')

# 添加标题
ax.set_title(
    f' Cross-Section\n'
    f'Start: {start_point}°N/°E → End: {end_point}°N/°E\n'
    f'Time: {ds.Time.dt.strftime("%Y-%m-%d %H:%M").values}'
)

# 7. 显示图形
plt.tight_layout()
plt.show()