基于MSWEP数据集计算全球坐标水文年降水总量问题求助

背景 #

处理合并后的NetCDF文件 precip_combined.nc，该文件包含1979-2020年全球月降水数据（3.79GB），xarray数据集结构如下：

<xarray.Dataset>
Dimensions:        (lon: 3600, lat: 1800, time: 504)
Coordinates:
  * lon            (lon) float32 -179.9 -179.8 -179.8 ... 179.8 179.9 179.9
  * lat            (lat) float32 89.95 89.85 89.75 ... -89.75 -89.85 -89.95
  * time           (time) datetime64[ns] 1979-01-01 1979-02-01 ... 2020-12-01
Data variables:
    precipitation  (time, lat, lon) float32 ...
Attributes:
    history:          Created on 2022-10-09 13:23
    input_data_hash:  1ff7f07166074240030172bae541afd3ecbc27ef941e50a65b04a26...

目标 #

1. 对每个经纬度坐标（lat, lon），计算水文年（当年10月至次年9月）的降水总和；
2. 将结果存入NetCDF文件，新增变量 precipitation_sums，保留经纬度维度，支持区域数据提取。

原代码问题 #

原代码中计算求和时错误指定了维度：

precipitation_sum = hydro_year_data['precipitation'].sum(dim=['lon', 'lat'])

这会对每个时间点的全球经纬度数据求和，最终得到仅含time维度的变量，丢失了lat和lon维度，无法提取特定区域数据。

解决方案 #

修正思路 #

核心是对时间维度求和，而非经纬度维度，保留lat和lon空间维度。同时注意原数据时间范围（1979-2020年），最后一个完整水文年是2019年10月-2020年9月，2020年10月-2021年9月的数据不完整，需排除。

方法1：修正循环逻辑 #

import xarray as xr
import pandas as pd

# 打开数据集
data = xr.open_dataset('.../precip_combined.nc')

start_month = 10
end_month = 9

precip_sums_list = []
hydro_year_labels = []

# 循环处理每个完整水文年（1979/1980 到 2019/2020）
for year in range(1979, 2020):
    start_date = f'{year}-{start_month:02d}-01'
    end_date = f'{year+1}-{end_month:02d}-01'
    # 筛选该水文年的降水数据
    hydro_data = data.sel(time=slice(start_date, end_date))
    # 对时间维度求和，保留lat和lon
    sum_hydro = hydro_data['precipitation'].sum(dim='time')
    precip_sums_list.append(sum_hydro)
    # 添加水文年标签（如"1979/1980"）
    hydro_year_labels.append(f"{year}/{year+1}")

# 合并结果，创建水文年维度
precip_sums = xr.concat(
    precip_sums_list,
    dim=pd.Index(hydro_year_labels, name='hydrological_year')
)

# 保存为新NetCDF文件
precip_sums.to_netcdf('.../precip_hydro_year_sums.nc')

方法2：使用groupby高效处理（推荐） #

无需循环，利用xarray的groupby功能批量处理：

import xarray as xr

data = xr.open_dataset('.../precip_combined.nc')

# 定义函数：给每个时间点分配对应的水文年
def assign_hydro_year(time):
    # 当年10月及以后属于下一年的水文年
    return time.year + 1 if time.month >= 10 else time.year

# 新增水文年坐标
data['hydro_year'] = xr.DataArray(
    data.time.map(assign_hydro_year),
    coords=data.time.coords,
    dims='time'
)

# 按水文年分组，对时间维度求和
precip_hydro_sums = data['precipitation'].groupby('hydro_year').sum(dim='time')

# 保存结果
precip_hydro_sums.to_netcdf('.../precip_hydro_year_sums.nc')

结果说明 #

生成的precip_hydro_year_sums.nc中，降水总和变量的维度为(hydro_year, lat, lon)，可直接提取任意区域（如伦敦）的水文年降水总和数据。

内容的提问来源于stack exchange，提问作者Cristallin