cartesian product of coordinates and using it to index / fill empty dataset

For a given empty dataset with only coordinates

import xarray as xr
import numpy as np
data = xr.Dataset(coords={'x': np.linspace(-1, 1), 'y': np.linspace(0, 10), 'a': 1, 'b': 5})

I'd like to iterate over the product of coordinates, in a similar way as it can be done for numpy.arrays

data = np.zeros((10, 5, 10))
for (i, j, k), _ in np.ndenumerate(data):
    data[i, j, k] = some_function(i, j, k)

to fill the data with values of some function.

Also I'd like to extend this to the cases of functions that are multi-valued, i.e. they return a numpy.array.

Is there an easy way to do so? I was unable to find anything similar in the docs.

Start with the empty xarray.Dataset example and compare its coordinate and indexing behavior with NumPy's np.ndenumerate example. Define the expected Cartesian-product iteration and handling of functions returning arrays, then verify that both scalar and multi-valued functions can fill the dataset as requested.