用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')