Python自动化处理GRACE-FO水文数据的完整指南

1. GRACE-FO数据概述与应用场景

GRACE-FO（Gravity Recovery and Climate Experiment Follow-On）是NASA与德国地球科学研究中心（GFZ）联合开展的卫星任务，通过监测地球重力场变化来追踪水循环过程。对于水文研究者而言，GRACE-FO提供的Level-3数据（如GFZ发布的NetCDF格式文件）能够反映陆地水储量（Terrestrial Water Storage, TWS）的时空变化，包括：

• 地下水储量变化
• 土壤湿度波动
• 积雪和冰川质量变化
• 地表水体（湖泊、河流）蓄水量

典型应用场景包括：

• 流域尺度水资源评估
• 干旱与洪水监测预警
• 冰川消融对海平面上升的贡献研究
• 跨区域地下水开采影响分析

# 示例：查看NetCDF文件结构
import netCDF4 as nc

ds = nc.Dataset('GFZ_GRACE-FO_TWS.nc')
print(ds.variables.keys())  # 输出：['time', 'lat', 'lon', 'tws', 'std_tws', ...]

2. 数据获取与预处理自动化

2.1 自动下载GRACE-FO数据

主流数据源包括：

• JPL Mascon数据：https://podaac-tools.jpl.nasa.gov/drive/files/allData/gracefo/L3/mascon/
• GFZ Level-3产品：ftp://rz-vm152.gfz-potsdam.de/gracefo/Level-3/
• CSR解决方案：https://www2.csr.utexas.edu/grace/RL06_mascons/

import ftplib
import os

def download_gracefo_ftp(remote_dir, local_dir):
    """通过FTP自动下载GRACE-FO数据"""
    ftp = ftplib.FTP('rz-vm152.gfz-potsdam.de')
    ftp.login()  # 匿名登录
    ftp.cwd(remote_dir)
    
    filenames = ftp.nlst('*.nc')  # 筛选NetCDF文件
    for filename in filenames:
        local_path = os.path.join(local_dir, filename)
        with open(local_path, 'wb') as f:
            ftp.retrbinary(f'RETR {filename}', f.write)
    ftp.quit()

提示：对于大规模下载，建议使用wget或curl工具配合断点续传功能，避免网络中断导致重复下载。

2.2 数据质量检查与清洗

原始数据常见问题处理：

问题类型	解决方法	Python实现
数据缺失	线性插值	scipy.interpolate.griddata
异常值	IQR过滤	numpy.percentile
坐标偏移	重投影	pyproj.transform
单位不一致	单位转换	pint库

import xarray as xr

def clean_gracefo_data(filepath):
    """数据清洗流程"""
    ds = xr.open_dataset(filepath)
    
    # 处理缺失值
    ds['tws'] = ds['tws'].interpolate_na(dim='time', method='linear')
    
    # 去除异常值（基于3σ原则）
    mean = ds['tws'].mean()
    std = ds['tws'].std()
    ds['tws'] = ds['tws'].where((ds['tws'] - mean).abs() < 3*std)
    
    return ds

3. 空间分析与区域掩膜应用

3.1 创建自定义流域掩膜

以长江流域为例，可通过以下步骤生成空间掩膜：

1. 从HydroSHEDS获取流域边界Shapefile
2. 使用GDAL将矢量转为栅格
3. 与GRACE-FO数据空间对齐

import geopandas as gpd
import rasterio
from rasterio import features

def create_basin_mask(shp_path, template_nc):
    """生成流域二值掩膜"""
    basin = gpd.read_file(shp_path).to_crs('EPSG:4326')
    with xr.open_dataset(template_nc) as ds:
        transform = rasterio.transform.from_bounds(
            west=ds.lon.min(), south=ds.lat.min(),
            east=ds.lon.max(), north=ds.lat.max(),
            width=len(ds.lon), height=len(ds.lat)
        )
        mask = features.geometry_mask(
            basin.geometry,
            out_shape=(len(ds.lat), len(ds.lon)),
            transform=transform,
            invert=True
        )
    return mask

3.2 区域统计分析

应用掩膜后，可计算流域内水储量变化的统计指标：

def basin_stats(tws_data, mask):
    """计算掩膜区域统计量"""
    masked_data = tws_data.where(mask)
    return {
        'mean': masked_data.mean(dim=('lat', 'lon')),
        'std': masked_data.std(dim=('lat', 'lon')),
        'trend': masked_data.polyfit(dim='time', deg=1)['polyfit_coefficients'][0]
    }

4. 时间序列分析与可视化

4.1 趋势分解与信号提取

GRACE-FO数据包含多种信号成分：

• 长期趋势（气候变化/人类活动）
• 季节周期（年际/半年际）
• 残余噪声（测量误差）

from statsmodels.tsa.seasonal import STL

def decompose_tws(ts):
    """时间序列分解"""
    res = STL(ts, period=12).fit()  # 假设月度数据
    return {
        'trend': res.trend,
        'seasonal': res.seasonal,
        'resid': res.resid
    }

4.2 动态可视化实现

使用matplotlib和cartopy创建专业级图表：

import matplotlib.pyplot as plt
import cartopy.crs as ccrs

def plot_tws_map(tws_slice, vmin=-20, vmax=20):
    """绘制TWS空间分布图"""
    fig = plt.figure(figsize=(12, 6))
    ax = fig.add_subplot(111, projection=ccrs.PlateCarree())
    
    # 绘制填色图
    mesh = ax.pcolormesh(tws_slice.lon, tws_slice.lat, tws_slice,
                         cmap='BrBG', vmin=vmin, vmax=vmax)
    
    # 添加地理要素
    ax.coastlines()
    ax.add_feature(cartopy.feature.BORDERS, linestyle=':')
    plt.colorbar(mesh, label='Equivalent Water Height (cm)')
    
    return fig

进阶技巧：创建交互式动态可视化

import plotly.express as px

def interactive_tws_animation(ds):
    """生成交互式时间动画"""
    fig = px.imshow(ds['tws'], 
                   animation_frame='time',
                   color_continuous_scale='BrBG',
                   range_color=[-30, 30],
                   labels={'color': 'TWS (cm)'})
    fig.update_layout(title='GRACE-FO Terrestrial Water Storage Animation')
    return fig

5. 完整数据处理流程示例

以下展示从原始数据到趋势分析的端到端流程：

# 1. 数据获取
download_gracefo_ftp('/gracefo/Level-3/GFZ/RL06/', './data')

# 2. 数据清洗
cleaned_ds = clean_gracefo_data('./data/GFZ_GRACE-FO_2020-2023.nc')

# 3. 区域分析
yangtze_mask = create_basin_mask('./boundaries/yangtze.shp', './data/GFZ_GRACE-FO_2020-2023.nc')
yangtze_tws = cleaned_ds['tws'].where(yangtze_mask)
stats = basin_stats(yangtze_tws, yangtze_mask)

# 4. 可视化
plot_tws_map(yangtze_tws.mean(dim='time'))
plt.savefig('yangtze_mean_tws.png')

# 5. 时间序列分析
ts_decomp = decompose_tws(stats['mean'])
ts_decomp['trend'].plot(title='Yangtze Basin TWS Trend')

注意：实际应用中应考虑数据间隙（如GRACE与GRACE-FO任务间的空档期）和误差校正（如冰川均衡调整GIA修正）。

通过上述方法，研究人员可以构建自动化流水线，将原始GRACE-FO数据转化为可直接用于科学研究的可视化成果。这种工作流程特别适合需要长期监测大区域水储量变化的团队，如流域管理机构或气候变化研究小组。