用Python分析全球水资源变化：基于WaterGAP模型月数据（1901-2019）的完整流程-编程实验室

用Python分析全球水资源变化：基于WaterGAP模型月数据（1901-2019）的完整流程

当我们需要研究全球水资源的长期变化趋势时，WaterGAP模型提供的月尺度数据无疑是一个宝贵的资源。这份跨越119年的数据集，记录了从1901年到2019年间全球水循环的详细变化，包括地表水、地下水、土壤水等多种水储量指标。本文将带领读者从数据获取开始，一步步完成数据处理、分析和可视化的全过程，最终生成专业的水资源变化趋势图。

1. 环境准备与数据获取

在开始分析之前，我们需要搭建一个合适的工作环境。推荐使用Anaconda创建独立的Python环境，确保依赖包的版本一致性。以下是核心依赖包及其作用：

# 创建conda环境 conda create -n watergap python=3.9 conda activate watergap # 安装核心包 pip install xarray dask netCDF4 pandas numpy matplotlib cartopy

WaterGAP数据可以从Pangaea数据仓库获取。数据集通常以NetCDF格式存储，这种格式特别适合存储多维科学数据。我们可以使用xarray库高效地处理这些数据，它提供了类似pandas的接口，但专门为多维数组设计。

import xarray as xr # 示例数据加载 ds = xr.open_dataset('watergap_22d_gswp3-w5e5_histsoc_tws_monthly_1901_2019.nc4')

2. 数据预处理与质量控制

原始数据往往需要经过预处理才能用于分析。WaterGAP数据虽然已经过质量控制，但我们仍需进行一些基本检查和处理。

2.1 数据完整性检查

首先检查数据的时间覆盖范围和空间分辨率：

print(f"时间范围: {ds.time.min().values} 到 {ds.time.max().values}") print(f"空间分辨率: {ds.lon[1].values - ds.lon[0].values} 度") print(f"可用变量: {list(ds.data_vars)}")

2.2 缺失值处理

WaterGAP数据中的缺失值通常用特定值标记（如-9999），我们需要将其替换为NaN：

ds['tws'] = ds['tws'].where(ds['tws'] != -9999)

2.3 时间一致性检查

确保时间轴连续且无跳跃：

import pandas as pd # 检查时间间隔是否一致 time_diffs = pd.Series(ds.time.values[1:]) - pd.Series(ds.time.values[:-1]) print(f"时间间隔是否一致: {all(time_diffs == time_diffs[0])}")

3. 区域选择与时间序列分析

针对特定区域的分析是水资源研究的常见需求。下面以长江流域为例，展示如何提取区域数据并进行分析。

3.1 定义区域边界

长江流域的大致经纬度范围：

yangtze_bbox = { 'lon_min': 90, 'lon_max': 122, 'lat_min': 24, 'lat_max': 35 }

3.2 区域数据提取

使用xarray的sel方法提取区域数据：

yangtze_ds = ds.sel( lon=slice(yangtze_bbox['lon_min'], yangtze_bbox['lon_max']), lat=slice(yangtze_bbox['lat_max'], yangtze_bbox['lat_min']) )

3.3 时间序列聚合

计算区域平均时间序列：

yangtze_ts = yangtze_ds['tws'].mean(dim=['lon', 'lat'])

4. 趋势分析与可视化

长期趋势分析是水资源研究的核心内容。下面介绍几种常用的分析方法。

4.1 年际变化趋势

首先将月数据聚合为年平均值：

yearly_mean = yangtze_ts.groupby('time.year').mean()

使用线性回归计算趋势：

from scipy.stats import linregress years = yearly_mean.year.values values = yearly_mean.values slope, intercept, r_value, p_value, std_err = linregress(years, values) print(f"趋势: {slope*10:.2f} mm/10年, p值: {p_value:.4f}")

4.2 可视化展示

使用matplotlib绘制时间序列和趋势线：

import matplotlib.pyplot as plt import matplotlib.dates as mdates plt.figure(figsize=(12, 6)) plt.plot(yearly_mean.year, yearly_mean, label='年平均值') plt.plot(years, intercept + slope*years, 'r', label='趋势线') plt.title('长江流域总水储量年际变化 (1901-2019)') plt.xlabel('年份') plt.ylabel('总水储量 (mm)') plt.legend() plt.grid() plt.show()

5. 空间分布与变化模式

除了时间序列分析，空间分布特征也是理解水资源变化的重要方面。

5.1 多年平均空间分布

计算整个时期的平均水储量：

mean_tws = ds['tws'].mean(dim='time')

5.2 变化趋势的空间分布

