当前位置： 首页 > article >正文

【WRF后处理】提取某要素数据并绘制地图

article 2024/11/15 13:26:26

【WRF后处理】提取某要素数据并绘制地图

    • 根据Domain提取文件
    • 提取某要素数据并绘制地图
  • 参考

根据Domain提取文件

为了满足根据需求提取不同区域（例如 D01、D02、D03）的文件，可以将区域作为函数的一个参数传入，并根据传入的区域自动匹配相应的文件。使用正则表达式来动态生成匹配模式，从而提取不同区域的 wrfout 文件。

Python代码如下：

import os
import re

def get_wrfout_files_by_domain(directory, domain):
    """
    从指定目录中提取指定域 (D01, D02, D03) 的 wrfout 文件。
    
    参数:
    directory (str): 文件夹路径
    domain (str): 需要提取的域 (例如 "d01", "d02", "d03")
    
    返回:
    list: 包含所有符合条件的 wrfout 文件的绝对路径
    """
    # 检查 domain 参数是否正确
    if domain.lower() not in ["d01", "d02", "d03"]:
        raise ValueError("域必须是 'd01'、'd02' 或 'd03'")
    
    # 用于存储符合条件的文件
    wrfout_files = []
    
    # 定义 wrfout 文件的正则表达式模式，动态匹配域
    pattern = re.compile(fr"^wrfout_{domain}_\d{{4}}-\d{{2}}-\d{{2}}_\d{{2}}_\d{{2}}_\d{{2}}$")
    
    # 遍历目录中的所有文件
    for filename in os.listdir(directory):
        # 检查文件名是否匹配 wrfout_{domain} 的格式
        if pattern.match(filename):
            # 将符合条件的文件的绝对路径添加到列表中
            wrfout_files.append(os.path.join(directory, filename))
    
    return wrfout_files

# 示例用法
directory = "/path/to/your/wrfout/files"  # 替换为实际的文件夹路径

# 提取 D03 区域文件
wrfout_d03_files = get_wrfout_files_by_domain(directory, "d03")
for file in wrfout_d03_files:
    print("D03区域文件:", file)

# 提取 D01 区域文件
wrfout_d01_files = get_wrfout_files_by_domain(directory, "d01")
for file in wrfout_d01_files:
    print("D01区域文件:", file)

# 提取 D02 区域文件
wrfout_d02_files = get_wrfout_files_by_domain(directory, "d02")
for file in wrfout_d02_files:
    print("D02区域文件:", file)

提取某要素数据并绘制地图

在这里插入图片描述
Python代码如下：

import os
import re
from netCDF4 import Dataset
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from wrf import getvar, to_np, get_cartopy, latlon_coords
import geopandas as gpd
from shapely.geometry import Polygon

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

def get_wrfout_files_by_domain(directory, domain):
    """
    从指定目录中提取指定域 (D01, D02, D03) 的 wrfout 文件。

    参数:
    directory (str): 文件夹路径
    domain (str): 需要提取的域 (例如 "d01", "d02", "d03")

    返回:
    list: 包含所有符合条件的 wrfout 文件的绝对路径
    """
    # 检查 domain 参数是否正确
    if domain.lower() not in ["d01", "d02", "d03"]:
        raise ValueError("域必须是 'd01'、'd02' 或 'd03'")

    # 用于存储符合条件的文件
    wrfout_files = []

    # 定义 wrfout 文件的正则表达式模式，动态匹配域
    pattern = re.compile(fr"^wrfout_{domain}_\d{{4}}-\d{{2}}-\d{{2}}_\d{{2}}_\d{{2}}_\d{{2}}$")

    # 遍历目录中的所有文件
    for filename in os.listdir(directory):
        # 检查文件名是否匹配 wrfout_{domain} 的格式
        if pattern.match(filename):
            # 将符合条件的文件的绝对路径添加到列表中
            wrfout_files.append(os.path.join(directory, filename))

    return wrfout_files


def list_variables_in_wrfout(file_path):
    """
    列出给定 wrfout 文件中的所有 WRF 变量。

    参数:
    file_path (str): wrfout 文件的路径
    """
    # 使用 netCDF4 打开 wrfout 文件
    dataset = Dataset(file_path)

    print(f"\n文件: {file_path}")
    print("包含的变量:")

    # 使用 wrf-python 的 getvar 函数获取所有可用变量
    for varname in dataset.variables:
        print(varname)

    # 关闭 NetCDF 文件
    dataset.close()


