一刻也没有为圣诞的结束而悲伤，下一刻赶来的是ERA5数据计算700hPa水汽通量散度

一刻也没有为圣诞的结束而悲伤，下一刻赶来的是ERA5数据计算700hPa水汽通量散度

ERA5数据计算700hPa水汽通量散度

温馨提示

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

前言

在小编为论文伏案工作的时候，想必读者们在欢乐的海洋中畅游吧。

言归正传，一个读者朋友说不知ERA5怎么计算水汽通量散度，而且单位老是弄不明白，那咱们好好说道说道。

在气象领域，水汽通量散度（Moisture Flux Divergence，简称MFD）是描述大气中水汽输送特性的重要物理量之一。它反映了某一地区水汽的净流入或流出情况，与降水、蒸发等水文循环过程密切相关。在本篇文章中，我们将基于 ERA5 再分析数据，利用 Python 代码计算 700 hPa 水汽通量散度，并绘制相关图形。

通过这篇文章，你不仅能够理解水汽通量散度的物理意义，还能够掌握如何使用 Python 和相关气象库计算和可视化这一关键变量。

项目目标

本项目的主要目标是：

1. 计算水汽通量散度 使用 ERA5 再分析数据，结合 metpy 和 xarray 等 Python 库，计算大气 500 hPa 层的水汽通量散度。
2. 绘制水汽通量散度分布图 将计算结果可视化，展示水汽通量散度的空间分布特征。
3. 逐步分析水汽通量散度的单位 对每一个计算步骤的物理单位进行剖析，确保物理意义清晰且计算结果准确。

项目功能

1. 数据读取与处理

首先我们读取 ERA5 气象数据，包括比湿（q）、经向风速（u）、纬向风速（v）等变量。以下是代码示例：

代码语言：javascript代码运行次数：0运行复制

import xarray as xr
import numpy as np
from metpy.units import units

# 读取 ERA5 数据
ds = _dataset('/home/mw/input/era58091/ERA5-202-08_')
q = ds['q'].sel(level=700, time='202-08-02T00:00:00.000000000') * units('kg/kg')
u = ds['u'].sel(level=700, time='202-08-02T00:00:00.000000000') * units('m/s')
v = ds['v'].sel(level=700, time='202-08-02T00:00:00.000000000') * units('m/s')
lat = ds['latitude']
lon = ds['longitude']

单位分析：

• 比湿（q）：单位为 kg/kg，表示单位空气质量中的水汽质量。
• 风速（u 和 v）：单位为 m/s，分别为经向和纬向风速。

水汽通量（qx, qy）的单位由 q 和风速的乘积决定：

• 水汽通量单位： qx 或 qy = q × 风速 /g 单位为 g•(cm•hPa•s) -1 ，表示单位时间内通过单位面积的水汽质量。

具体推算请看

此处小编采用了国家气象局的单位，即采用了g和cm作为单位，因此仅需乘10

国家气象局的定义：.html

2. 计算水汽通量散度

接下来，我们利用 提供的工具计算水汽通量散度。水汽通量散度公式为：

国家气象局给出的单位是g•(cm2•hPa•s) -1

代码如下：

代码语言：javascript代码运行次数：0运行复制

import  as mpcalc
g= 9.81
# 计算水汽通量 ,单位转为
qu = q * u /g 
qv = q * v /g 

# 计算网格间距
dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat)

# 计算散度
div_q = mpcalc.divergence(qu, qv, dx=dx, dy=dy)
div_q = div_q.rename('moisture_flux_divergence') 
print(div_q )

代码语言：javascript代码运行次数：0运行复制

<xarray.DataArray 'moisture_flux_divergence' (latitude: 241, longitude: 61)>
<Quantity([[ .44509074e-09 -1.61719059e-09 -.68101855e-09 ... -1.72974808e-09
  -1.7512951e-09 -1.8992010e-09]
 [-1.02897068e-08 -1.12989667e-08 -1.21819860e-08 ... -.286579e-10
   1.0769980e-09  .81845815e-09]
 [-1.81800527e-08 -1.77221528e-08 -1.817400e-08 ...  .8684615e-09
   .6089579e-09  .548750e-09]
 ...
 [-4.78795414e-09 -1.2207051e-08 -1.666407e-08 ...  .1557467e-09
   4.996780e-09  1.7627656e-09]
 [-1.7789446e-08 -1.59844e-08 -1.458801e-08 ...  .22479922e-09
  -2.46608056e-09 -1.81507760e-08]
 [-1.76600870e-08 -1.17116196e-08 -6.222871e-09 ...  .94095656e-09
  -8.81956850e-09 -2.09485875e-08]], '1 / second')>