计算每个格点的线性趋势：

from scipy.stats import linregress def calc_trend(ts): years = np.arange(len(ts)) slope, _, _, _, _ = linregress(years, ts) return slope * len(years) # 总变化量 trend = xr.apply_ufunc( calc_trend, ds['tws'].chunk({'lat': 10, 'lon': 10}), input_core_dims=[['time']], output_core_dims=[[]], vectorize=True )

5.3 空间可视化

使用cartopy绘制空间分布图：

import cartopy.crs as ccrs import cartopy.feature as cfeature proj = ccrs.PlateCarree() fig = plt.figure(figsize=(15, 8)) ax = fig.add_subplot(111, projection=proj) # 添加地理要素 ax.add_feature(cfeature.LAND) ax.add_feature(cfeature.COASTLINE) ax.add_feature(cfeature.BORDERS, linestyle=':') # 绘制趋势 im = trend.plot(ax=ax, transform=proj, cmap='coolwarm', cbar_kwargs={'label': '水储量变化 (mm/世纪)'}) ax.set_title('全球总水储量变化趋势 (1901-2019)') plt.show()

6. 多变量分析与相关性研究

水资源系统各分量之间存在复杂的相互作用，多变量分析有助于理解这些关系。

6.1 数据整合

加载多个相关变量：

variables = ['tws', 'groundwstor', 'soilmoist', 'swe'] datasets = {} for var in variables: filename = f'watergap_22d_gswp3-w5e5_histsoc_{var}_monthly_1901_2019.nc4' datasets[var] = xr.open_dataset(filename)[var]

6.2 区域平均时间序列

计算长江流域各变量的年平均值：

yearly_means = {} for var, da in datasets.items(): region_da = da.sel( lon=slice(yangtze_bbox['lon_min'], yangtze_bbox['lon_max']), lat=slice(yangtze_bbox['lat_max'], yangtze_bbox['lat_min']) ) yearly_means[var] = region_da.groupby('time.year').mean()

6.3 相关性分析

计算各变量间的相关系数：

import pandas as pd df = pd.DataFrame(yearly_means) corr_matrix = df.corr() print("变量间相关系数矩阵:") print(corr_matrix)

7. 高级分析与应用

在前面的基础上，我们可以进行更深入的分析，为水资源管理提供科学依据。

7.1 干旱事件识别

定义基于TWS的干旱指标：

# 计算标准化异常 tws_mean = yearly_means['tws'].mean() tws_std = yearly_means['tws'].std() tws_anomaly = (yearly_means['tws'] - tws_mean) / tws_std # 识别严重干旱年份 drought_years = tws_anomaly[tws_anomaly < -1.5] print(f"严重干旱年份: {list(drought_years.index)}")

7.2 变化点检测

使用Pettitt检验检测突变点：

from pyhomogeneity import pettitt_test result = pettitt_test(yearly_means['tws'].values) print(f"突变点年份: {years[result.cp]}, 显著性: {result.p:.3f}")

7.3 未来情景预测

基于历史趋势的简单预测：

future_years = np.arange(2020, 2051) future_tws = intercept + slope * future_years plt.figure(figsize=(10, 5)) plt.plot(years, yearly_means['tws'], label='观测') plt.plot(future_years, future_tws, '--', label='预测') plt.title('长江流域总水储量预测') plt.xlabel('年份') plt.ylabel('总水储量 (mm)') plt.legend() plt.grid() plt.show()

8. 性能优化与大数据处理

当处理全球长时间序列数据时，性能优化至关重要。以下是几种有效的优化策略。

8.1 使用Dask进行并行计算

import dask.array as da # 分块加载数据 ds = xr.open_mfdataset('watergap_*.nc4', chunks={'time': 120, 'lat': 100, 'lon': 100}) # 并行计算 mean_tws = ds['tws'].mean(dim='time').compute()

8.2 内存优化技巧

对于大型计算，可以采取以下策略：

使用xarray的chunk方法控制内存使用
及时删除不再需要的中间变量
使用dask的persist方法缓存常用数据

# 示例：分块计算 chunked = ds['tws'].chunk({'time': 120, 'lat': 50, 'lon': 50}) result = chunked.groupby('time.year').mean().compute()

8.3 结果存储优化

将中间结果保存为Zarr格式，提高后续读取效率：

# 保存为Zarr格式 ds.to_zarr('watergap_tws.zarr') # 从Zarr读取 ds_zarr = xr.open_zarr('watergap_tws.zarr')

用Python分析全球水资源变化：基于WaterGAP模型月数据（1901-2019）的完整流程