使用 dask.array.map_overlap 时的 dask 输出问题

Question

我想用dask.array.map_overlap来处理scipy插值函数。但是，我总是遇到无法理解的错误，希望有人能回答我。

这是我收到的错误消息，如果我想运行 .compute().

ValueError: could not broadcast input array from shape (1070,0) into shape (1045,0)

为了解决这个问题，我开始使用.to_delayed()来检查每个分区的输出，这是我发现的。

以下是我的 python 代码。

步骤1.通过Xarray加载netCDF文件，然后输出到dask.array，块大小为(400,400)

df = xr.open_dataset('./Brazil Sentinal2 Tile/' + data_file +'.nc')
lon, lat = df['lon'].data, df['lat'].data
slon = da.from_array(df['lon'], chunks=(400,400))
slat = da.from_array(df['lat'], chunks=(400,400))
data = da.from_array(df.isel(band=0).__xarray_dataarray_variable__.data, chunks=(400,400))

第 2 步。为 da.map_overlap use

声明一个函数

def sumsum2(lon,lat,data,  hex_res=10):
    hex_col = 'hex' + str(hex_res)
    lon_max, lon_min = lon.max(), lon.min()
    lat_max, lat_min = lat.max(), lat.min()
    
    b = box(lon_min, lat_min, lon_max, lat_max, ccw=True)
    b = transform(lambda x, y: (y, x), b)
    b = mapping(b)
    
    target_df = pd.DataFrame(h3.polyfill( b, hex_res), columns=[hex_col])    

    target_df['lat'] = target_df[hex_col].apply(lambda x: h3.h3_to_geo(x)[0])
    target_df['lon'] = target_df[hex_col].apply(lambda x: h3.h3_to_geo(x)[1])
    tlon, tlat = target_df[['lon','lat']].values.T    

    abc = lNDI(points=(lon.ravel(), lat.ravel()), 
               values= data.ravel())(tlon,tlat)
    target_df['out'] = abc
    print(np.stack([tlon, tlat, abc],axis=1).shape)
    return np.stack([tlon, tlat, abc],axis=1)

第 3 步。应用 da.map_overlap

b = da.map_overlap(sumsum2, slon[:1200,:1200], slat[:1200,:1200], data[:1200,:1200], depth=10, trim=True, boundary=None, align_arrays=False, dtype='float64', 
                  )

第 4 步。使用 to_delayed() 测试输出形状

print(b.to_delayed().flatten()[0].compute().shape, )
print(b.to_delayed().flatten()[1].compute().shape)

(1065, 3)
(1045, 0)
(1090, 3)
(1070, 0)

表示 da.map_overlap 的输出仅输出一维维度（即 (1045,0) 和 (1070,0) ），而在 da.map_overlap 中，我准备的输出是二维维度（即 (1065,3) 和 (1090,3)）。

此外，如果我关闭 trim 参数，即

c = da.map_overlap(sumsum2, 
                   slon[:1200,:1200], 
                   slat[:1200,:1200], 
                   data[:1200,:1200], 
                   depth=10,
                   trim=False,
                   boundary=None,
                   align_arrays=False,
                   dtype='float64', 
                  )

print(c.to_delayed().flatten()[0].compute().shape, )
print(c.to_delayed().flatten()[1].compute().shape)

输出变为

(1065, 3)
(1065, 3)
(1090, 3)
(1090, 3)

这是说trim=True的时候我把所有东西都剪掉了？

因为...

#-- print out the values 
b.to_delayed().flatten()[0].compute()[:10,:]

(1065, 3)
array([], shape=(1045, 0), dtype=float64)

而...

#-- print out the values
c.to_delayed().flatten()[0].compute()[:10,:]

array([[ -47.83683837, -18.98359832, 1395.01848583],
[ -47.8482856 , -18.99038681, 2663.68391094],
[ -47.82800624, -18.99207069, 1465.56517187],
[ -47.81897323, -18.97919009, 2769.91556363],
[ -47.82066663, -19.00712956, 1607.85927095],
[ -47.82696896, -18.97167714, 2110.7516765 ],
[ -47.81562653, -18.98302933, 2662.72112163],
[ -47.82176881, -18.98594465, 2201.83205114],
[ -47.84567 , -18.97512514, 1283.20631652],
[ -47.84343568, -18.97270783, 1282.92117225]])

对此有什么想法吗？

谢谢。

Answer 1

我想我得到了答案。有错请指教

我不允许使用 trim=True 是因为我改变了输出数组的形状（上网后，我注意到输出数组的形状应该与输入数组的形状相同）。由于我改变了形状，dask 不知道如何处理它所以它 return 对我来说是空数组（奇怪）。
而不是使用trim=False，因为我没有问cutting-out缓冲区，所以现在可以了输出 return 值。（虽然我仍然不知道为什么dask不能连接分块数组，但相信也与形状有关）
解决方案是在da.concatenate上使用delayed函数，即

delayed(da.concatenate)([e.to_delayed().flatten()[idx] for idx in range(len(e.to_delayed().flatten()))])

在这种情况下，我们不依赖 map_overlap 中的 concat 函数，而是使用我们自己的 concat 来组合我们想要的输出。

使用 dask.array.map_overlap 时的 dask 输出问题

question for dask output when using dask.array.map_overlap

dask