thanks a lot @fujiisoup, your suggestion does help getting rid of the necessity to build the ds['_y'] variable. Here is the updated apply_ufunc solution: ``` x = np.arange(10100) y = np.arange(20100)

ds = xr.Dataset(coords={'x': x, 'y': y})

ds = ds.chunk({'x': 1, 'y':1}) # does not change anything

let's say each experiment outputs 5 statistical diagnostics

Nstats = 5 some_exp = lambda x, y: np.ones((Nstats,))

out = xr.apply_ufunc(some_exp, ds.x, ds.y, dask='parallelized', vectorize=True, output_dtypes=[float], output_sizes={'stats': Nstats}, output_core_dims=[['stats']]) ``` An inspection of the dask dashboard indicates that the computation is not distributed among workers though. How could I make sure this happens?