def add_shapefile_to_ax(ax, shapefile_path):
    """
    将 shapefile 添加到指定的 cartopy 绘图轴上。

    参数:
    ax (geoaxes.GeoAxesSubplot): cartopy 投影的绘图轴
    shapefile_path (str): shapefile 的路径
    """
    # 使用 geopandas 读取 shapefile
    gdf = gpd.read_file(shapefile_path)

    # 检查 shapefile 是否有坐标系
    if gdf.crs is not None:
        print(f"Shapefile {shapefile_path} 的原始坐标系: {gdf.crs}")
    else:
        raise ValueError(f"Shapefile {shapefile_path} 没有坐标系！")

    # 将 shapefile 转换为 PlateCarree 投影，这是 cartopy 常用的地理坐标系
    gdf = gdf.to_crs(ccrs.PlateCarree().proj4_init)

    # 将 GeoDataFrame 中的几何体逐一添加到 ax 中
    for geom in gdf.geometry:
        ax.add_geometries([geom], crs=ccrs.PlateCarree(), edgecolor='black', facecolor='none', linewidth=1)


def extract_and_plot_t2(file_path, output_dir):
    """
    从 wrfout 文件中提取 2米温度 (T2) 并绘制地图。

    参数:
    file_path (str): wrfout 文件的路径
    output_dir (str): 保存地图的目录路径
    """
    # 打开 wrfout 文件
    dataset = Dataset(file_path)

    # 使用 wrf-python 提取 2米温度 (T2)
    t2 = getvar(dataset, "T2")

    # 将温度从开尔文转换为摄氏度
    t2_celsius = t2 - 273.15

    # 提取经纬度坐标
    lats, lons = latlon_coords(t2)

    # 获取 cartopy 投影
    cart_proj = get_cartopy(t2)

    # 创建绘图
    fig = plt.figure(figsize=(12, 10))
    ax = plt.axes(projection=cart_proj)

    # 绘制T2的值（即2米温度，单位为摄氏度）
    t2_contour = ax.contourf(to_np(lons), to_np(lats), to_np(t2_celsius), 10, transform=ccrs.PlateCarree(), cmap="coolwarm")

    # 添加海岸线和州界
    # ax.coastlines('50m', linewidth=0.8)
    # ax.add_feature(cfeature.BORDERS, linewidth=0.5)

    # 加载并绘制 SHP 文件（大湾区和中国边界）
    GBAshp_file_path = "C:/PycharmProjects/GBA_Data_Process/shpFile/GBA/GBA.shp"
    Chinashp_file_path = "C:/PycharmProjects/GBA_Data_Process/shpFile/China/中国.shp"

    # 使用转换后的绘制方法
    add_shapefile_to_ax(ax, GBAshp_file_path)
    add_shapefile_to_ax(ax, Chinashp_file_path)

    # 添加经纬度网格和标签
    gl = ax.gridlines(draw_labels=True, linewidth=0.5, color='gray', alpha=0.5, linestyle='--')
    gl.top_labels = False  # 关闭顶部标签
    gl.right_labels = False  # 关闭右侧标签
    gl.bottom_labels = False  # 关闭底部标签（稍后手动添加）

    # 设置经纬度刻度的字体大小
    gl.xlabel_style = {'size': 14, 'color': 'black'}
    gl.ylabel_style = {'size': 14, 'color': 'black'}

    # 手动添加经度标签，确保它们显示在底部外侧
    lon_labels = np.arange(-180, 181, 1)  # 经度范围 [-180, 180]，间隔 1 度
    for lon in lon_labels:
        ax.text(lon, -89.5, f"{lon}°", va='top', ha='center', fontsize=12,
                transform=ccrs.PlateCarree())  # 这里使用 -89.5 代替 -90

    # 添加自定义的 X 和 Y 轴标签（经度和纬度的名称）
    ax.text(0.5, -0.02, "Longitude (°E)", va='top', ha='center', fontsize=16, transform=ax.transAxes)
    ax.text(-0.10, 0.5, "Latitude (°N)", va='center', ha='right', fontsize=16, rotation='vertical', transform=ax.transAxes)

    # 设置经度和纬度的刻度在外侧显示
    ax.xaxis.set_ticks_position('bottom')  # 经度刻度放在下侧
    ax.xaxis.set_label_position('bottom')  # 经度标签放在下侧
    ax.yaxis.set_ticks_position('left')    # 纬度刻度放在左侧
    ax.yaxis.set_label_position('left')    # 纬度标签放在左侧

    # 将刻度设置为指向图形外部
    ax.tick_params(axis='both', which='both', direction='out')

    # 添加颜色栏（竖直放置，放置在图形右侧）
    cbar = plt.colorbar(t2_contour, ax=ax, orientation="vertical", pad=0.05, aspect=20, shrink=0.96)
    cbar.set_label("Temperature (°C)", fontsize=14)

    # 修改 colorbar 刻度的字体大小
    cbar.ax.tick_params(labelsize=12)  # 将字体大小设为 12

    # 设置标题
    plt.title(f"2-meter Temperature (T2) in Celsius from {os.path.basename(file_path)}", fontsize=14)

    # 保存图像
    output_file = os.path.join(output_dir, f"T2_map_{os.path.basename(file_path)}.png")
    plt.savefig(output_file, dpi=600)
    plt.close()

    print(f"保存地图: {output_file}")

    # 关闭 NetCDF 文件
    dataset.close()

