1. 项目概述：从气象数据到直观洞察

在气象分析和气候研究中，我们常常需要理解大气中水分的“来龙去脉”。水分是能量传输和天气过程（如暴雨、台风）的核心载体，但仅知道某个高度上的湿度或风速是远远不够的。我们需要一个能综合反映整层大气水汽输送强度和汇聚、辐散情况的物理量，这就是整层水汽通量和整层水汽通量散度。前者告诉你“有多少水汽在流动”，后者则揭示“水汽在哪里堆积或流失”，这对于暴雨落区预报、干旱监测乃至气候模式评估都至关重要。

然而，从原始的格点数据（比如再分析资料ERA5、NCEP/NCAR）到一张清晰明了的分析图，中间隔着一条由公式、编程和可视化构成的“鸿沟”。很多刚入门的研究生或业务人员会被困在数据处理和编程调试的细节里。这个项目的目的，就是手把手带你用Python这条“捷径”，打通从理论公式到科研绘图的全流程。我们将基于最常见、最权威的ERA5再分析数据，详细拆解每一个计算步骤，并最终用Cartopy和Matplotlib绘制出专业级的空间分布图。无论你是大气科学、水文气象专业的学生，还是从事相关领域的工程师，这篇内容都将为你提供一个可直接复现的“工具箱”，让你能独立完成从数据下载到成果展示的完整分析。

2. 核心概念与数据准备

在动手写代码之前，我们必须把物理概念和数据处理流程理清楚。这一步基础打得牢，后面的计算和绘图才能顺畅无误。

2.1 物理概念深度解析

整层水汽通量，其物理意义是单位时间内、通过单位宽度垂直大气柱所输送的水汽质量。它不是一个单层的数据，而是从地面到大气顶的垂直积分。公式表示为：

[ \vec{Q} = \frac{1}{g} \int_{p_{top}}^{p_{sfc}} q \vec{V} , dp ]

这里，(\vec{Q}) 是整层水汽通量矢量（通常包含Q_u和Q_v两个分量），单位是 (kg \cdot m^{-1} \cdot s^{-1})。(g) 是重力加速度，(q) 是比湿（特定湿度），(\vec{V}) 是水平风矢量（包含u和v分量），(p) 是气压，积分从大气顶气压(p_{top})到地表气压(p_{sfc})。简单理解，它就是把每一层“随风而动”的水汽累加起来，得到一个能代表整层大气水汽输送总效果的矢量。箭头方向代表输送方向，箭头长度代表输送强度。

整层水汽通量散度，则是这个通量矢量的散度，公式为：

[ \nabla \cdot \vec{Q} = \frac{\partial Q_u}{\partial x} + \frac{\partial Q_v}{\partial y} ]

其单位是 (kg \cdot m^{-2} \cdot s^{-1})。它的物理意义更为关键：散度为负（辐合），表示该区域水汽收入大于支出，水汽在此汇聚，是降水发生的有利条件；散度为正（辐散），表示水汽净流出，不利于降水。因此，水汽通量散度场是预报员寻找暴雨潜在落区的关键参考图之一。

注意：在实际计算中，我们处理的是离散的格点数据。因此，“垂直积分”转化为对各个气压层的求和，“水平散度”则转化为对相邻格点值的差分计算。理解这个离散化的过程，是避免计算结果出现物理上不合理现象（如虚假的强辐散辐合中心）的基础。

2.2 数据源选择与下载

工欲善其事，必先利其器。数据质量直接决定分析结果的可靠性。对于大尺度气象分析，欧洲中期天气预报中心（ECMWF）的ERA5再分析资料是目前综合质量最高、应用最广的数据集之一。它提供了高时空分辨率、多要素且物理一致的数据。

数据下载策略（以ERA5为例）：

1. 访问平台：通过ECMWF的CDS（Climate Data Store）API进行下载。你需要先在官网注册账号并获取API密钥。
2. 确定要素：计算水汽通量，我们需要以下变量：
   • specific_humidity（比湿 q）
   • u_component_of_wind（纬向风 u）
   • v_component_of_wind（经向风 v）
   • surface_pressure（地表气压，用于确定积分下界）
