ERA5水汽通量散度剖面计算与绘图

用户11172986

发布于 2025-01-14 20:53:09

13000

代码可运行

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

运行总次数：0

代码可运行

ERA5水汽通量散度剖面计算与绘图

温馨提示

由于可视化代码过长隐藏，可点击运行Fork查看若没有成功加载可视化图，点击运行可以查看 ps:隐藏代码在【代码已被隐藏】所在行，点击所在行，可以看到该行的最右角，会出现个三角形，点击查看即可

前言

之前的文章中，有朋友提出水汽通量散度剖面图怎么画，那么我们来探索一下

项目目标

本项目旨在通过 Python 编程语言，结合气象数据处理库（如 xarray、metpy）和可视化工具（如 matplotlib、cartopy），实现以下目标：

计算整层水汽通量散度：基于气象数据（如 ERA5 再分析数据），计算从地表到特定高度范围内的水汽通量散度。

绘制水汽通量散度剖面图：通过剖面图直观展示水汽通量散度在垂直方向上的分布特征。

优化可视化效果：通过调整坐标轴、颜色条、地图投影等参数，提升剖面图的可读性和美观性。

嵌入小地图：在剖面图中嵌入小地图，显示剖面路径和地理信息，增强图的实用性。

项目功能

计算数据存为nc

以下是实现整层水汽通量散度剖面图的核心代码：

import xarray as xr
import numpy as np
from metpy.units import units
import metpy.calc as mpcalc
import cartopy.crs as ccrs

def calculate_and_save_moisture_flux_divergence(ds, output_path, time_idx=0):
    """
    计算各层水汽通量散度并保存为 NetCDF 文件
    
    Parameters:
    -----------
    ds : xarray.Dataset
        包含气象数据的数据集
    output_path : str
        保存 NetCDF 文件的路径
    time_idx : int
        时间索引
    """
    # 获取需要的气压层
    levels = [1000, 925, 850, 700, 600, 500, 400, 300]
    
    # 提取所需变量并添加坐标信息
    q = ds['q'].sel(level=levels, time=ds.time[time_idx]) * units('kg/kg')
    u = ds['u'].sel(level=levels, time=ds.time[time_idx]) * units('m/s')
    v = ds['v'].sel(level=levels, time=ds.time[time_idx]) * units('m/s')
    
    # 确保数据具有必要的坐标属性
    q = q.metpy.assign_crs(ccrs.PlateCarree())
    u = u.metpy.assign_crs(ccrs.PlateCarree())
    v = v.metpy.assign_crs(ccrs.PlateCarree())
    
    # 计算水汽通量
    qx = q * u / 9.81
    qy = q * v / 9.81
    
    # 计算网格间距
    dx, dy = mpcalc.lat_lon_grid_deltas(ds['longitude'], ds['latitude'])
    
    # 逐层计算散度
    div_q_list = []
    for level in levels:
        qx_level = qx.sel(level=level)
        qy_level = qy.sel(level=level)
        div_q_level = mpcalc.divergence(qx_level, qy_level, dx=dx, dy=dy)
        div_q_list.append(div_q_level)
    
    # 将各层散度拼接在一起
    div_q = xr.concat(div_q_list, dim='level')
    div_q['level'] = levels
    
    # 将 DataArray 转换为 Dataset 并指定变量名
    div_q_ds = div_q.to_dataset(name='moisture_flux_divergence')
    
    # 保存为 NetCDF 文件
    div_q_ds.to_netcdf(output_path)
    print(f"Moisture flux divergence saved to {output_path}")

# 使用示例：
ds = xr.open_dataset('/home/mw/input/era58091/ERA5-2023-08_pl.nc')
output_path = 'moisture_flux_divergence.nc'
calculate_and_save_moisture_flux_divergence(ds, output_path)

Moisture flux divergence saved to moisture_flux_divergence.nc

使用metpy函数获取剖面

import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
from metpy.units import units
from metpy.interpolate import cross_section
import cartopy.crs as ccrs


nc_path = 'moisture_flux_divergence.nc'
start_point = (45, 100)  # (纬度, 经度)
end_point = (20, 130)    # (纬度, 经度)
div_q = xr.open_dataset(nc_path)
div_q = div_q.metpy.parse_cf().squeeze()
# 获取剖面
cross_section_data = cross_section(div_q, start_point, end_point, steps=100)
print(cross_section_data)

<xarray.Dataset> Size: 9kB
Dimensions:                   (level: 8, index: 100)
Coordinates:
    time                      datetime64[ns] 8B 2023-08-02
  * level                     (level) int64 64B 1000 925 850 700 600 500 400 300
    metpy_crs                 object 8B Projection: latitude_longitude
    longitude                 (index) float64 800B 100.0 100.4 ... 129.8 130.0
    latitude                  (index) float64 800B 45.0 44.79 ... 20.28 20.0
  * index                     (index) int64 800B 0 1 2 3 4 5 ... 95 96 97 98 99
Data variables:
    moisture_flux_divergence  (level, index) float64 6kB 2.086e-08 ... -3.127...

简单绘图

import matplotlib.pyplot as plt
# 创建图形
plt.figure(figsize=(10, 6))
# 绘制水汽通量散度
plot = cross_section_data['moisture_flux_divergence'].plot(
    x='longitude',  # x 轴为经度
    y='level',      # y 轴为气压层
    cmap='RdBu_r',  # 使用红蓝渐变色
    levels=21,      # 等值线数量
    extend='both',  # 扩展颜色条范围
    add_colorbar=True,  # 添加颜色条
    cbar_kwargs={'label': 'Moisture Flux Divergence (10⁻⁵ g kg⁻¹ s⁻¹)'}  # 颜色条标签
)

# 设置 y 轴反向（气压递减）
plt.gca().invert_yaxis()

# 添加标题和标签
plt.title('Moisture Flux Divergence Cross Section', fontsize=14)
plt.xlabel('Longitude (degrees east)', fontsize=12)
plt.ylabel('Pressure (hPa)', fontsize=12)

# 添加网格
plt.grid(True, linestyle='--', alpha=0.5)

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