【WRF后处理】基于wrf-python处理wrf运行结果wrfout_d01

WRF的模拟结果是按照指定的时间间隔和模拟时间段依次输出的。

wrf-python概述

wrf-python 是一个 Python 库包，专门用于处理和分析 WRF（Weather Research and Forecasting）模型输出的数据。它提供了一系列函数，帮助用户更轻松地从 WRF 输出的 NetCDF 文件中提取和计算常见的气象变量，并进行相关的气象分析和绘图。

wrf-python 提供了许多方便的工具来处理 WRF 模型的输出，例如：
1、提取和计算常见气象变量：
例如，温度、湿度、风速、风向、相对湿度、位势温度、气压等。
它封装了很多 WRF 模型中的复杂计算，用户可以直接使用内置的函数计算这些变量，而不必手动处理复杂的数学公式。
2、插值和网格处理：
提供了水平和垂直插值的功能，例如将数据从 eta 坐标插值到等压面。
支持将 WRF 的 Arakawa-C 网格数据插值到常规的网格上，便于可视化。
3、地形和地图投影支持：
wrf-python 支持处理 WRF 的地图投影、地形高度等信息，便于与地图和其他地理数据结合使用。
4、集成 Numpy 和 Xarray：
wrf-python 与 xarray、numpy 等库集成良好，能够方便地处理大规模的气象数据，并进行并行计算和分析。
5、可视化支持：
wrf-python 可以与 matplotlib 等可视化库配合使用，方便绘制等高线图、风场图、气压图等各种气象图。

wrf-python安装

1、可以通过 Python 包管理工具 pip 安装 wrf-python。安装命令如下：

pip install wrf-python

2、如果你使用的是 conda 环境，也可以通过 conda-forge 安装：

conda install wrf-python
conda install -c conda-forge wrf-python

3、检查是否安装成功，命令如下：

conda list wrf-python

wrf-python主要函数

wrf-python 提供了许多常见的函数，以下列举了一些常用的：

1、getvar()：从 WRF 文件中提取常见的气象变量，例如温度、气压、风速等。

temp = getvar(wrf_file, "T2")  # 提取 2 米温度

2、interplevel()：将数据从模型层插值到常见的大气层面（例如等压面）。

temp_500 = interplevel(temp, pressure, 500)  # 插值到 500 hPa 等压面

3、latlon_coords()：提取 WRF 网格的纬度和经度坐标。

4、get_basemap()：获取与 WRF 文件相关联的 Basemap 地图投影对象，便于绘图。

wrf-python和NCL总结

wrf-python和NCL都有对应的wrf变量提取、插值、诊断量计算的函数，其中常用函数总结如下表：

WRF后处理（未使用wrf-python库）

批量添加.nc后缀

因为一般WRF 默认输出文件的文件名后缀没有.nc，无法直接使用xarray进行读取，也就用不了concat函数。所以首先需要给所有的输出文件批量添加后缀名".nc"。

参考Python代码如下：

#导入库
import numpy as np
import xarray as xr
import os
from netCDF4 import Dataset

# 设置 WRF 模型输出文件的路径
# 如果 WRF_output_path 为空，则使用当前工作目录
WRF_output_path = ""
if WRF_output_path == "":
    path = os.getcwd()  # 获取当前工作目录
else:
    path = WRF_output_path  # 使用指定的路径

# 批量修改文件名：为文件加上 .nc 后缀，方便使用 xarray 处理
file_names = os.listdir(path)  # 列出路径中的所有文件
for file in file_names:
    # 仅处理以 "wrfout_d01" 开头的文件，排除其他无关文件
    if file.startswith('wrfout_d01'):
        base_name = os.path.basename(file)  # 获取文件名
        # 如果文件名已经有 .nc 后缀，跳过重命名
        if not base_name.endswith('.nc'):
            new_name = base_name + '.nc'  # 添加 .nc 后缀
            os.rename(os.path.join(path, file), os.path.join(path, new_name))  # 执行重命名

# 选取路径下所有前缀名为 'wrfout_d01' 的 .nc 文件
list_names = []
for ncfile in os.listdir(path):
    # 仅选择前缀为 'wrfout_d01' 且后缀为 '.nc' 的文件
    if ncfile.startswith('wrfout_d01') and ncfile.endswith('.nc'):
        list_names.append(ncfile)