3. 时空范围与层次：
   • 时间：选择你需要分析的日期和时间（例如，一次暴雨过程期间）。ERA5提供逐小时数据。
   • 空间：选定经纬度范围。为了计算水平散度，区域范围应比最终成图范围稍大一圈，避免边界缺值。
   • 层次：选择气压层数据。通常需要从高层（如100 hPa）到近地面（如1000 hPa）的多层数据。ERA5提供了从1000 hPa到1 hPa的多个标准层。对于水汽通量积分，积分到100 hPa通常已足够，因为高层水汽含量极低。
4. 使用CDS API工具：推荐使用cdsapi这个Python库。你需要编写一个请求脚本，指定上述参数。一个典型的请求字典如下所示：

   import cdsapi c = cdsapi.Client() c.retrieve('reanalysis-era5-pressure-levels', { 'product_type': 'reanalysis', 'variable': ['specific_humidity', 'u_component_of_wind', 'v_component_of_wind'], 'pressure_level': ['1000', '925', '850', '700', '600', '500', '400', '300', '250', '200', '150', '100'], 'year': '2023', 'month': '07', 'day': ['20', '21', '22'], 'time': ['00:00', '06:00', '12:00', '18:00'], 'area': [50, 70, 15, 140], # 北纬，西经，南纬，东经 'format': 'netcdf', }, 'era5_data_pl.nc') c.retrieve('reanalysis-era5-single-levels', { 'product_type': 'reanalysis', 'variable': 'surface_pressure', 'year': '2023', 'month': '07', 'day': ['20', '21', '22'], 'time': ['00:00', '06:00', '12:00', '18:00'], 'area': [50, 70, 15, 140], 'format': 'netcdf', }, 'era5_data_sl.nc')

   这里我们将气压层数据和地表气压数据分开下载，因为它们在CDS中属于不同的数据集。

实操心得：下载数据可能是最耗时的一步，尤其是高时空分辨率、长时序的数据。建议首次调试代码时，先下载一个小范围、短时段的数据进行测试，待整个流程跑通后，再扩展下载完整数据。另外，注意CDS有排队系统，请求可能不会立即完成，需要等待。

3. 计算流程详解与Python实现

有了数据，我们就可以进入核心的计算环节。这一部分将把理论公式一步步转化为可执行的Python代码，并解释每一个操作背后的气象学和数学考量。

3.1 环境配置与库导入

首先，确保你的Python环境安装了必要的科学计算和地理绘图库。推荐使用Anaconda管理环境。

# 创建一个新的conda环境（可选） conda create -n meteorology python=3.9 conda activate meteorology # 安装核心库 conda install -c conda-forge xarray dask netCDF4 cartopy matplotlib numpy scipy

xarray是处理NetCDF格式气象数据的利器，它能够优雅地处理带标签的多维数组。cartopy则是专业的地图投影和地理绘图库。

计算开始前，先导入所有需要的库：

import xarray as xr import numpy as np import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature from scipy import constants import warnings warnings.filterwarnings('ignore') # 忽略一些不影响计算的警告

3.2 数据读取与预处理

使用xarray打开下载好的NetCDF文件。

# 读取数据 ds_pl = xr.open_dataset('era5_data_pl.nc') # 气压层数据 ds_sl = xr.open_dataset('era5_data_sl.nc') # 单层（地表气压）数据 # 查看数据结构 print(ds_pl) print(ds_sl)

你会看到数据集的维度（如time,level,latitude,longitude）和变量。我们需要从中提取出计算所需的变量：比湿q，纬向风u，经向风v，以及地表气压sp。

关键预处理步骤：

1. 单位统一：ERA5的比湿单位是kg/kg，风分量单位是m/s，地表气压单位是Pa。这些单位与公式要求一致，通常无需转换。但务必确认。
2. 坐标对齐：确保两个数据集（气压层和单层）的time，latitude，longitude坐标完全一致。xarray的运算会自动对齐坐标，但为保险起见，可以用.sel或.interp进行精确匹配。
3. 处理缺失值：再分析数据通常质量很高，但沿海或复杂地形处可能有缺测值。可以用.fillna(0)或插值方法处理，但需谨慎，因为这会改变物理场。对于水汽计算，海洋上的缺测可以填充为0（无陆地），但陆地上的处理需要根据实际情况判断。

3.3 整层水汽通量的垂直积分计算

这是最核心的一步。我们采用梯形法对气压进行垂直积分，这是气象上处理层结数据积分的常用且精度较高的方法。

# 定义常数 g = constants.g # 从scipy.constants获取标准重力加速度，约9.80665 m/s^2 # 从数据集中提取变量，假设变量名就是标准的‘q’， ‘u’， ‘v’， ‘sp’ q = ds_pl['q'] # 比湿 [kg/kg] u = ds_pl['u'] # 纬向风 [m/s] v = ds_pl['v'] # 经向风 [m/s] sp = ds_sl['sp'] # 地表气压 [Pa] # 获取气压层坐标（单位：Pa）。ERA5气压层单位是hPa，需要转换为Pa。 pressure_levels = ds_pl['level'] * 100 # 转换为Pa # 注意：ERA5数据中，level坐标是从高到低（1000, 925,...），这符合气压从地面向高空递减的物理事实，也方便积分。 # 计算各层的水汽通量分量 (q*u/g) 和 (q*v/g) # 这里先不积分，只计算被积函数 integrand_u = q * u / g integrand_v = q * v / g # 关键：垂直积分（梯形法） # 思路：对每一层，其上下气压差（dp）作为权重。对于最顶层和最底层需要特殊处理。 # 我们假设积分从最顶层气压（p_top）到地表气压（sp），但sp是变化的（随格点和时间）。 # 因此，我们需要构建一个与q形状相同的压力矩阵，其中每一层的“层中气压”用于差分计算。 # 扩展pressure_levels的维度，使其与q的维度对齐（除了level维） # pressure_levels是一个一维数组，我们需要将其扩展为多维 pressures = xr.DataArray(pressure_levels.values, dims=['level'], coords={'level': pressure_levels}) # 现在pressures是一个带‘level’坐标的一维DataArray # 计算各层之间的气压差 dp # 使用diff函数，但注意diff会减少一个维度。我们需要在顶层和底层进行外推。 dp = pressures.diff('level') # 这是层与层之间的气压差 # 为了进行梯形积分，我们需要每个层“节点”对应的“层厚”。 # 一种常见方法是：对于内部层，dp取前后两半的平均；对于顶层和底层，dp取相邻层差的一半。 # 更精确的做法是构建一个与q同维度的dp_full数组。 # 构建全层的dp权重数组 dp_values = np.zeros_like(pressure_levels.values, dtype=float) levs = pressure_levels.values # 内部层：使用梯形法的权重 (p_{i+1} - p_{i-1}) / 2 for i in range(1, len(levs)-1): dp_values[i] = (levs[i-1] - levs[i+1]) / 2.0 # 顶层（i=0）：假设其上界气压为0？不对。更合理的是用 (levs[0] - levs[1])，但这是从levs[0]到levs[1]的整层厚度。 # 在梯形法中，对于最顶层，我们通常用 (levs[0] - levs[1]) / 2 作为该层“节点”的代表厚度？这比较复杂。 # 实际操作中，一个更稳健且物理意义明确的方法是： # 1. 将地表气压sp作为积分下界。 # 2. 对于每个格点，找出sp所在的气压层区间。 # 3. 只对从最顶层到sp所在层之间的层次进行积分。 # 鉴于上述复杂性，一个在科研中广泛采用的简化且有效的方法是： # 假设我们选取的标准气压层（如1000, 925,...100 hPa）足够密集，能够代表大气垂直结构。 # 我们直接对所有层次进行积分，积分下界统一取为最底层的气压（如1000 hPa）。 # 但这样忽略了地表气压变化的影响。对于高原地区（地表气压可能只有700 hPa），用1000 hPa作为下界会错误地积分不存在的空气柱。 # 因此，必须考虑真实的地表气压。 # 下面展示一个考虑真实地表气压的积分方法： # 步骤：对于每一个格点、每一个时次，动态确定需要积分的层次。 # 由于涉及复杂循环，效率较低。我们可以利用xarray的向量化运算进行优化。 def integrate_vertical(integrand, pressure, surface_pressure): """ 使用梯形法则进行垂直积分，积分上限为各层气压，下限为地表气压。 参数： integrand: DataArray，被积函数，维度需包含‘level’ pressure: DataArray，气压层坐标（单位：Pa），一维 surface_pressure: DataArray，地表气压（单位：Pa），维度需与integrand除level外的维度一致 返回： integral: DataArray，垂直积分结果，维度不含‘level’ """ # 确保pressure是递减的（从地面向高空） if pressure[0] < pressure[-1]: pressure = pressure[::-1] integrand = integrand.sel(level=pressure) # 重新索引，确保对应关系正确 # 将地表气压扩展到与integrand相同的维度结构（除了level） # 这里surface_pressure需要与integrand的‘time’， ‘latitude’， ‘longitude’对齐 sp_expanded = surface_pressure.broadcast_like(integrand.isel(level=0)) # 创建一个掩膜，标识出气压值大于地表气压的层次（即在地表以下的层次，不参与积分） # 因为气压向下增加，所以“大于”地表气压的层次是无效的。 mask = pressure > sp_expanded # 将无效层次的被积函数设为0 integrand_masked = integrand.where(~mask, 0) # 对level维进行梯形积分 # xarray的integrate方法可以方便地进行梯形积分 integral = integrand_masked.integrate(coord='level') return integral # 应用函数计算整层水汽通量 Q_u = integrate_vertical(integrand_u, pressures, sp) Q_v = integrate_vertical(integrand_v, pressures, sp) # 给计算结果添加单位属性 Q_u.attrs['units'] = 'kg m-1 s-1' Q_v.attrs['units'] = 'kg m-1 s-1'

这段代码实现了一个考虑真实地表气压的、健壮的垂直积分函数。它避免了在高原地区积分“虚假”气柱的问题，计算结果更可靠。

3.4 整层水汽通量散度的计算

计算散度需要用到水平差分。在球坐标系（经纬网格）上，计算散度的公式需要考虑地球曲率和纬圈变化。

对于经纬网格上的矢量场(A_x, A_y)（这里A_x沿经线方向，A_y沿纬线方向），其散度公式为：

[ \nabla \cdot \vec{A} = \frac{1}{R \cos \phi} \left[ \frac{\partial A_x}{\partial \lambda} + \frac{\partial (A_y \cos \phi)}{\partial \phi} \right] ]

其中，(R)是地球半径（约6371 km），(\phi)是纬度（弧度），(\lambda)是经度（弧度）。在我们的计算中，(A_x = Q_u)（纬向通量，沿纬圈方向，对应经度变化），(A_y = Q_v)（经向通量，沿经线方向，对应纬度变化）。注意：气象上通常定义u为东西风（向东为正），对应经度方向；v为南北风（向北为正），对应纬度方向。因此，公式中的(A_x)对应(Q_u)，(A_y)对应(Q_v)。但有些文献定义可能相反，务必与风场分量的定义保持一致。

因此，水汽通量散度(D)的计算公式为：

[ D = \nabla \cdot \vec{Q} = \frac{1}{R \cos \phi} \left[ \frac{\partial Q_u}{\partial \lambda} + \frac{\partial (Q_v \cos \phi)}{\partial \phi} \right] ]

在离散的格点上，我们使用中心差分来计算偏导数。

def calculate_divergence(Q_u, Q_v, lat, lon, radius=6371000.): """ 计算球面上矢量场（Q_u, Q_v）的散度。 参数： Q_u, Q_v: DataArray，矢量场的两个分量。假设Q_u对应东西方向（经度λ），Q_v对应南北方向（纬度φ）。 lat: DataArray，纬度坐标（单位：度） lon: DataArray，经度坐标（单位：度） radius: 地球半径，单位米 返回： div: DataArray，散度场 """ # 将经纬度转换为弧度 lat_rad = np.deg2rad(lat) lon_rad = np.deg2rad(lon) # 计算纬度的余弦（每个格点一个值） cos_lat = np.cos(lat_rad) # 计算经度方向上的差分 d(Q_u)/dλ # 使用中心差分，注意处理周期边界（经度360度循环） dQ_u_dlambda = Q_u.differentiate('longitude') # xarray的differentiate方法计算中心差分 # differentiate 计算的是 ΔQ_u / Δlon (单位：度)，我们需要 ΔQ_u / Δλ (单位：弧度) # Δlon (度) 转换为 Δλ (弧度): Δλ = Δlon * (π/180) # 所以 d(Q_u)/dλ = (ΔQ_u/Δlon) * (180/π) dQ_u_dlambda = dQ_u_dlambda * (180.0 / np.pi) # 因为 differentiate('longitude') 给出的是 per degree， 乘以 (180/π) 得到 per radian？ 等等，这里需要仔细推导。 # 实际上：d(Q_u)/dλ = [Q_u(λ+Δλ) - Q_u(λ-Δλ)] / (2 * Δλ) # 而 xarray 的 differentiate 是：d(Q_u)/d(lon) = [Q_u(lon+Δlon) - Q_u(lon-Δlon)] / (2 * Δlon)，其中lon是度。 # 因为 λ = lon * (π/180)， 所以 Δλ = Δlon * (π/180) # 因此，d(Q_u)/dλ = d(Q_u)/d(lon) * (d(lon)/dλ) = d(Q_u)/d(lon) * (180/π) # 所以上面的转换是正确的。 # 计算纬度方向上的差分 d(Q_v cosφ)/dφ # 先计算 Q_v_cos = Q_v * cosφ Q_v_cos = Q_v * cos_lat # 对纬度求差分 dQ_v_cos_dphi = Q_v_cos.differentiate('latitude') # 单位：每度 # 同理， φ = lat * (π/180), Δφ = Δlat * (π/180) # d(Q_v_cos)/dφ = d(Q_v_cos)/d(lat) * (180/π) dQ_v_cos_dphi = dQ_v_cos_dphi * (180.0 / np.pi) # 计算散度 D = 1/(R cosφ) * [ dQ_u/dλ + d(Q_v_cos)/dφ ] divergence = (dQ_u_dlambda + dQ_v_cos_dphi) / (radius * cos_lat) divergence.attrs['units'] = 'kg m-2 s-1' return divergence # 获取经纬度坐标 lat = Q_u.latitude lon = Q_u.longitude # 计算散度 div_Q = calculate_divergence(Q_u, Q_v, lat, lon)

这个函数严格遵循了球坐标下的散度公式，并使用xarray内置的.differentiate()方法进行中心差分计算，代码简洁且高效。注意，.differentiate()会自动处理网格间距（即Δlon和Δlat），前提是你的经纬度坐标是等间距的。ERA5数据通常是规则网格，所以适用。

注意事项：散度计算对数据噪声非常敏感，尤其是在小尺度上。计算结果可能会出现一些小的、不规则的极值点。在绘图前，通常需要进行适度的平滑处理（如使用高斯滤波或简单滑动平均），以突出大尺度的辐合辐散特征，这更符合天气尺度分析的目的。可以使用scipy.ndimage.gaussian_filter进行平滑，但要注意平滑强度不宜过大，以免失真。

4. 使用Cartopy进行专业气象绘图

计算得到Q_u,Q_v,div_Q这些物理量场后，如何将它们清晰、美观、专业地呈现出来，是科研成果表达的关键一步。我们将使用Cartopy和Matplotlib来创建包含地图背景、风矢量和散度填色的综合图表。

4.1 绘图基础与地图投影

Cartopy的核心是地图投影。气象上最常用的是兰勃特正形圆锥投影（Lambert Conformal Conic）和极射赤面投影（Polar Stereographic），它们在中高纬度地区变形小，适合分析温带天气系统。对于中国区域，兰勃特投影是标准选择。

# 设置图形和地图投影 fig = plt.figure(figsize=(14, 10)) # 创建兰勃特投影，中心点通常设在中国中部（如105°E， 35°N） proj = ccrs.LambertConformal(central_longitude=105, central_latitude=35, standard_parallels=(25, 47)) ax = fig.add_subplot(1, 1, 1, projection=proj) # 设置地图范围（西经， 东经， 南纬， 北纬） map_extent = [70, 140, 15, 55] ax.set_extent(map_extent, crs=ccrs.PlateCarree()) # 注意：set_extent的坐标范围需用PlateCarree经纬度表示 # 添加地理特征 ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.8) ax.add_feature(cfeature.BORDERS.with_scale('50m'), linewidth=0.5, linestyle=':') ax.add_feature(cfeature.OCEAN, color='lightcyan', alpha=0.6) ax.add_feature(cfeature.LAND, color='wheat', alpha=0.3) # 添加省界（需要中国省界shapefile，这里用河流替代示意，实际应用中可加载自定义数据） ax.add_feature(cfeature.RIVERS.with_scale('50m'), color='blue', linewidth=0.5, alpha=0.5) # 添加经纬网格线 gl = ax.gridlines(draw_labels=True, dms=True, x_inline=False, y_inline=False, linewidth=0.5, color='gray', alpha=0.5, linestyle='--') gl.top_labels = False # 不显示顶部标签 gl.right_labels = False # 不显示右侧标签 gl.xlabel_style = {'size': 10} gl.ylabel_style = {'size': 10}

4.2 绘制水汽通量散度填色图

散度场通常用填色图表示，并使用红-蓝发散色系（如RdBu_r），其中蓝色表示辐合（负值，水汽汇聚），红色表示辐散（正值，水汽流失）。

# 选择一个时次进行绘图，例如第一个时次 time_idx = 0 div_plot = div_Q.isel(time=time_idx) # 定义散度填色的等级和色标 # 先计算数据的近似范围，以确定合理的色阶 div_max = np.nanmax(np.abs(div_plot.values)) levels = np.linspace(-div_max, div_max, 21) # 生成从-最大绝对值到+最大绝对值的21个等级 # 或者手动设置一个固定范围，例如对于夏季强降水过程，散度量级可能在 -50 到 50 * 1e-5 kg m-2 s-1 之间 # levels = np.arange(-50, 51, 5) * 1e-5 # 绘制填色图 # 注意：绘图时数据坐标是经纬度（PlateCarree），但地图投影是Lambert。 # 需要使用 transform 参数告诉 cartopy 数据的坐标系。 cf = ax.contourf(lon, lat, div_plot, levels=levels, cmap='RdBu_r', transform=ccrs.PlateCarree(), extend='both') # extend='both' 表示色标向两端延伸 # 添加色标 cbar = plt.colorbar(cf, ax=ax, orientation='horizontal', pad=0.05, aspect=40, shrink=0.8) cbar.set_label('Integrated Water Vapor Flux Divergence (kg m$^{-2}$ s$^{-1}$)', fontsize=12)

4.3 叠加水汽通量矢量箭头（风羽）

水汽通量矢量（Q_u, Q_v）用箭头（风羽）表示，可以直观看到水汽输送的方向和强度。由于通量矢量通常很大，直接画箭头会过于密集，需要降采样。

# 降采样，每隔几个格点画一个箭头 stride = 5 # 根据你的格点密度调整，密度大则stride大一些 Q_u_plot = Q_u.isel(time=time_idx) Q_v_plot = Q_v.isel(time=time_idx) lon_sub = lon[::stride] lat_sub = lat[::stride] Q_u_sub = Q_u_plot[::stride, ::stride] Q_v_sub = Q_v_plot[::stride, ::stride] # 计算箭头的大小（速度），用于归一化箭头长度，避免因量级过大导致箭头过长 speed = np.sqrt(Q_u_sub**2 + Q_v_sub**2) # 对箭头进行缩放，使得图形美观。scale参数需要反复调试。 scale = 3e6 # 这是一个经验值，需要根据你的数据量级调整 # 如果箭头太密或太长，增大scale；如果箭头太稀疏或太短，减小scale。 # 绘制箭头 quiver = ax.quiver(lon_sub.values, lat_sub.values, Q_u_sub.values, Q_v_sub.values, speed.values, # 用速度着色箭头 cmap='plasma', scale=scale, scale_units='inches', transform=ccrs.PlateCarree(), width=0.002, headwidth=3, headlength=4) # 为风羽图添加一个独立的色标（可选，表示通量强度） # cbar2 = plt.colorbar(quiver, ax=ax, orientation='vertical', pad=0.1, shrink=0.8) # cbar2.set_label('IVT Magnitude (kg m$^{-1}$ s$^{-1}$)', fontsize=12)

4.4 添加标题与修饰

# 添加标题 plot_time = pd.to_datetime(str(div_plot.time.values)).strftime('%Y-%m-%d %H:%M UTC') ax.set_title(f'Integrated Vapor Transport and Divergence\n{plot_time}', fontsize=16, fontweight='bold', pad=20) # 添加文本标注，如研究区域或说明 ax.text(0.02, 0.98, 'Blue: Convergence (Moisture Sink)\nRed: Divergence (Moisture Source)', transform=ax.transAxes, verticalalignment='top', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.8), fontsize=10) plt.tight_layout() # 保存图像 plt.savefig('IVT_Divergence.png', dpi=300, bbox_inches='tight') plt.show()

实操心得：绘图的调试往往比计算更耗时。关键点在于：1)投影选择要适合你的分析区域；2)颜色映射要符合气象学惯例（如散度用红蓝，降水用蓝白红）；3)箭头缩放（scale参数）需要多次尝试才能达到疏密适中、长度合适的效果；4) 图形要素（海岸线、网格线、标签）要清晰但不喧宾夺主。建议将绘图代码封装成函数，便于对不同时次或个例进行批量出图。

