1. WaterGAP全球水文模型(WGHM)简介WaterGAP全球水文模型(WGHM)是当前最权威的全球水资源评估工具之一由德国法兰克福大学开发维护。这个模型最厉害的地方在于它能模拟全球陆地水循环的各个关键环节包括降水、蒸发、径流、土壤湿度、地下水等水文要素。我最早接触这个模型是在2015年做全球水资源评估项目时当时就被它精细的空间分辨率(0.5°×0.5°)和长时间序列(1901年至今)所震撼。模型输出的数据格式主要是NetCDF这种格式特别适合存储多维科学数据。在实际应用中我们最常用的变量包括TWS (Terrestrial Water Storage)陆地水储量GW (Groundwater)地下水储量SM (Soil Moisture)土壤湿度ET (Evapotranspiration)蒸散发这些数据在气候变化研究、水资源管理、干旱监测等领域都有广泛应用。比如去年我们团队就用WGHM数据分析了中亚地区近20年地下水变化趋势发现了一些很有意思的规律。2. 数据获取与准备2.1 官方数据源定位找对数据源是第一步也是最容易踩坑的地方。根据我的经验最靠谱的数据来源是法兰克福大学官网和Pangaea地球科学数据库。这里分享一个实用技巧直接搜索WaterGAP v2.2d data download就能找到最新版本的数据集。数据下载时要注意几个关键点选择正确的驱动数据版本GSWP3-W5E5或WFDEI-GPCC确认时间范围1901-2016或2017-2019选择需要的水文变量TWS、GW等我建议第一次使用时下载watergap_22d_WFDEI-GPCC_histsoc_tws_monthly_1901_2016.nc4这个示例文件练手文件大小约1.5GB包含了完整的陆地水储量月数据。2.2 本地环境配置处理这种大尺寸的NetCDF文件我强烈推荐使用MATLABNetCDF库的组合。安装时要注意MATLAB版本最好在2018b以上必须安装Mapping Toolbox和NetCDF支持包如果是Linux系统需要提前安装libnetcdf-dev这里分享一个我常用的环境检查脚本% 检查必要工具箱是否安装 if ~license(test,map_toolbox) error(需要安装Mapping Toolbox); end if ~exist(ncread,file) error(需要安装NetCDF支持包); end3. 数据读取与预处理3.1 NetCDF文件解析读取NetCDF文件看似简单但有几个坑我不得不提醒。首先是坐标系问题WGHM数据采用0.5度网格经度从0.25°到359.75°纬度从-89.75°到89.75°。这个和常见的-180°到180°经度范围不同需要特别注意。下面是我优化过的数据读取代码file watergap_22d_WFDEI-GPCC_histsoc_tws_monthly_1901_2016.nc4; lon ncread(file,lon); % 0.25:0.5:359.75 lat ncread(file,lat); % -89.75:0.5:89.75 tws ncread(file,tws); % 720×360×1392 % 数据旋转和裁剪 tws_rot rot90(tws(:,:,1:12)); % 取1901年数据3.2 缺失值处理WGHM数据中的缺失值处理是个技术活。原始数据用-9999表示缺失值但直接替换为0可能会引入偏差。我的经验是采用空间插值法处理% 识别缺失值 missing_val -9999; tws_rot(tws_rot missing_val) NaN; % 简单空间插值 for i 2:size(tws_rot,1)-1 for j 2:size(tws_rot,2)-1 if isnan(tws_rot(i,j)) tws_rot(i,j) nanmean(tws_rot(i-1:i1,j-1:j1),all); end end end对于边缘区域可以考虑使用邻近格网的平均值或者直接保留NaN值具体取决于后续分析需求。4. 数据重构与分析4.1 时间序列处理WGHM数据的时间维度处理需要特别注意。原始数据是连续月数据但存储方式比较特殊。这里分享我整理的时间轴重建方法% 构建时间轴 start_date datetime(1901,1,1); end_date datetime(2016,12,31); date_vec start_date:calmonths(1):end_date; % 计算月度异常值 tws_monthly zeros(360,720,12); for m 1:12 monthly_idx month(date_vec) m; tws_monthly(:,:,m) mean(tws_rot(:,:,monthly_idx),3); end4.2 空间重采样有时候我们需要将数据重采样到其他分辨率比如1°×1°。这里给出一个保持水量守恒的重采样方法% 定义新网格 new_lon 0.5:1:359.5; new_lat -89.5:1:89.5; % 面积加权重采样 tws_resampled zeros(length(new_lat),length(new_lon)); for i 1:length(new_lat) for j 1:length(new_lon) lat_mask lat new_lat(i)-0.5 lat new_lat(i)0.5; lon_mask lon new_lon(j)-0.5 lon new_lon(j)0.5; tws_resampled(i,j) mean(tws_rot(lat_mask,lon_mask),all,omitnan); end end5. 可视化与成果输出5.1 使用GRACE工具箱制图冯伟老师的GRACE_Matlab_Toolbox确实是神器但使用时有几个技巧值得分享。首先是配色方案的选择我推荐使用jet或balance色标更适合表示水文异常值。% 加载工具箱 addpath(GRACE_Matlab_Toolbox-master); % 准备数据 temp squeeze(tws_monthly(:,:,1)); % 取1月数据 temp flipud(temp); % 转置并翻转 % 绘制全球图 figure(Position,[100 100 1200 600]); gmt_grid2map(temp*0.1, -30, 30, 0, 0, cm, WGHM TWS 1901.01); caxis([-50 50]); % 统一色标范围 colormap(jet(256)); colorbar(southoutside);5.2 自定义高级可视化对于发表级图表我通常会用更精细的控制参数。下面是一个制作多子图对比的示例% 创建1×2子图 figure(Position,[100 100 1400 600]); subplot(1,2,1); % 绘制1月数据 gmt_grid2map(squeeze(tws_monthly(:,:,1))*0.1, -30, 30, 0, 0, cm, January); caxis([-50 50]); subplot(1,2,2); % 绘制7月数据 gmt_grid2map(squeeze(tws_monthly(:,:,7))*0.1, -30, 30, 0, 0, cm, July); caxis([-50 50]); % 添加统一色标 cb colorbar(southoutside); set(cb,Position,[0.3 0.05 0.4 0.02]);6. 实际应用案例去年我们团队用WGHM数据做了一个有趣的研究分析全球主要流域2000-2016年的水储量变化趋势。这里分享部分关键代码% 加载流域边界 basin_mask ncread(global_basins.nc,mask); % 计算流域平均 unique_basins unique(basin_mask); basin_tws zeros(length(unique_basins),1392); for b 1:length(unique_basins) mask basin_mask unique_basins(b); for t 1:1392 tmp tws(:,:,t); basin_tws(b,t) mean(tmp(mask),omitnan); end end % 计算趋势 trends zeros(length(unique_basins),1); for b 1:length(unique_basins) p polyfit(1:1392,basin_tws(b,:),1); trends(b) p(1)*12; % 转换为年变化率 end这个分析揭示了印度恒河平原、华北平原等地区显著的地下水下降趋势成果最终发表在了水文领域顶级期刊上。7. 常见问题解决在使用WGHM数据过程中我遇到过各种奇怪的问题。这里整理几个最有代表性的内存不足错误处理全时间序列数据时1392个时间步很容易导致内存溢出。解决方案是分块处理chunk_size 120; % 10年数据 num_chunks ceil(1392/chunk_size); for c 1:num_chunks idx (c-1)*chunk_size1 : min(c*chunk_size,1392); chunk_data tws(:,:,idx); % 处理代码... end坐标偏移问题有时候绘图会发现大陆轮廓对不齐这通常是因为经度范围设置不当。正确的做法是% 调整经度范围 lon_adj lon; lon_adj(lon_adj 180) lon_adj(lon_adj 180) - 360;异常值处理除缺失值外有时会遇到异常大的正值或负值。我通常会用标准差过滤mean_val mean(tws_rot,all,omitnan); std_val std(tws_rot,0,all,omitnan); tws_rot(tws_rot mean_val3*std_val | tws_rot mean_val-3*std_val) NaN;8. 性能优化技巧处理全球尺度的水文数据效率很重要。经过多次优化我总结出几个提升速度的方法使用内存映射对于超大型NetCDF文件ncid netcdf.open(file,NOWRITE); varid netcdf.inqVarID(ncid,tws); tws netcdf.getVar(ncid,varid,double); netcdf.close(ncid);并行计算利用MATLAB的parfor加速循环parfor t 1:12 monthly_data(:,:,t) mean(tws(:,:,t:12:end),3); end预处理数据将常用数据保存为MAT文件save(WGHM_monthly.mat,monthly_data,-v7.3);这些技巧在处理长时间序列数据时可以将运行时间从几小时缩短到几分钟。