1. 项目概述从“水汽收支”到一张图做气象、水文或者地理研究的朋友对“区域水汽收支”这个概念肯定不陌生。简单来说它就像给一个特定的地理区域比如一个流域、一个省份甚至整个青藏高原算一笔“水汽账”在一定时间内有多少水汽通过大气流动平流和蒸发进入这个区域又有多少通过降水和流出离开。算清楚这笔账是理解区域水候变化、极端降水事件、干旱成因乃至水资源评估的物理基础。“区域净水汽收支计算并绘图”这个项目标题听起来很学术但拆解开来就是一个非常典型的“数据获取-科学计算-可视化呈现”的科研或业务分析流程。它的核心目标是把抽象的、多维的大气物理量通过严谨的计算最终浓缩成一张或一套直观的、信息量丰富的图表。这张图可能就是一篇论文的核心证据或一份气候评估报告的关键结论。我处理过不少类似的项目从全球尺度到城市群尺度都有。实话说这个过程里80%的精力可能都花在了数据预处理和计算逻辑的校准上真正画图反而是最后水到渠成的一步。但正是前面的“脏活累活”决定了最终结果的可靠性和说服力。接下来我就以一个典型的研究场景——比如计算“长江中下游地区夏季的净水汽收支”——为例把整个流程掰开揉碎了讲清楚包括每一步的原理、实操中的“坑”以及如何高效地得到一张漂亮的成果图。2. 核心原理与数据基础算的是什么用什么算在动手写代码之前我们必须彻底搞清楚要计算的对象到底是什么以及需要哪些数据来支撑这个计算。这能避免后续很多无效劳动。2.1 水汽收支方程拆解大气中的水汽收支通常基于一个“气柱”模型来考虑。对于一个给定的区域其大气柱在单位时间内的净水汽收支Net Moisture Budget, NMB可以表示为NMB (E - P) ∇·Q这里每个符号都有明确的物理意义E: 蒸发Evaporation。从地表海洋、湖泊、土壤、植被进入大气的水汽。P: 降水Precipitation。从大气中降落到地面的水分。∇·Q: 水汽通量散度Moisture Flux Divergence。这是最关键的一项代表了由于风场携带水汽进出该区域所造成的水汽净收入或支出。Q是整层积分的水汽通量矢量Q (1/g) ∫ qV dp其中 q 是比湿V 是水平风矢量g 是重力加速度积分从地面到大气顶。注意这个公式是“诊断”公式。在实际应用中我们通常从再分析资料或模式输出中直接获取 E、P 以及各层的 q、u纬向风、v经向风然后计算 ∇·Q。理论上如果数据完美公式左右应该平衡NMB ≈ 0但由于观测误差、再分析资料的同化过程以及计算本身的近似总会存在一个残差。这个残差的大小也是检验我们计算过程和数据质量的一个参考。2.2 数据源选择再分析资料的“门派之争”计算水汽收支我们几乎不可能自己去布设全球探空站最现实、最常用的数据源是大气再分析资料。这就好比做菜食材数据的好坏直接决定成品结果的档次。目前主流的选择有几个ERA5欧洲中期天气预报中心当前绝对的“顶流”。时空分辨率高0.25°×0.25°逐小时变量齐全且不断更新。对于区域尺度研究ERA5通常是首选。它的数据一致性很好官网提供的计算工具CDS Toolbox也能部分简化工作。MERRA-2美国宇航局同样非常优秀尤其在气溶胶和水文循环方面有特色。它的“水文学”变量输出可能更直接符合我们的需求。JRA-55日本气象厅时间序列长质量稳定在亚洲区域表现良好是长期气候趋势分析的常用数据集。NCEP/NCAR Reanalysis 1老牌经典数据集分辨率较粗约2.5°但时间序列很长1948年至今适合做初步探索或对计算资源要求不高的长时间序列分析。我的选型心得对于2000年以后的、需要精细刻画区域过程的研究无脑推荐ERA5。它的数据易获取性、社区支持度和精度综合得分最高。如果研究涉及更早的年代比如1979年前可能需要考虑JRA-55或NCEP。在做对比研究或不确定性分析时同时使用2-3套再分析资料是非常好的做法。2.3 关键变量清单确定了数据源比如ERA5接下来就要明确需要下载哪些变量。以月平均数据为例计算季节或年尺度收支常用核心变量包括地表变量单层evaporation或e(ERA5中名为evaporation)total_precipitation或tp可选runoff,surface_pressure用于交叉验证或更精细计算。高空变量多层specific_humidity(q)u_component_of_wind(u)v_component_of_wind(v)geopotential用于计算气压层厚度或位势高度通常需要从1000 hPa到至少100 hPa的多层数据。数据获取实操提示通过ECMWF的CDS APIPython库cdsapi下载ERA5数据是最佳途径。你需要先注册获取API key。下载时务必注意时间范围精确到年月日。区域area下载比你目标区域稍大一圈的数据因为后续计算散度需要空间差分边界会损失一圈格点。例如研究区域是[20°N-35°N, 110°E-125°E]下载时可以设定为[18°N-37°N, 108°E-127°E]。格式NetCDF是标准选择兼容性好。3. 计算流程全解析从格点数据到收支数字拿到NetCDF数据文件后真正的计算就开始了。这个过程可以分解为几个清晰的步骤。3.1 数据读取与预处理我习惯用xarray库来处理NetCDF数据它对于这种多维网格数据简直是“神器”。import xarray as xr import numpy as np import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature # 读取数据 ds_q xr.open_dataset(specific_humidity.nc) ds_u xr.open_dataset(u_wind.nc) ds_v xr.open_dataset(v_wind.nc) ds_evap xr.open_dataset(evaporation.nc) ds_precip xr.open_dataset(precipitation.nc) # 确保所有数据的时间、层次、经纬度坐标对齐 # 通常再分析资料坐标一致但保险起见可以检查 print(ds_q.coords)预处理的关键点单位统一ERA5的比湿q单位是kg/kg风场u/v是m/s蒸发和降水单位是m液态水当量深度但需要注意其累积量属性。月平均的tp和e通常是“月累积量”单位为米。计算收支时我们通常需要转换为通量形式即单位时间单位面积的水量例如 kg m⁻² s⁻¹ 或 mm day⁻¹。时间处理如果下载了逐月数据可能需要按季节如JJA代表夏季或年份进行平均。区域裁剪使用sel方法裁剪到目标区域。# 定义目标区域 lat_range [20, 35] lon_range [110, 125] # 裁剪数据 (假设经度是0-360格式需要转换到-180-180或直接使用) # ERA5经度默认是0-360我们可以这样处理 def convert_lon(ds): ds ds.assign_coords(longitude(((ds.longitude 180) % 360) - 180)) ds ds.sortby(longitude) return ds ds_q convert_lon(ds_q) ds_q_region ds_q.sel(latitudeslice(lat_range[1], lat_range[0]), # 注意纬度降序 longitudeslice(lon_range[0], lon_range[1])) # 对其他变量做同样处理...3.2 核心计算水汽通量散度这是整个项目最核心、也最容易出错的部分。计算整层积分的水汽通量散度 ∇·Q。步骤1计算各层的水汽通量在气压坐标下单位气柱在某一层的水汽通量向量为(qu/g, qv/g)其中除以重力加速度g是为了将质量单位统一。更常见的做法是直接计算q*u和q*v单位kg kg⁻¹ m s⁻¹然后在垂直积分时考虑气压差。步骤2垂直积分我们需要将各层的水汽通量从地面到大气顶进行积分。在离散的气压层数据上这相当于一个加权求和Q Σ (q * V * Δp / g)对所有气压层求和。 其中 Δp 是两层之间的气压差单位Pa。g 取 9.80665 m/s²。# 假设我们已经裁剪并对齐了 ds_q, ds_u, ds_v它们都有维度 (time, level, lat, lon) # 获取气压层数据 (单位Pa)ERA5的level坐标就是hPa需要转换 pressure ds_q.level * 100.0 # 转换为Pa # 计算气压厚度 Δp注意边界处理。这里采用中心差分顶层和底层用单边差分。 dp np.zeros_like(pressure.values) dp[0] pressure[0] - pressure[1] # 顶层 dp[-1] pressure[-2] - pressure[-1] # 底层 for i in range(1, len(pressure)-1): dp[i] (pressure[i-1] - pressure[i1]) / 2.0 # 将dp转换为xarray DataArray便于广播计算 dp_da xr.DataArray(dp, dims[level], coords{level: ds_q.level}) # 计算整层积分的水汽通量 Q_u 和 Q_v g 9.80665 Q_u (ds_q * ds_u * dp_da / g).sum(dimlevel) # 对level维度求和 Q_v (ds_q * ds_v * dp_da / g).sum(dimlevel) # 现在 Q_u, Q_v 的维度是 (time, lat, lon)单位是 kg m⁻¹ s⁻¹步骤3计算水平散度在球坐标系下水平散度的计算公式为∇·Q 1/(a cosφ) * [∂(Q_λ)/∂λ ∂(Q_φ cosφ)/∂φ]其中 a 是地球半径约6371 kmφ 是纬度弧度λ 是经度弧度。Q_λ 是纬向通量Q_uQ_φ 是经向通量Q_v。我们可以利用xarray的差分方法.differentiate或者手动使用np.gradient来计算偏导数。这里使用np.gradient更直观# 获取经纬度并转换为弧度 lat np.deg2rad(Q_u.latitude.values) lon np.deg2rad(Q_u.longitude.values) a 6371000.0 # 地球半径单位米 # 创建经纬度网格 lon_grid, lat_grid np.meshgrid(lon, lat) # 计算经纬向的格距弧度 dlon np.gradient(lon_grid, axis1) # 经度方向差分 dlat np.gradient(lat_grid, axis0) # 纬度方向差分 # 将Q_u, Q_v转换为numpy数组以便计算假设是单个时次 Q_u_np Q_u.isel(time0).values # 取第一个时次 Q_v_np Q_v.isel(time0).values # 计算偏导数 dQu_dlon np.gradient(Q_u_np, axis1) / dlon # 注意∂(Q_φ cosφ)/∂φ 的计算 Qv_coslat Q_v_np * np.cos(lat_grid) dQvcos_dlat np.gradient(Qv_coslat, axis0) / dlat # 计算散度 divQ (1/(a * np.cos(lat_grid))) * (dQu_dlon dQvcos_dlat) # 单位: kg m⁻² s⁻¹ # 将结果存回xarray DataArray divQ_da xr.DataArray(divQ, dims[latitude, longitude], coords{latitude: Q_u.latitude, longitude: Q_u.longitude})重要提醒散度计算对边界很敏感。上述方法在区域边界会产生不可靠的值。因此我们通常先计算全局或大区域的散度场然后再裁剪出目标区域或者直接忽略边界一圈格点。3.3 净收支计算与单位转换有了散度场divQ以及区域的蒸发E和降水P我们就可以计算每个格点的净水汽收支变化率∂W/∂t E - P - ∇·Q其中∂W/∂t代表大气柱水汽含量的局地变化率通常对于月或更长时间尺度平均这项接近于零因此有E - P ≈ ∇·Q。我们更关心的是区域的面积积分平均值。计算区域平均对目标区域内的所有格点的(E-P)和∇·Q进行面积加权平均。# 计算每个格点的面积权重球面梯形面积 lat_rad np.deg2rad(divQ_da.latitude) dlat np.gradient(lat_rad) # 纬度间隔弧度 dlon np.deg2rad(2.5) # 假设经度分辨率为2.5度需根据实际数据调整 area_weight a**2 * np.cos(lat_rad[:, np.newaxis]) * dlat[:, np.newaxis] * dlon area_weight_da xr.DataArray(area_weight, dims[latitude, longitude], coords{latitude: divQ_da.latitude, longitude: divQ_da.longitude}) # 读取并裁剪E和P数据并转换为相同单位如 mm/day # ERA5的月平均蒸发和降水是累积量m需要转换为通量 # 假设 ds_evap 和 ds_precip 是月累积量m/month days_in_month 30.44 # 近似值最好根据具体月份计算 E_flux ds_evap.evaporation * 1000 / days_in_month # 转换为 mm/day P_flux ds_precip.tp * 1000 / days_in_month # 转换为 mm/day E_minus_P E_flux - P_flux # 计算区域面积加权平均 def area_weighted_mean(field): return (field * area_weight_da).sum(dim[latitude, longitude]) / area_weight_da.sum() mean_E_minus_P area_weighted_mean(E_minus_P_region) mean_divQ area_weighted_mean(divQ_da) # divQ单位是 kg m⁻² s⁻¹ # 将 divQ 单位转换为 mm/day 以便比较1 kg m⁻² s⁻¹ 86400 mm/day mean_divQ_mmday mean_divQ * 86400计算净收支对于整个区域净水汽收入正值表示水汽净流入为Net Moisture Inflow - mean_divQ_mmday因为散度为正表示净流出。 理论上它应该约等于mean_E_minus_P。两者的差值残差可以用来评估计算闭合度。4. 可视化绘图让数据自己说话计算出一堆数字和格点场后最终要通过图表呈现。好的可视化能瞬间传达核心信息。4.1 空间分布图水汽通量散度场这是最直观的图可以展示水汽在哪里汇聚负散度降水潜势大、在哪里辐散正散度可能对应干旱区。import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature # 创建地图投影 fig plt.figure(figsize(12, 8)) ax plt.axes(projectionccrs.PlateCarree()) ax.set_extent([lon_range[0], lon_range[1], lat_range[0], lat_range[1]], crsccrs.PlateCarree()) # 添加地理特征 ax.add_feature(cfeature.COASTLINE, linewidth0.8) ax.add_feature(cfeature.BORDERS, linewidth0.5, linestyle:) ax.add_feature(cfeature.RIVERS, linewidth0.5, edgecolorblue, alpha0.5) ax.add_feature(cfeature.LAKES, linewidth0.5, edgecolorblue, facecolornone) # 绘制填色图 # divQ_da 是之前计算的水汽通量散度场单位已转换为 mm/day im ax.pcolormesh(divQ_da.longitude, divQ_da.latitude, divQ_da.values * 86400, # 假设之前是kg m-2 s-1 cmapRdBu_r, vmin-5, vmax5, transformccrs.PlateCarree()) # 添加等值线可选 # cs ax.contour(divQ_da.longitude, divQ_da.latitude, divQ_da.values*86400, # levelsnp.arange(-4, 4.5, 1), colorsk, linewidths0.5, transformccrs.PlateCarree()) # ax.clabel(cs, inlineTrue, fontsize9, fmt%1.0f) # 添加色标 cbar plt.colorbar(im, axax, orientationhorizontal, pad0.05, aspect40, shrink0.8) cbar.set_label(Moisture Flux Divergence (mm day$^{-1}$), fontsize12) # 负值蓝色表示水汽汇聚净流入正值红色表示水汽辐散净流出 # 添加标题和网格 ax.set_title(f夏季平均水汽通量散度 (∇·Q)\n{region_name}, fontsize14, fontweightbold) ax.gridlines(draw_labelsTrue, dmsTrue, x_inlineFalse, y_inlineFalse, linestyle--, alpha0.5) plt.tight_layout() plt.savefig(moisture_flux_divergence.png, dpi300, bbox_inchestight) plt.show()绘图技巧色标选择RdBu_r红蓝反转是经典选择蓝色冷色通常表示负值水汽汇聚/流入红色暖色表示正值辐散/流出符合直觉。叠加矢量箭头可以叠加平均的整层水汽通量矢量(Q_u, Q_v)的箭头直观显示水汽输送的主流路径。使用ax.quiver注意需要对矢量场进行稀疏化采样否则箭头会太密。突出研究区域用矩形框或阴影突出显示你计算平均值的具体子区域。4.2 时间序列图收支分量的年际变化如果你计算了多年数据绘制各分量E, P, ∇·Q的时间序列或气候态月变化图非常有价值。# 假设我们已经计算了1981-2020年每年夏季JJA的区域平均 mean_E, mean_P, mean_divQ_mmday years np.arange(1981, 2021) mean_E ... # 形状为 (40,) 的数组 mean_P ... # 形状为 (40,) 的数组 mean_divQ ... # 形状为 (40,) 的数组单位 mm/day已取负号表示净流入 fig, axes plt.subplots(2, 1, figsize(14, 10)) # 子图1各分量年际变化 ax1 axes[0] ax1.bar(years, mean_E, width0.6, labelEvaporation (E), colorskyblue, edgecolorblack) ax1.bar(years, mean_P, width0.6, labelPrecipitation (P), bottommean_E, colorlightcoral, edgecolorblack, alpha0.7) ax1.plot(years, mean_divQ, k-o, linewidth2, markersize5, label-∇·Q (Net Inflow)) ax1.axhline(y0, colorgrey, linestyle-, linewidth0.5) ax1.set_ylabel(Water Flux (mm day$^{-1}$), fontsize12) ax1.set_title(f{region_name} 夏季水汽收支分量年际变化 (1981-2020), fontsize14, fontweightbold) ax1.legend(locupper left, fontsize11) ax1.grid(True, alpha0.3) # 子图2净收支E-P与 -∇·Q 的对比及残差 ax2 axes[1] net_balance mean_E - mean_P residual net_balance - mean_divQ # 残差 (E-P) - (-∇·Q) E-P∇·Q理论上应为0 width 0.35 x np.arange(len(years)) ax2.bar(x - width/2, net_balance, width, labelE - P, colordarkorange, edgecolorblack) ax2.bar(x width/2, mean_divQ, width, label-∇·Q (Net Inflow), colorforestgreen, edgecolorblack) ax2.plot(x, residual, s-, colorpurple, labelResidual (E-P∇·Q), linewidth1.5) ax2.axhline(y0, colorgrey, linestyle-, linewidth0.5) ax2.set_xticks(x[::5]) # 每5年显示一个刻度 ax2.set_xticklabels(years[::5]) ax2.set_ylabel(Water Flux (mm day$^{-1}$), fontsize12) ax2.set_xlabel(Year, fontsize12) ax2.set_title(净水汽收支平衡对比与残差, fontsize13) ax2.legend(locupper right, fontsize11) ax2.grid(True, alpha0.3) plt.tight_layout() plt.savefig(moisture_budget_timeseries.png, dpi300, bbox_inchestight) plt.show()这种对比图能清晰展示蒸发、降水、净水汽输送之间的平衡关系以及计算闭合度残差大小。理想情况下残差线应该在0附近小幅波动。4.3 综合剖面或箱线图对于更深入的分析你还可以绘制垂直剖面图展示水汽通量散度或水汽输送的垂直结构。箱线图比较不同子区域如上、中、下游或不同年代际水汽收支统计特征的差异。空间趋势图用Sen‘s斜率或Mann-Kendall检验计算 ∇·Q 场的空间变化趋势并绘图。5. 常见问题、排查技巧与心得这个项目看似流程清晰但实操中会遇到各种“坑”。下面是我总结的一些典型问题和解决方法。5.1 数据与计算问题问题1计算出的散度场量级异常大或空间图案杂乱无章。可能原因1单位混乱。这是最常见的问题。确保q(kg/kg),u/v(m/s),dp(Pa),g(m/s²) 单位正确。检查垂直积分后的Q_u/Q_v单位是否为 kg m⁻¹ s⁻¹。最终散度单位应为 kg m⁻² s⁻¹转换为 mm/day 需要乘以 86400。可能原因2经纬度单位未转换为弧度。np.gradient计算偏导数时分母必须是弧度差。务必使用np.deg2rad()转换。可能原因3风场或比湿数据存在缺测值NaN。在计算q*u前用.fillna(0)或插值方法处理缺测值但需理解其气象合理性。排查方法先在一个格点上进行手算验证。选取一个点提取其上下几层的 q, u, v, p手动计算该点的垂直积分和散度与程序结果对比。问题2区域平均的 (E-P) 与 -∇·Q 相差很大残差过大。可能原因1数据时间不匹配。确保蒸发、降水、高空变量是同一时间段、同一时间分辨率如都是月平均的数据。可能原因2区域平均未做面积加权。高纬度格点实际面积小直接算术平均会引入偏差。必须使用cos(纬度)进行面积加权。可能原因3再分析资料本身的闭合误差。不同再分析资料的水文循环闭合程度不同。ERA5相对较好但仍有误差。可以计算多年平均的残差如果存在系统性偏差可能是资料本身的特性。可能原因4垂直积分不完整。计算∇·Q时积分上限应达到足够高的气压层如50 hPa或100 hPa确保包含了绝大部分水汽约99%的水汽在对流层。检查你下载的数据是否包含了所有必要层次。问题3绘图时地图投影扭曲或数据错位。可能原因Cartopy 投影与数据坐标不匹配。确保绘图时transformccrs.PlateCarree()参数正确设置且数据的经纬度坐标是十进制度。如果数据是0-360°经度务必先转换为-180°-180°。5.2 效率与优化技巧分块计算与并行如果处理全球多年数据内存和计算量巨大。可以使用dask与xarray结合进行分块和并行计算。在打开数据集时使用chunks参数。ds xr.open_mfdataset(era5_*.nc, chunks{time: 10, level: 5, latitude: 100, longitude: 100})预计算与存储垂直积分Q_u/Q_v的计算比较耗时。可以将其作为中间结果计算出来并保存为NetCDF文件后续计算散度和绘图时直接读取避免重复计算。使用气候学工具库对于标准计算可以考虑使用像windspharm计算球谐函数散度、xMIP、climpred等专业气候诊断工具包它们内置了经过验证的散度、涡度等计算函数更可靠高效。5.3 科学解读要点算出结果并画出漂亮的图只是第一步更重要的是正确解读符号意义∇·Q 0表示水汽通量辐散该地区大气是水汽的“源”通常对应净蒸发大于降水 (E P)。∇·Q 0表示水汽通量辐合该地区大气是水汽的“汇”通常对应净降水大于蒸发 (P E)。长期变化分析净水汽收支时间序列的趋势可以揭示区域干湿状况的长期变化并与大型环流指数如季风指数、ENSO指数做相关分析探究其驱动机制。空间格局结合地形、海陆分布、主要天气系统如副高、急流、低涡来解释水汽通量散度的空间分布特征。最后分享一个我个人的深刻体会永远不要相信第一次跑出来的结果。一定要用多种方式进行验证比如用不同的再分析资料ERA5 vs MERRA-2做对比用更简单的“盒子法”计算区域四边界的净水汽通量输入输出来验证散度法的结果或者与已有的权威文献中的图表进行定性对比。只有当多个独立的数据源和方法指向一致的结论时你的计算结果才是坚实可信的。这个过程虽然繁琐但正是科研工作从“做出图表”到“获得认知”的关键一跃。