Coordinates:
  * longitude  (longitude) float2 90.0 90.25 90.5 90.75 ... 179.5 179.8 180.0
  * latitude   (latitude) float2 70.0 69.75 69.5 69.25 ... 10.5 10.25 10.0
    level      int2 700
    time       datetime64[ns] 202-08-02

单位分析：

水汽通量散度的单位由水汽通量（kg/m²/s）和空间导数（1/m）共同决定：

• 空间导数单位：1/m
• 水汽通量散度单位： ∇·Q = 水汽通量 × 空间导数 此时单位为 **kg/(m²·s·pa)**，表示单位时间内单位面积的水汽质量净变化率。 需要将其转为g•(cm2•hPa•s) -1，那么需要div_q需要乘10

代码语言：javascript代码运行次数：0运行复制

div_q = div_q*10

. 绘制水汽通量散度图

将计算结果可视化，使用 cartopy 库绘制水汽通量散度的空间分布图。代码示例如下：

代码语言：javascript代码运行次数：0运行复制

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

# 绘制散度图
def plot_moisture_flux_divergence(div_q):
    fig = plt.figure(figsize=(8, 6), dpi=00)
    ax = plt.axes(projection=ccrs.PlateCarree())

    # 绘制散度场
    im = ax.pcolormesh(div_q.longitude, div_q.latitude, div_q * 1e6, cmap='RdBu', vmin=-1, vmax=1)
    
    # 添加标
    cbar = (im, ax=ax, orientation='horizontal', pad=0.05)
    cbar.set_label('Divergence 10 -6 *(g•(cm2•hPa•s) -1)', fontsize=12)
    
    # 添加地图特征
    ax.add_feature(cfeature.COASTLIE)
    ax.set_title('Moisture Flux Divergence at 700 hPa', fontsize=14)
    ax.set_xlabel('Longitude', fontsize=12)
    ax.set_ylabel('Latitude', fontsize=12)
    
    # 保存图片
    plt.savefig('moisture_flux_divergence.png', bbox_inches='tight')

# 调用绘图函数
plot_moisture_flux_divergence(div_q)

代码语言：javascript代码运行次数：0运行复制

/opt/conda/lib/python.9/site-packages/cartopy/mpl/geoaxes.py:1797: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensi as C is deprecated since ..  Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading'].  This will become an error two minor releases later.
  result = matplotlib.axes.Axes.pcolormesh(self, *args, **kwargs)

.5 绘图升级版

如需查看代码请按照温馨提示操作

4. 完整计算绘图代码（不含升级版）代码语言：javascript代码运行次数：0运行复制

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

# 计算水汽通量散度
def calculate_moisture_flux_divergence(filepath, level, time):
    ds = _dataset(filepath)
    q = ds['q'].sel(level=level, time=time) * units('kg/kg')
    u = ds['u'].sel(level=level, time=time) * units('m/s')
    v = ds['v'].sel(level=level, time=time) * units('m/s')
    lat = ds['latitude']
    lon = ds['longitude']
    
    # 计算水汽通量
    qx = q * u
    qy = q * v
    
    # 计算网格间距
    dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat)
    
    # 计算散度
    div_q = mpcalc.divergence(qx, qy, dx=dx, dy=dy)
    div_q = div_q.rename('moisture_flux_divergence')
    return div_q

# 绘制水汽通量散度图
def plot_moisture_flux_divergence(div_q):
    fig = plt.figure(figsize=(8, 6), dpi=00)
    ax = plt.axes(projection=ccrs.PlateCarree())
    im = ax.pcolormesh(div_q.longitude, div_q.latitude, div_q * 1e6, cmap='RdBu', vmin=-5, vmax=5)
    cbar = (im, ax=ax, orientation='horizontal', pad=0.05)
    cbar.set_label('Divergence (10⁻⁶ kg/m²/s)', fontsize=12)
    ax.add_feature(cfeature.COASTLIE)
    ax.set_title('Moisture Flux Divergence at 700 hPa', fontsize=14)
    plt.savefig('moisture_flux_divergence.png', bbox_inches='tight')

# 主程序
filepath = '/home/mw/input/era58091/ERA5-202-08_'
level = 500# hPa
time = '202-08-02T00:00:00.000000000'
div_q = calculate_moisture_flux_divergence(filepath, level, time)
plot_moisture_flux_divergence(div_q)

小结

通过本项目，我们达成了如下目标：

1. 计算了 ERA5 数据中 500 hPa 大气层的水汽通量散度，并对物理单位进行了详细分析。
2. 绘制了水汽通量散度的空间分布图，直观展示了水汽的辐合与散发区域。
3. 优化代码结构，提升了代码的通用性与可读性。

本文参与 腾讯云自媒体同步曝光计划，分享自。原始发表：2024-12-26，如有侵权请联系 cloudcommunity@tencent 删除importtime变量可视化数据