一瞬又一瞬，累积起来便是一生 | ERA5数据计算垂直积分整层水汽通量散度

用户11172986

发布于 2024-12-30 15:14:22

72700

代码可运行

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

运行总次数：0

代码可运行

ERA5数据计算垂直积分整层水汽通量散度

温馨提示

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

前言

嗨大家好，这次带来的是上期的续集：关于整层的水汽通量散度的计算

在气象学中，整层水汽通量散度（Vertically Integrated Moisture Flux Convergence）是评估大气柱内水分输送及变化的重要指标。它综合考虑了从地表到特定高度范围内的水汽输送情况，对研究降水形成机制和气候变化有着重要意义。本文将基于ERA5再分析数据，利用Python编程语言计算整层水汽通量散度，并进行相关图形绘制。

通过这篇文章的学习，您不仅能理解整层水汽通量散度的物理意义，还将掌握如何使用Python及其生态中的工具来处理和分析气象数据。

项目目标

本项目的主要目标包括：

计算整层水汽通量散度：使用ERA5再分析数据，结合metpy、xarray等库，计算从地面至一定高度（例如1000 hPa至300 hPa）之间的整层水汽通量散度。
绘制整层水汽通量散度分布图：将计算结果可视化，以图形方式展示整层水汽通量散度的空间分布特征。
详细解析单位转换：确保每个步骤中的物理单位清晰明了，保证计算结果的准确性。

项目功能

1. 数据读取与处理

首先，我们需要读取ERA5数据集，提取比湿（q）、经向风速（u）、纬向风速（v）以及气压水平（level）。以下为代码示例：

import xarray as xr
from metpy.units import units

# 读取 ERA5 数据
ds = xr.open_dataset('/home/mw/input/era58091/ERA5-2023-08_pl.nc')
levels = [1000, 925, 850, 700, 600, 500, 400, 300]  # 定义需要整合的高度层次
times = ds['time']  # 获取时间维度

# 提取所需变量
q = ds['q'].sel(level=levels,time=times[0]) * units('kg/kg')  # 比湿
u = ds['u'].sel(level=levels,time=times[0]) * units('m/s')    # 经向风速
v = ds['v'].sel(level=levels,time=times[0]) * units('m/s')    # 纬向风速

2. 整层水汽通量计算

根据所选的高度层次，我们接下来要计算各层的水汽通量，并对其进行垂直积分得到整层水汽通量。

import numpy as np

# 计算各层的水汽通量
qx = q * u / 9.81  # 经向水汽通量
qy = q * v / 9.81  # 纬向水汽通量

# 进行垂直积分获得整层水汽通量
integrated_qx = qx.integrate(coord='level')
integrated_qy = qy.integrate(coord='level')

3. 计算整层水汽通量散度

使用metpy.calc.divergence函数计算整层水汽通量散度。

import metpy.calc as mpcalc

# 计算网格间距
dx, dy = mpcalc.lat_lon_grid_deltas(ds['longitude'], ds['latitude'])

# 计算散度
div_integrated_q = mpcalc.divergence(integrated_qx, integrated_qy, dx=dx, dy=dy)
div_integrated_q = div_integrated_q.rename('integrated_moisture_flux_divergence')*10

4. 绘制整层水汽通量散度图

最后，我们将计算结果绘制成地图形式，以便直观地观察整层水汽通量散度的空间分布。

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

def plot_integrated_moisture_flux_divergence(div_q, extent=None):
    fig = plt.figure(figsize=(10, 8), dpi=300)
    ax = plt.axes(projection=ccrs.PlateCarree())
    
    # 如果提供了extent参数，则设置地图范围
    if extent isnotNone:
        ax.set_extent(extent, crs=ccrs.PlateCarree())

    im = ax.pcolormesh(div_q.longitude, div_q.latitude, div_q * 1e5, 
                       cmap='RdBu', vmin=-30, vmax=30, transform=ccrs.PlateCarree())
    
    # 添加色标
    cbar = fig.colorbar(im, ax=ax, orientation='horizontal', pad=0.05)
    cbar.set_label('Divergence 10-5 *(g•(cm²•s)-1)', fontsize=12)
    
    # 添加地图特征
    ax.add_feature(cfeature.COASTLINE)
    ax.set_xlabel('Longitude', fontsize=12)
    ax.set_ylabel('Latitude', fontsize=12)
    ax.set_title('Integrated Moisture Flux Divergence', fontsize=16)

    plt.show()

# 定义需要显示的地图范围 [最小经度, 最大经度, 最小纬度, 最大纬度]
extent = [90, 140, 15, 55]  # 例如：东亚地区

# 调用绘图函数，并传入extent参数
plot_integrated_moisture_flux_divergence(div_integrated_q, extent=extent)

单位什么就不多说了，就是上期的单位去掉高度层的单位而已

5. 完整计算绘图代码

import xarray as xr
from metpy.units import units
import numpy as np
import metpy.calc as mpcalc
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

# 读取 ERA5 数据
ds = xr.open_dataset('/home/mw/input/era58091/ERA5-2023-08_pl.nc')
levels = [1000, 925, 850, 700, 600, 500, 400, 300]  # 定义需要整合的高度层次
times = ds['time']  # 获取时间维度

# 提取所需变量
q = ds['q'].sel(level=levels, time=times[0]) * units('kg/kg')  # 比湿
u = ds['u'].sel(level=levels, time=times[0]) * units('m/s')    # 经向风速
v = ds['v'].sel(level=levels, time=times[0]) * units('m/s')    # 纬向风速

# 计算各层的水汽通量
qx = q * u / 9.81# 经向水汽通量
qy = q * v / 9.81# 纬向水汽通量

# 进行垂直积分获得整层水汽通量
integrated_qx = qx.integrate(dim='level')
integrated_qy = qy.integrate(dim='level')

# 计算网格间距
dx, dy = mpcalc.lat_lon_grid_deltas(ds['longitude'], ds['latitude'])

# 计算散度
div_integrated_q = mpcalc.divergence(integrated_qx, integrated_qy, dx=dx, dy=dy)
div_integrated_q = div_integrated_q.rename('integrated_moisture_flux_divergence')

def plot_integrated_moisture_flux_divergence(div_q, extent=None):
    fig = plt.figure(figsize=(10, 8), dpi=300)
    ax = plt.axes(projection=ccrs.PlateCarree())
    
    # 如果提供了extent参数，则设置地图范围
    if extent isnotNone:
        ax.set_extent(extent, crs=ccrs.PlateCarree())

    im = ax.pcolormesh(div_q.longitude, div_q.latitude, div_q * 1e6, 
                       cmap='RdBu', vmin=-30, vmax=30, transform=ccrs.PlateCarree())
    
    # 添加色标
    cbar = fig.colorbar(im, ax=ax, orientation='horizontal', pad=0.05)
    cbar.set_label('Divergence 10⁻⁶ *(g•(cm²•s)⁻¹)', fontsize=12)
    
    # 添加地图特征
    ax.add_feature(cfeature.COASTLINE)
    ax.set_xlabel('Longitude', fontsize=12)
    ax.set_ylabel('Latitude', fontsize=12)
    ax.set_title('Integrated Moisture Flux Divergence', fontsize=16)

    plt.show()

# 定义需要显示的地图范围 [最小经度, 最大经度, 最小纬度, 最大纬度]
extent = [90, 140, 15, 55]  

# 调用绘图函数，并传入extent参数
plot_integrated_moisture_flux_divergence(div_integrated_q, extent=extent)