def extract_and_plot_precipitation(file_path, output_dir):
    """
    从 wrfout 文件中提取累计降水 (RAINC + RAINNC) 并绘制地图。

    参数:
    file_path (str): wrfout 文件的路径
    output_dir (str): 保存地图的目录路径
    """
    # 打开 wrfout 文件
    dataset = Dataset(file_path)

    # 使用 wrf-python 提取累计对流降水 (RAINC) 和累计非对流降水 (RAINNC)
    rainc = getvar(dataset, "RAINC")  # 累计对流降水
    rainnc = getvar(dataset, "RAINNC")  # 累计非对流降水

    # 累计降水量 = RAINC + RAINNC
    total_precip = rainc + rainnc

    # 提取经纬度坐标
    lats, lons = latlon_coords(total_precip)

    # 获取 cartopy 投影 - 改为使用 T2 变量或其他原始变量来获取投影信息
    t2 = getvar(dataset, "T2")  # 使用 T2 提取投影
    cart_proj = get_cartopy(t2)

    # 创建绘图
    fig = plt.figure(figsize=(12, 10))
    ax = plt.axes(projection=cart_proj)

    # 设置图形范围与 WRF 数据一致
    ax.set_extent([np.min(to_np(lons)), np.max(to_np(lons)), np.min(to_np(lats)), np.max(to_np(lats))], crs=ccrs.PlateCarree())

    # 绘制累计降水量（单位：毫米）
    precip_contour = ax.contourf(to_np(lons), to_np(lats), to_np(total_precip), 10, transform=ccrs.PlateCarree(), cmap="Blues")

    # 添加海岸线和州界
    # ax.coastlines('50m', linewidth=0.8)
    # ax.add_feature(cfeature.BORDERS, linewidth=0.5)

    # 加载并绘制 SHP 文件（大湾区和中国边界）
    GBAshp_file_path = "C:/PycharmProjects/GBA_Data_Process/shpFile/GBA/GBA.shp"
    Chinashp_file_path = "C:/PycharmProjects/GBA_Data_Process/shpFile/China/中国.shp"

    # 使用转换后的绘制方法
    add_shapefile_to_ax(ax, GBAshp_file_path)
    add_shapefile_to_ax(ax, Chinashp_file_path)

    # 添加经纬度网格和标签
    gl = ax.gridlines(draw_labels=True, linewidth=0.5, color='gray', alpha=0.5, linestyle='--')
    gl.top_labels = False  # 关闭顶部标签
    gl.right_labels = False  # 关闭右侧标签
    gl.bottom_labels = False  # 关闭底部标签（稍后手动添加）

    # 设置经纬度刻度的字体大小
    gl.xlabel_style = {'size': 14, 'color': 'black'}
    gl.ylabel_style = {'size': 14, 'color': 'black'}

    # 手动添加经度标签，确保它们显示在底部外侧
    lon_labels = np.arange(-180, 181, 1)  # 经度范围 [-180, 180]，间隔 1 度
    for lon in lon_labels:
        ax.text(lon, -89.5, f"{lon}°", va='top', ha='center', fontsize=12,
                transform=ccrs.PlateCarree())  # 这里使用 -89.5 代替 -90

    # 添加自定义的 X 和 Y 轴标签（经度和纬度的名称）
    ax.text(0.5, -0.02, "Longitude (°E)", va='top', ha='center', fontsize=16, transform=ax.transAxes)
    ax.text(-0.10, 0.5, "Latitude (°N)", va='center', ha='right', fontsize=16, rotation='vertical', transform=ax.transAxes)

    # 设置经度和纬度的刻度在外侧显示
    ax.xaxis.set_ticks_position('bottom')  # 经度刻度放在下侧
    ax.xaxis.set_label_position('bottom')  # 经度标签放在下侧
    ax.yaxis.set_ticks_position('left')    # 纬度刻度放在左侧
    ax.yaxis.set_label_position('left')    # 纬度标签放在左侧

    # 将刻度设置为指向图形外部
    ax.tick_params(axis='both', which='both', direction='out')

    # 添加颜色栏（竖直放置，放置在图形右侧）
    cbar = plt.colorbar(precip_contour, ax=ax, orientation="vertical", pad=0.05, aspect=20, shrink=0.96)
    cbar.set_label("Cumulative Precipitation (mm)", fontsize=14)

    # 修改 colorbar 刻度的字体大小
    cbar.ax.tick_params(labelsize=12)  # 将字体大小设为 12

    # 设置标题
    plt.title(f"Cumulative Precipitation from {os.path.basename(file_path)}", fontsize=14)

    # 保存图像
    output_file = os.path.join(output_dir, f"Precipitation_map_{os.path.basename(file_path)}.png")
    plt.savefig(output_file, dpi=600)
    plt.close()

    print(f"保存降水地图: {output_file}")

    # 关闭 NetCDF 文件
    dataset.close()