# 将文件名按照时间进行排序
# 一般 WRF 文件名格式为 wrfout_d01_YYYY-MM-DD_HH:MM:SS
list_names_sort = np.sort(list_names)

# 打印排序后的文件名列表，便于检查
print("排序后的文件名列表：")
for name in list_names_sort:
    print(name)

提取单个变量：以降水为例

以单个变量P为例（可按需更改），基于concat函数，按照时间顺序进行合并。

参考Python代码如下：

# 导入必要的库
import numpy as np
import xarray as xr
import os  

# 假设已经有排序后的文件列表 list_names_sort

# 创建一个空的列表用于存储变量 P
file_list = []

# 当前目录路径（或你自定义的 WRF 输出路径）
path = os.getcwd()

# 遍历排序后的文件名，读取每个文件中的变量 P
for i in list_names_sort:
    file_path = os.path.join(path, i)  # 构建文件的完整路径
    print(f"正在读取文件: {file_path}")  # 打印当前读取的文件
    ds = xr.open_dataset(file_path)  # 打开 NetCDF 文件
    file_list.append(ds['P'])  # 提取变量 'P' 并添加到列表中

# 按时间维度合并所有文件中的 P 变量
# 如果数据集中没有明确的时间维度，可以在 concat 时显式指定
data = xr.concat(file_list, dim="Time")

# 保存合并后的数据为新的 NetCDF 文件
output_file = 'wrf_data_P.nc'  # 输出文件名
print(f"保存合并后的数据到 {output_file}")
data.to_netcdf(output_file)

# 提示完成
print("数据合并并保存完成！")

提取所有变量

参考Python代码如下：

# 导入必要的库
import numpy as np
import xarray as xr
import os  

# 假设已经有排序后的文件列表 list_names_sort

# 使用 xarray 的 open_mfdataset 进行多文件合并，按时间维度拼接
# 使用 combine='nested' 并设置 concat_dim='Time' 按时间维度合并
data = xr.open_mfdataset([os.path.join(os.getcwd(), f) for f in list_names_sort], 
                         concat_dim="Time", combine="nested", engine="netcdf4")

# 保存合并后的所有变量为新的 NetCDF 文件
output_file = 'wrf_all_variables.nc'
print(f"保存合并后的数据到 {output_file}")
data.to_netcdf(output_file)

# 提示完成
print("所有变量合并并保存完成！")

计算气压（pressure）

气压计算原理

在 WRF 模型中，气压（pressure）是通过将微扰气压（P）和基础状态气压（PB）相加来得到的：

pressure = P + PB

其中：

• P：微扰气压，单位为帕斯卡（Pa）。这是气压相对于基础状态的偏差，代表了由于天气系统、地形影响等因素造成的气压波动。这个变量是 WRF 模型中求解的主要变量之一，描述了大气中相对于基础状态的动态变化。
• PB：基础状态气压，单位为帕斯卡（Pa）。这是一个理想化的大气静力平衡气压，通常是由模型的静力方程给出的。它代表了一个平衡状态下的大气压力分布，通常是一个水平分布的静态气压。
• 合成气压（P + PB）
  通过将两者相加，得到的是完整的、真实的大气气压（pressure），即模型中的总气压。这个总气压表示了在模型每个网格点处的大气压力，包括了基础静力气压和由于天气系统等动态引起的气压变化。

参考Python代码如下：

# 导入必要的库
import numpy as np
import xarray as xr
import os  

# 假设已经有排序后的文件列表 list_names_sort

# 使用 xarray 的 open_mfdataset 进行多文件合并，按时间维度拼接
# 使用 combine='nested' 并设置 concat_dim='Time' 按时间维度合并
data = xr.open_mfdataset([os.path.join(os.getcwd(), f) for f in list_names_sort], 
                         concat_dim="Time", combine="nested", engine="netcdf4")

# 计算气压：P + PB
# P 是微扰气压，PB 是基础状态气压，计算结果存储在新的变量 'pressure' 中
data['pressure'] = data['P'] + data['PB']

# 为新变量 'pressure' 添加元数据（属性）
data['pressure'].attrs['FieldType'] = '104'
data['pressure'].attrs['MemoryOrder'] = 'XYZ'
data['pressure'].attrs['description'] = 'Full Model Pressure'
data['pressure'].attrs['units'] = 'Pa'
data['pressure'].attrs['stagger'] = ' '  # 无 stagger

