如何对内存超限的大型地理空间xarray数据集进行降采样？

我有一个尺寸为86400×43200的大型地理空间xarray数据集，包含经纬度（lat, lon）坐标与高程数据，希望将其分辨率降至原维度的1/20（即4320×2160，对应1/12度）。

最初直接使用interp方法时，该方法会加载整个数据集，导致需要28GB内存，超出我的内存分配能力。

当前采用的方法步骤如下：

• 创建填充为0的降采样新数据集
• 将大型数据集切分为36个均等切片
• 加载并插值单个切片
• 将降采样后的切片添加到新数据集

但结果全为NaN值，步骤2可正常执行，仅步骤3存在问题。请问如何正确对这类大型数据集进行降采样，或正确填充新数据集？

附上当前代码：

import xarray as xr
import numpy as np  # 原代码遗漏的导入

def format_spatial_resolution(
    xarray, res_lat=1 / 12, res_lon=1 / 12, sample_array=None
):

    new_lon = np.arange(xarray.lon[0], xarray.lon[-1], res_lon)
    new_lat = np.arange(xarray.lat[0], xarray.lat[-1], res_lat)
    xarray_new = xarray.interp(
        lat=new_lat,
        lon=new_lon,
        # 外插以避免边缘出现NaN值
        kwargs={
            "fill_value": "extrapolate",
        },
    )
    return xarray_new


xarray = xr.open_dataset(filename)
# 设置分辨率
res_lat = 1/12 # 1/12度
res_lon = 1/12
# 创建填充为0的降采样xarray
new_lon = np.arange(xarray.lon[0], xarray.lon[-1], res_lon)
new_lat = np.arange(xarray.lat[0], xarray.lat[-1], res_lat)

xarray_downsampled = xr.Dataset(
    dict(elevation=(["lat", "lon"], np.zeros((len(new_lat), len(new_lon))))),  # 修正原代码括号错误
    coords=dict(lon=new_lon, lat=new_lat),
)
# 插值每个切片并更新降采样xarray
for i in range(-180, 180, 10):
    xarray_downsampled.update(format_spatial_resolution(xarray.sel(lon=slice(i, i + 10))))

解决方案 #

1. 当前代码的核心问题 #

使用xarray_downsampled.update()是错误的：该方法会用新切片的插值结果完全覆盖原有变量，后续循环会不断替换前面的结果，最终只有最后一个切片的有效数据，其余区域因未被正确赋值而变成NaN（原0值会被覆盖）。

正确做法是通过坐标索引匹配，将每个切片的插值结果精准赋值到目标数据集的对应位置。

2. 修正后的代码 #

import xarray as xr
import numpy as np

def format_spatial_resolution(xarray, res_lat=1/12, res_lon=1/12):
    # 直接使用全局目标坐标进行插值，确保结果维度匹配
    return xarray.interp(
        lat=new_lat,
        lon=new_lon,
        kwargs={"fill_value": "extrapolate"}
    )

# 延迟加载原数据集，不一次性读入内存
xarray = xr.open_dataset(filename)
res_lat = 1/12
res_lon = 1/12

# 生成全局目标分辨率的经纬度坐标
new_lon = np.arange(xarray.lon[0], xarray.lon[-1], res_lon)
new_lat = np.arange(xarray.lat[0], xarray.lat[-1], res_lat)

# 创建NaN初始化的目标数据集，后续赋值覆盖有效区域
xarray_downsampled = xr.Dataset(
    dict(elevation=(["lat", "lon"], np.full((len(new_lat), len(new_lon)), np.nan))),
    coords=dict(lon=new_lon, lat=new_lat),
)

# 循环处理每个经度切片
for i in range(-180, 180, 10):
    # 提取当前经度范围的切片
    slice_ds = xarray.sel(lon=slice(i, i + 10))
    # 插值到目标分辨率
    interpolated_slice = format_spatial_resolution(slice_ds)
    # 将插值结果赋值到目标数据集的对应经度区域
    xarray_downsampled["elevation"].loc[dict(lon=slice(i, i + 10))] = interpolated_slice["elevation"]

# 保存降采样结果
xarray_downsampled.to_netcdf("downsampled_elevation.nc")

3. 更高效的降采样方案 #

如果你的数据是规则网格，优先使用coarsen或interp_like，比手动切片更高效且内存占用更低：

• coarsen聚合降采样：适合整数倍降采样（你的场景正好是20倍），直接对网格做均值/最大值等聚合：

# 按20倍因子聚合经纬度，自动处理分块避免内存溢出
downsampled = xarray.coarsen(lon=20, lat=20, boundary="trim").mean()
downsampled.to_netcdf("coarsened_elevation.nc")

• interp_like分块插值：利用xarray延迟计算自动分块，无需手动切片：

# 创建目标分辨率的模板数据集
template = xr.Dataset(coords=dict(lon=new_lon, lat=new_lat))
# 自动分块插值，内存占用可控
downsampled = xarray.interp_like(template, kwargs={"fill_value": "extrapolate"}).compute()

4. 额外注意事项 #

• 确保原数据集经纬度坐标连续且规则，否则插值/聚合可能出现异常。
• 处理大型数据集时，保留xarray默认的延迟加载特性，避免一次性读入全量数据。
• 若切片后仍有内存压力，可给sel添加chunks参数启用分块加载：

slice_ds = xarray.sel(lon=slice(i, i + 10)).chunk({"lat": 1000, "lon": 1000})

内容的提问来源于stack exchange，提问作者a-doering