3D网格化Xarray数据集1D插值遇块大小为0问题求助

核心问题分析 #

你遇到的0维度数组问题，主要由以下几个错误导致：

1. 插值函数的返回位置错误

interp1d_chunk函数中，return new_chunk被放在了**内层循环（j循环）**里，导致函数只处理第一个经纬度网格点(i=0,j=0)就直接返回，不仅未完成所有网格的插值，还会让Dask无法正确识别输出数组的完整形状。

2. Dask map_blocks的参数传递错误

调用da.map_blocks时，直接传入了Xarray对象test_stemp，但map_blocks需要处理的是Dask数组而非Xarray DataArray；同时参数传递格式不符合要求，chunk=这类命名参数会被忽略，应直接传递数组作为位置参数。

3. 测试Xarray创建的语法错误

test_stemp的创建代码存在多处问题：

• xr.DataArray的参数格式错误，第二个参数直接写'lat':test_lats不符合函数签名
• 使用了未定义的变量lats、lat_var、lon_var，应替换为已定义的test_lats、'lat'、'lon'

修复后的完整代码 #

第一步：修正插值函数

将return new_chunk移到循环外部，确保所有经纬度网格都完成插值：

import numpy as np
import xarray as xr
import scipy.interpolate
import dask.array as da

def interp1d_chunk(chunk, new_depths, depths):
    nlat, nlon, _ = chunk.shape
    new_chunk = np.empty((nlat, nlon, len(new_depths)), dtype=chunk.dtype)
    
    for i in range(nlat):
        for j in range(nlon):
            f = scipy.interpolate.interp1d(
                depths, 
                chunk[i, j, :], 
                bounds_error=False, 
                fill_value="extrapolate"
            )
            new_chunk[i, j, :] = f(new_depths)
    return new_chunk

第二步：修正测试Xarray的创建

test_array = np.asarray([[[ 9.984375  ,  9.315826  ,  6.753296  , -0.71829224],
    [ 9.812378  ,  9.15155   ,  6.6003723 , -0.7545166 ],
    [ 9.925293  ,  9.266205  ,  6.729767  , -0.67056274],
    [ 9.925293  ,  9.266205  ,  6.729767  , -0.67056274]],

   [[10.201508  ,  9.524597  ,  6.9415283 , -0.6557312 ],
    [ 9.812378  ,  9.15155   ,  6.6003723 , -0.7545166 ],
    [10.083374  ,  9.394531  ,  6.7566833 , -0.7775879 ],
    [ 9.925293  ,  9.266205  ,  6.729767  , -0.67056274]],

   [[10.201508  ,  9.524597  ,  6.9415283 , -0.6557312 ],
    [ 9.812378  ,  9.15155   ,  6.6003723 , -0.7545166 ],
    [10.083374  ,  9.394531  ,  6.7566833 , -0.7775879 ],
    [ 9.925293  ,  9.266205  ,  6.729767  , -0.67056274]],

   [[10.201508  ,  9.524597  ,  6.9415283 , -0.6557312 ],
    [10.109436  ,  9.4236145 ,  6.80542   , -0.7413635 ],
    [10.083374  ,  9.394531  ,  6.7566833 , -0.7775879 ],
    [10.083374  ,  9.394531  ,  6.7566833 , -0.7775879 ]]],dtype=np.float32)

test_lats = [60.275, 60.325, 60.375, 60.425]
test_lons = [140.75, 140.8, 140.85, 140.9]
depths = np.asarray([3.5,17.5,64,194.5])
new_depths = np.arange(0,289.1,0.5)

# 修正后的Xarray创建
test_stemp = xr.DataArray(
    test_array,
    coords={
        'lat': test_lats,
        'lon': test_lons,
        'depth': depths
    },
    dims=['lat','lon','depth']
).rename('Tsoil').chunk(
    chunks={
        'lat': len(test_lats),
        'lon': len(test_lons),
        'depth': 4
    }
)

第三步：正确调用Dask map_blocks

传递底层的Dask数组，并指定正确的输出块形状：

# 获取Xarray的底层Dask数组
dask_arr = test_stemp.data

# 调用map_blocks
interp_dask_arr = da.map_blocks(
    interp1d_chunk,
    dask_arr,
    new_depths=new_depths,
    depths=depths,
    dtype=test_stemp.dtype,
    chunks=(len(test_lats), len(test_lons), len(new_depths))
)

# 将插值后的Dask数组转回Xarray
stemp_interp = xr.DataArray(
    interp_dask_arr,
    coords={
        'lat': test_stemp.lat,
        'lon': test_stemp.lon,
        'depth': new_depths
    },
    dims=['lat','lon','depth'],
    name='Tsoil'
)

# 验证结果形状
print(stemp_interp.shape)  # 输出应为(4,4,579)

额外优化建议 #

• 移除函数中的print语句，避免Dask任务执行时产生大量冗余输出
• 可以改用xr.apply_ufunc（Xarray通用函数接口），更贴合Xarray使用习惯，无需手动处理Dask块：

stemp_interp = xr.apply_ufunc(
    lambda x, old_d, new_d: scipy.interpolate.interp1d(old_d, x, axis=-1, bounds_error=False, fill_value="extrapolate")(new_d),
    test_stemp,
    depths,
    new_depths,
    input_core_dims=[['depth'], [], []],
    output_core_dims=[['depth']],
    output_sizes={'depth': len(new_depths)},
    vectorize=True,
    dask='parallelized',
    output_dtypes=[test_stemp.dtype]
)

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