5. 常见问题、排查技巧与性能优化

在实际操作中，你几乎一定会遇到各种报错和不如预期的结果。这里汇总了一些典型问题及其解决方案。

5.1 计算相关的问题

问题1：垂直积分结果量级异常大或小。

• 可能原因：
  1. 单位错误：检查比湿q的单位是否为kg/kg（无量纲）。ERA5数据通常是正确的。检查重力加速度g的单位是m/s²。
  2. 积分上下限错误：确认积分是从高压（地面）向低压（高空）进行。检查pressure_levels数组是否按降序排列。确保surface_pressure参与积分下界的判断。
  3. 地表气压处理不当：在高原地区，如果直接用1000 hPa作为积分下界，会严重高估水汽通量。务必使用integrate_vertical函数中实现的、基于真实地表气压的掩膜方法。
• 排查方法：打印出单个格点（例如一个海洋点和一个高原点）的pressure_levels、sp值，以及积分过程中被积函数和掩膜的情况，进行人工核对。

问题2：散度场出现棋盘状噪声或极端值。

• 可能原因：
  1. 微分对噪声敏感：原始风场和湿度场本身有小尺度噪声，微分会将其放大。
  2. 网格非均匀：虽然ERA5是规则网格，但如果你使用的数据是变分辨率或跳点采样的，直接使用中心差分公式会不准。
  3. 边界效应：在计算区域边界，中心差分缺少一侧的格点，可能导致异常值。