# 示例用法
directory = "C:/Database/wrfout"                   # 替换为实际的文件夹路径
output_dir = "C:/Database/wrfout-Outputs/Figures"  # 替换为保存地图的路径

# 检查输出目录是否存在，如果不存在则创建
if not os.path.exists(output_dir):
    os.makedirs(output_dir)

# 提取 D03 区域文件
wrfout_d03_files = get_wrfout_files_by_domain(directory, "d03")
for file in wrfout_d03_files:
    print("D03区域文件:", file)

    # 遍历 D03 文件，列出每个文件的变量
    # list_variables_in_wrfout(file)

    # 提取T2并绘制地图
    # extract_and_plot_t2(file, output_dir)

    # 提取降水并绘制地图
    extract_and_plot_precipitation(file, output_dir)

参考


http://www.kler.cn/a/394443.html

相关文章：

  • UAC2.0 speaker——同时支持 16bit，24bit 和 32bit
  • Spring Cloud Gateway（分发请求）
  • Vue2：组件
  • 华为大变革？仓颉编程语言会代替ArkTS吗？
  • java八股-jvm入门-程序计数器，堆，元空间，虚拟机栈，本地方法栈，类加载器，双亲委派，类加载执行过程
  • C++模板特化实战：在使用开源库boost::geometry::index::rtree时，用特化来让其支持自己的数据类型
  • 基于Java Springboot剧本杀管理系统
  • webSocket的使用文档
  • MySQL【四】
  • 【.GetConnectionTimeoutException的2种情况分析】
  • 打包python代码为exe文件
  • Flutter：Widget生命周期
  • Spring MVC进阶
  • R语言基础| 机器学习
  • 改扩配系列：浪潮英政服务器CS5280H2、IR5280H2——后置SATA、NVME硬盘安装
  • SpringBoot实战：AI大模型+亮数据代理高效获取视频资源
  • 【Apache Paimon】-- 1 -- Apache Paimon 是什么？
  • Python Pandas 结构之 Series 和 DataFrame
  • NFS存储基础操作
  • PostgreSQL 行转列实现
  • 存储大挑战：如何在可靠性与大容量之间玩转平衡术？
  • LabVIEW 使用 Snippet
  • 【Excel】数据透视表分析方法大全
  • 【C++进阶实战】基于linux的天气预报系统
  • CTF攻防世界小白刷题自学笔记15
  • 【Golang】golang框架，为什么选择GoFrame, GoFrame使用心得