# 将合并后的数据保存为新的 NetCDF 文件
output_file = 'wrf_data_with_pressure.nc'
print(f"保存合并后的数据到 {output_file}")
data.to_netcdf(output_file)

# 提示完成
print("所有变量合并并保存完成，气压计算已添加！")

WRF后处理-基于wrf-python

提取气压数据并绘制图像，图形如下：

相应Python代码如下：

import netCDF4 as nc
from wrf import getvar, interplevel, to_np, get_basemap, latlon_coords
import matplotlib.pyplot as plt
import geopandas as gpd

# 设置字体为 Times New Roman
plt.rcParams['font.family'] = 'Times New Roman'

# 打开 WRF 模型的 NetCDF 文件
WRF_output_path = "C:/Database/wrfout/wrfout_d01_2020-07-06_00_00_00"
wrf_file = nc.Dataset(WRF_output_path)

# 提取地面气压（注意 getvar 自动处理 P + PB）
pressure = getvar(wrf_file, "pressure")

# 提取经纬度
lats, lons = latlon_coords(pressure)

# 检查数据维度
print(f"Pressure data dimensions: {pressure.shape}")

# 假设 pressure 是 3D（Time, south_north, west_east），我们选择第一个时间步
# 选择第一个时间步
pressure_2d = pressure[0, :, :]  # 第一个时间步

# 打印选择后的气压数据形状
print(f"2D Pressure data shape: {pressure_2d.shape}")

# 提取纬度和经度的范围
min_lat = to_np(lats).min()
max_lat = to_np(lats).max()
min_lon = to_np(lons).min()
max_lon = to_np(lons).max()

# 打印经纬度的范围，便于检查
print(f"Latitude range: {min_lat} to {max_lat}")
print(f"Longitude range: {min_lon} to {max_lon}")

# 绘制等高线图，将经纬度作为 x 和 y 轴，确保地理位置正确
plt.figure(figsize=(10, 8))
contourf = plt.contourf(to_np(lons), to_np(lats), to_np(pressure_2d), cmap='coolwarm', levels=10)

# 设置经纬度范围与 WRF 模型一致
plt.xlim([min_lon, max_lon])
plt.ylim([min_lat, max_lat])

# 加载并绘制 SHP 文件
GBAshp_file_path = "C:/PycharmProjects/GBA_Data_Process/shpFile/GBA/GBA.shp"
Chinashp_file_path = "C:/PycharmProjects/GBA_Data_Process/shpFile/China/中国.shp"
gdf_GBA = gpd.read_file(GBAshp_file_path)
gdf_China = gpd.read_file(Chinashp_file_path)

# 绘制 SHP 文件，将边界绘制在图上
gdf_GBA.boundary.plot(ax=plt.gca(), color='black', linewidth=0.5)
gdf_China.boundary.plot(ax=plt.gca(), color='black', linewidth=1)

# 添加颜色条
cbar = plt.colorbar(contourf)

# 设置颜色条标签和字体大小
cbar.set_label("Pressure (Pa)", fontsize=14)
cbar.ax.tick_params(labelsize=12)

# 设置坐标轴标签
plt.xlabel('Longitude (°)', fontsize=16)
plt.ylabel('Latitude (°)', fontsize=16)
plt.title("Surface Pressure", fontsize=16)

# 设置坐标轴刻度值字体大小
plt.tick_params(axis='both', labelsize=14)

# 自定义坐标轴刻度值，格式化为带小数点和度数符号
def format_ticks(x, pos):
    return f'{x:.1f}°'

# 应用自定义格式化函数
plt.gca().xaxis.set_major_formatter(plt.FuncFormatter(format_ticks))
plt.gca().yaxis.set_major_formatter(plt.FuncFormatter(format_ticks))

# 导出图形
output_file_path = "C:/PycharmProjects/20241111 WRF_Postprocess/Figures/Surface Pressure.png"  # 修改为所需的输出路径
plt.savefig(output_file_path, dpi=300, bbox_inches='tight')

# 显示图像
plt.show()

参考

1、CSDN博客-WRF后处理总结：wrf-python与NCL在WRF后处理中的基本应用——变量提取、计算与可视化
2、用Python批处理指定数据-以WRF输出结果为例演示按照指定维度合并(附示例代码)