• 解决方案：
  1. 平滑滤波：在计算散度前，对Q_u和Q_v进行轻微的高斯平滑。

     from scipy.ndimage import gaussian_filter Q_u_smooth = xr.apply_ufunc(gaussian_filter, Q_u, input_core_dims=[['latitude', 'longitude']], output_core_dims=[['latitude', 'longitude']], kwargs={'sigma': 1.0}) # sigma控制平滑强度 Q_v_smooth = ... # 同理

  2. 使用更稳健的差分方法：可以考虑使用二次精度的差分格式，或者直接调用气象专用库（如windspharm）来计算球面散度，后者经过充分测试，更为可靠。
  3. 剔除边界：计算散度后，将边界附近一圈格点的值设为NaN。

问题3：内存不足，处理大数据时卡死。

• 原因：ERA5高时空分辨率数据体积庞大，直接读入内存可能超过限制。
• 解决方案：
  1. 使用Dask进行分块计算：xarray与dask无缝集成。在打开数据集时使用chunks参数。

     ds_pl = xr.open_dataset('era5_data_pl.nc', chunks={'time': 1, 'level': 5, 'latitude': 100, 'longitude': 100})

     这样数据以“惰性”方式加载，计算任务会被图优化并分块执行，极大减少内存峰值。
  2. 分时次处理：如果不需要做时间平均，可以循环处理每个时次，处理完一个就保存或绘图，然后释放内存。
  3. 降低空间分辨率：对于大范围气候分析，可以先用.coarsen()或.interp()方法将数据插值到更粗的网格上。

