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为例访问平台通过ECMWF的CDSClimate Data StoreAPI进行下载。你需要先在官网注册账号并获取API密钥。确定要素计算水汽通量我们需要以下变量specific_humidity比湿 qu_component_of_wind纬向风 uv_component_of_wind经向风 vsurface_pressure地表气压用于确定积分下界时空范围与层次时间选择你需要分析的日期和时间例如一次暴雨过程期间。ERA5提供逐小时数据。空间选定经纬度范围。为了计算水平散度区域范围应比最终成图范围稍大一圈避免边界缺值。层次选择气压层数据。通常需要从高层如100 hPa到近地面如1000 hPa的多层数据。ERA5提供了从1000 hPa到1 hPa的多个标准层。对于水汽通量积分积分到100 hPa通常已足够因为高层水汽含量极低。使用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 python3.9 conda activate meteorology # 安装核心库 conda install -c conda-forge xarray dask netCDF4 cartopy matplotlib numpy scipyxarray是处理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。关键预处理步骤单位统一ERA5的比湿单位是kg/kg风分量单位是m/s地表气压单位是Pa。这些单位与公式要求一致通常无需转换。但务必确认。坐标对齐确保两个数据集气压层和单层的timelatitudelongitude坐标完全一致。xarray的运算会自动对齐坐标但为保险起见可以用.sel或.interp进行精确匹配。处理缺失值再分析数据通常质量很高但沿海或复杂地形处可能有缺测值。可以用.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, dtypefloat) levs pressure_levels.values # 内部层使用梯形法的权重 (p_{i1} - p_{i-1}) / 2 for i in range(1, len(levs)-1): dp_values[i] (levs[i-1] - levs[i1]) / 2.0 # 顶层i0假设其上界气压为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(levelpressure) # 重新索引确保对应关系正确 # 将地表气压扩展到与integrand相同的维度结构除了level # 这里surface_pressure需要与integrand的‘time’ ‘latitude’ ‘longitude’对齐 sp_expanded surface_pressure.broadcast_like(integrand.isel(level0)) # 创建一个掩膜标识出气压值大于地表气压的层次即在地表以下的层次不参与积分 # 因为气压向下增加所以“大于”地表气压的层次是无效的。 mask pressure sp_expanded # 将无效层次的被积函数设为0 integrand_masked integrand.where(~mask, 0) # 对level维进行梯形积分 # xarray的integrate方法可以方便地进行梯形积分 integral integrand_masked.integrate(coordlevel) 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, radius6371000.): 计算球面上矢量场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_longitude105, central_latitude35, standard_parallels(25, 47)) ax fig.add_subplot(1, 1, 1, projectionproj) # 设置地图范围西经 东经 南纬 北纬 map_extent [70, 140, 15, 55] ax.set_extent(map_extent, crsccrs.PlateCarree()) # 注意set_extent的坐标范围需用PlateCarree经纬度表示 # 添加地理特征 ax.add_feature(cfeature.COASTLINE.with_scale(50m), linewidth0.8) ax.add_feature(cfeature.BORDERS.with_scale(50m), linewidth0.5, linestyle:) ax.add_feature(cfeature.OCEAN, colorlightcyan, alpha0.6) ax.add_feature(cfeature.LAND, colorwheat, alpha0.3) # 添加省界需要中国省界shapefile这里用河流替代示意实际应用中可加载自定义数据 ax.add_feature(cfeature.RIVERS.with_scale(50m), colorblue, linewidth0.5, alpha0.5) # 添加经纬网格线 gl ax.gridlines(draw_labelsTrue, dmsTrue, x_inlineFalse, y_inlineFalse, linewidth0.5, colorgray, alpha0.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(timetime_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, levelslevels, cmapRdBu_r, transformccrs.PlateCarree(), extendboth) # extendboth 表示色标向两端延伸 # 添加色标 cbar plt.colorbar(cf, axax, orientationhorizontal, pad0.05, aspect40, shrink0.8) cbar.set_label(Integrated Water Vapor Flux Divergence (kg m$^{-2}$ s$^{-1}$), fontsize12)4.3 叠加水汽通量矢量箭头风羽水汽通量矢量Q_u, Q_v用箭头风羽表示可以直观看到水汽输送的方向和强度。由于通量矢量通常很大直接画箭头会过于密集需要降采样。# 降采样每隔几个格点画一个箭头 stride 5 # 根据你的格点密度调整密度大则stride大一些 Q_u_plot Q_u.isel(timetime_idx) Q_v_plot Q_v.isel(timetime_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, # 用速度着色箭头 cmapplasma, scalescale, scale_unitsinches, transformccrs.PlateCarree(), width0.002, headwidth3, headlength4) # 为风羽图添加一个独立的色标可选表示通量强度 # cbar2 plt.colorbar(quiver, axax, orientationvertical, pad0.1, shrink0.8) # cbar2.set_label(IVT Magnitude (kg m$^{-1}$ s$^{-1}$), fontsize12)4.4 添加标题与修饰# 添加标题 plot_time pd.to_datetime(str(div_plot.time.values)).strftime(%Y-%m-%d %H:%M UTC) ax.set_title(fIntegrated Vapor Transport and Divergence\n{plot_time}, fontsize16, fontweightbold, pad20) # 添加文本标注如研究区域或说明 ax.text(0.02, 0.98, Blue: Convergence (Moisture Sink)\nRed: Divergence (Moisture Source), transformax.transAxes, verticalalignmenttop, bboxdict(boxstyleround, facecolorwheat, alpha0.8), fontsize10) plt.tight_layout() # 保存图像 plt.savefig(IVT_Divergence.png, dpi300, bbox_inchestight) plt.show()实操心得绘图的调试往往比计算更耗时。关键点在于1)投影选择要适合你的分析区域2)颜色映射要符合气象学惯例如散度用红蓝降水用蓝白红3)箭头缩放scale参数需要多次尝试才能达到疏密适中、长度合适的效果4) 图形要素海岸线、网格线、标签要清晰但不喧宾夺主。建议将绘图代码封装成函数便于对不同时次或个例进行批量出图。5. 常见问题、排查技巧与性能优化在实际操作中你几乎一定会遇到各种报错和不如预期的结果。这里汇总了一些典型问题及其解决方案。5.1 计算相关的问题问题1垂直积分结果量级异常大或小。可能原因单位错误检查比湿q的单位是否为kg/kg无量纲。ERA5数据通常是正确的。检查重力加速度g的单位是m/s²。积分上下限错误确认积分是从高压地面向低压高空进行。检查pressure_levels数组是否按降序排列。确保surface_pressure参与积分下界的判断。地表气压处理不当在高原地区如果直接用1000 hPa作为积分下界会严重高估水汽通量。务必使用integrate_vertical函数中实现的、基于真实地表气压的掩膜方法。排查方法打印出单个格点例如一个海洋点和一个高原点的pressure_levels、sp值以及积分过程中被积函数和掩膜的情况进行人工核对。问题2散度场出现棋盘状噪声或极端值。可能原因微分对噪声敏感原始风场和湿度场本身有小尺度噪声微分会将其放大。网格非均匀虽然ERA5是规则网格但如果你使用的数据是变分辨率或跳点采样的直接使用中心差分公式会不准。边界效应在计算区域边界中心差分缺少一侧的格点可能导致异常值。解决方案平滑滤波在计算散度前对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 ... # 同理使用更稳健的差分方法可以考虑使用二次精度的差分格式或者直接调用气象专用库如windspharm来计算球面散度后者经过充分测试更为可靠。剔除边界计算散度后将边界附近一圈格点的值设为NaN。问题3内存不足处理大数据时卡死。原因ERA5高时空分辨率数据体积庞大直接读入内存可能超过限制。解决方案使用Dask进行分块计算xarray与dask无缝集成。在打开数据集时使用chunks参数。ds_pl xr.open_dataset(era5_data_pl.nc, chunks{time: 1, level: 5, latitude: 100, longitude: 100})这样数据以“惰性”方式加载计算任务会被图优化并分块执行极大减少内存峰值。分时次处理如果不需要做时间平均可以循环处理每个时次处理完一个就保存或绘图然后释放内存。降低空间分辨率对于大范围气候分析可以先用.coarsen()或.interp()方法将数据插值到更粗的网格上。5.2 绘图相关的问题问题1箭头quiver太密集或太稀疏看不清。解决调整stride参数和scale参数。stride控制跳过的格点数scale控制箭头的整体长度。两者配合调整。一个经验法则是让最强的箭头长度大约等于图上5个经度/纬度的距离。问题2填色图contourf在Cartopy投影上出现奇怪的空洞或扭曲。原因数据可能包含NaN值或者地图投影的转换在数据边界处出现问题。解决确保绘图数据在绘图区域内没有大量的NaN。可以使用.where()或.fillna()处理。尝试使用ax.contourf的transform参数并确保其与数据的实际坐标系一致我们用的是ccrs.PlateCarree()。对于极地投影等变形大的区域可以考虑先将数据插值到投影坐标系再绘图但这更复杂。对于兰勃特投影通常直接使用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, variablediv, **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, variablediv, save_pathoutput.png)将计算和绘图逻辑模块化不仅使主程序清晰也便于进行参数敏感性试验如改变积分上限、平滑参数等和批量处理多个天气个例。 Python实战：基于ERA5数据计算与可视化整层水汽通量及散度
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为例访问平台通过ECMWF的CDSClimate Data StoreAPI进行下载。你需要先在官网注册账号并获取API密钥。确定要素计算水汽通量我们需要以下变量specific_humidity比湿 qu_component_of_wind纬向风 uv_component_of_wind经向风 vsurface_pressure地表气压用于确定积分下界时空范围与层次时间选择你需要分析的日期和时间例如一次暴雨过程期间。ERA5提供逐小时数据。空间选定经纬度范围。为了计算水平散度区域范围应比最终成图范围稍大一圈避免边界缺值。层次选择气压层数据。通常需要从高层如100 hPa到近地面如1000 hPa的多层数据。ERA5提供了从1000 hPa到1 hPa的多个标准层。对于水汽通量积分积分到100 hPa通常已足够因为高层水汽含量极低。使用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 python3.9 conda activate meteorology # 安装核心库 conda install -c conda-forge xarray dask netCDF4 cartopy matplotlib numpy scipyxarray是处理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。关键预处理步骤单位统一ERA5的比湿单位是kg/kg风分量单位是m/s地表气压单位是Pa。这些单位与公式要求一致通常无需转换。但务必确认。坐标对齐确保两个数据集气压层和单层的timelatitudelongitude坐标完全一致。xarray的运算会自动对齐坐标但为保险起见可以用.sel或.interp进行精确匹配。处理缺失值再分析数据通常质量很高但沿海或复杂地形处可能有缺测值。可以用.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, dtypefloat) levs pressure_levels.values # 内部层使用梯形法的权重 (p_{i1} - p_{i-1}) / 2 for i in range(1, len(levs)-1): dp_values[i] (levs[i-1] - levs[i1]) / 2.0 # 顶层i0假设其上界气压为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(levelpressure) # 重新索引确保对应关系正确 # 将地表气压扩展到与integrand相同的维度结构除了level # 这里surface_pressure需要与integrand的‘time’ ‘latitude’ ‘longitude’对齐 sp_expanded surface_pressure.broadcast_like(integrand.isel(level0)) # 创建一个掩膜标识出气压值大于地表气压的层次即在地表以下的层次不参与积分 # 因为气压向下增加所以“大于”地表气压的层次是无效的。 mask pressure sp_expanded # 将无效层次的被积函数设为0 integrand_masked integrand.where(~mask, 0) # 对level维进行梯形积分 # xarray的integrate方法可以方便地进行梯形积分 integral integrand_masked.integrate(coordlevel) 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, radius6371000.): 计算球面上矢量场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_longitude105, central_latitude35, standard_parallels(25, 47)) ax fig.add_subplot(1, 1, 1, projectionproj) # 设置地图范围西经 东经 南纬 北纬 map_extent [70, 140, 15, 55] ax.set_extent(map_extent, crsccrs.PlateCarree()) # 注意set_extent的坐标范围需用PlateCarree经纬度表示 # 添加地理特征 ax.add_feature(cfeature.COASTLINE.with_scale(50m), linewidth0.8) ax.add_feature(cfeature.BORDERS.with_scale(50m), linewidth0.5, linestyle:) ax.add_feature(cfeature.OCEAN, colorlightcyan, alpha0.6) ax.add_feature(cfeature.LAND, colorwheat, alpha0.3) # 添加省界需要中国省界shapefile这里用河流替代示意实际应用中可加载自定义数据 ax.add_feature(cfeature.RIVERS.with_scale(50m), colorblue, linewidth0.5, alpha0.5) # 添加经纬网格线 gl ax.gridlines(draw_labelsTrue, dmsTrue, x_inlineFalse, y_inlineFalse, linewidth0.5, colorgray, alpha0.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(timetime_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, levelslevels, cmapRdBu_r, transformccrs.PlateCarree(), extendboth) # extendboth 表示色标向两端延伸 # 添加色标 cbar plt.colorbar(cf, axax, orientationhorizontal, pad0.05, aspect40, shrink0.8) cbar.set_label(Integrated Water Vapor Flux Divergence (kg m$^{-2}$ s$^{-1}$), fontsize12)4.3 叠加水汽通量矢量箭头风羽水汽通量矢量Q_u, Q_v用箭头风羽表示可以直观看到水汽输送的方向和强度。由于通量矢量通常很大直接画箭头会过于密集需要降采样。# 降采样每隔几个格点画一个箭头 stride 5 # 根据你的格点密度调整密度大则stride大一些 Q_u_plot Q_u.isel(timetime_idx) Q_v_plot Q_v.isel(timetime_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, # 用速度着色箭头 cmapplasma, scalescale, scale_unitsinches, transformccrs.PlateCarree(), width0.002, headwidth3, headlength4) # 为风羽图添加一个独立的色标可选表示通量强度 # cbar2 plt.colorbar(quiver, axax, orientationvertical, pad0.1, shrink0.8) # cbar2.set_label(IVT Magnitude (kg m$^{-1}$ s$^{-1}$), fontsize12)4.4 添加标题与修饰# 添加标题 plot_time pd.to_datetime(str(div_plot.time.values)).strftime(%Y-%m-%d %H:%M UTC) ax.set_title(fIntegrated Vapor Transport and Divergence\n{plot_time}, fontsize16, fontweightbold, pad20) # 添加文本标注如研究区域或说明 ax.text(0.02, 0.98, Blue: Convergence (Moisture Sink)\nRed: Divergence (Moisture Source), transformax.transAxes, verticalalignmenttop, bboxdict(boxstyleround, facecolorwheat, alpha0.8), fontsize10) plt.tight_layout() # 保存图像 plt.savefig(IVT_Divergence.png, dpi300, bbox_inchestight) plt.show()实操心得绘图的调试往往比计算更耗时。关键点在于1)投影选择要适合你的分析区域2)颜色映射要符合气象学惯例如散度用红蓝降水用蓝白红3)箭头缩放scale参数需要多次尝试才能达到疏密适中、长度合适的效果4) 图形要素海岸线、网格线、标签要清晰但不喧宾夺主。建议将绘图代码封装成函数便于对不同时次或个例进行批量出图。5. 常见问题、排查技巧与性能优化在实际操作中你几乎一定会遇到各种报错和不如预期的结果。这里汇总了一些典型问题及其解决方案。5.1 计算相关的问题问题1垂直积分结果量级异常大或小。可能原因单位错误检查比湿q的单位是否为kg/kg无量纲。ERA5数据通常是正确的。检查重力加速度g的单位是m/s²。积分上下限错误确认积分是从高压地面向低压高空进行。检查pressure_levels数组是否按降序排列。确保surface_pressure参与积分下界的判断。地表气压处理不当在高原地区如果直接用1000 hPa作为积分下界会严重高估水汽通量。务必使用integrate_vertical函数中实现的、基于真实地表气压的掩膜方法。排查方法打印出单个格点例如一个海洋点和一个高原点的pressure_levels、sp值以及积分过程中被积函数和掩膜的情况进行人工核对。问题2散度场出现棋盘状噪声或极端值。可能原因微分对噪声敏感原始风场和湿度场本身有小尺度噪声微分会将其放大。网格非均匀虽然ERA5是规则网格但如果你使用的数据是变分辨率或跳点采样的直接使用中心差分公式会不准。边界效应在计算区域边界中心差分缺少一侧的格点可能导致异常值。解决方案平滑滤波在计算散度前对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 ... # 同理使用更稳健的差分方法可以考虑使用二次精度的差分格式或者直接调用气象专用库如windspharm来计算球面散度后者经过充分测试更为可靠。剔除边界计算散度后将边界附近一圈格点的值设为NaN。问题3内存不足处理大数据时卡死。原因ERA5高时空分辨率数据体积庞大直接读入内存可能超过限制。解决方案使用Dask进行分块计算xarray与dask无缝集成。在打开数据集时使用chunks参数。ds_pl xr.open_dataset(era5_data_pl.nc, chunks{time: 1, level: 5, latitude: 100, longitude: 100})这样数据以“惰性”方式加载计算任务会被图优化并分块执行极大减少内存峰值。分时次处理如果不需要做时间平均可以循环处理每个时次处理完一个就保存或绘图然后释放内存。降低空间分辨率对于大范围气候分析可以先用.coarsen()或.interp()方法将数据插值到更粗的网格上。5.2 绘图相关的问题问题1箭头quiver太密集或太稀疏看不清。解决调整stride参数和scale参数。stride控制跳过的格点数scale控制箭头的整体长度。两者配合调整。一个经验法则是让最强的箭头长度大约等于图上5个经度/纬度的距离。问题2填色图contourf在Cartopy投影上出现奇怪的空洞或扭曲。原因数据可能包含NaN值或者地图投影的转换在数据边界处出现问题。解决确保绘图数据在绘图区域内没有大量的NaN。可以使用.where()或.fillna()处理。尝试使用ax.contourf的transform参数并确保其与数据的实际坐标系一致我们用的是ccrs.PlateCarree()。对于极地投影等变形大的区域可以考虑先将数据插值到投影坐标系再绘图但这更复杂。对于兰勃特投影通常直接使用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, variablediv, **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, variablediv, save_pathoutput.png)将计算和绘图逻辑模块化不仅使主程序清晰也便于进行参数敏感性试验如改变积分上限、平滑参数等和批量处理多个天气个例。