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

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

  • wrf-python概述
  • wrf-python安装
  • wrf-python主要函数
  • wrf-python和NCL总结
  • WRF后处理(未使用wrf-python库)
  • 批量添加.nc后缀
  • 提取单个变量:以降水为例
  • 提取所有变量
  • 计算气压(pressure)
  • 气压计算原理
  • WRF后处理-基于wrf-python
  • 参考
  • 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变量提取、插值、诊断量计算的函数,其中常用函数总结如下表:

    NCL WRF-python Value
    wrf_user_get_var wrf.getvar Extract and computes varaiables
    wrf_user_interp_level wrf.interplevel Interpolate a 3D field at the given vertical level(s)
    wrf_user_vert_cross wrf.vertcross Interpolates a 3D field through a user-specified horizontal line
    wrf_user_interp_line wrf.interpline Interpolates a 2D field to a user-specified line
    wrf_user_ll_to_xy wrf.ll_to_xy Lat/Lon <-> XY Routines

    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输出结果为例演示按照指定维度合并(附示例代码)

    作者:WW、forever

    物联沃分享整理
    物联沃-IOTWORD物联网 » 【WRF后处理】基于wrf-python处理wrf运行结果wrfout_d01

    发表回复