5.2 绘图相关的问题

问题1：箭头（quiver）太密集或太稀疏，看不清。

• 解决：调整stride参数和scale参数。stride控制跳过的格点数，scale控制箭头的整体长度。两者配合调整。一个经验法则是，让最强的箭头长度大约等于图上5个经度/纬度的距离。

问题2：填色图（contourf）在Cartopy投影上出现奇怪的空洞或扭曲。

• 原因：数据可能包含NaN值，或者地图投影的转换在数据边界处出现问题。
• 解决：
  1. 确保绘图数据在绘图区域内没有大量的NaN。可以使用.where()或.fillna()处理。
  2. 尝试使用ax.contourf的transform参数，并确保其与数据的实际坐标系一致（我们用的是ccrs.PlateCarree()）。
  3. 对于极地投影等变形大的区域，可以考虑先将数据插值到投影坐标系再绘图，但这更复杂。对于兰勃特投影，通常直接使用PlateCarree转换即可。

问题3：图形保存为PDF或SVG时，箭头或文字错位。

• 解决：这是一个常见的后端渲染问题。尝试在保存前调用plt.savefig(..., metadata={'Creator': ‘’, ‘Producer’: ‘’})，或者使用Agg后端（import matplotlib; matplotlib.use(‘Agg’)）在脚本开头非交互式地生成图形。

5.3 代码优化与封装建议

为了提升代码的复用性和可读性，建议将核心功能封装成函数或类：

class IVTCalculator: """整层水汽通量及散度计算器""" def __init__(self, data_path_pl, data_path_sl): self.ds_pl = xr.open_dataset(data_path_pl, chunks={'time': 1}) self.ds_sl = xr.open_dataset(data_path_sl) self.g = constants.g def calculate_ivt(self): """计算整层水汽通量""" # ... 集成上述计算步骤 ... return Q_u, Q_v def calculate_divergence(self, Q_u, Q_v): """计算水汽通量散度""" # ... 集成上述散度计算步骤 ... return div def plot_field(self, time_idx, variable='div', **plot_kwargs): """绘制指定变量场""" # ... 集成上述绘图步骤 ... pass # 使用示例 calc = IVTCalculator('era5_pl.nc', 'era5_sl.nc') Q_u, Q_v = calc.calculate_ivt() div = calc.calculate_divergence(Q_u, Q_v) calc.plot_field(0, variable='div', save_path='output.png')

将计算和绘图逻辑模块化，不仅使主程序清晰，也便于进行参数敏感性试验（如改变积分上限、平滑参数等）和批量处理多个天气个例。