issue_comments: 417802225

I started having a go at writing the second half of this - the "n-dimensional-concatenation" function which would accept a grid of xarray.DataSet/DataArray objects (assumed to be in the correct order along all dimensions), and return a single merged dataset.

However, I think there's an issue with using

something like a NumPy object array of xarray.Dataset/DataArray objects.

My plan was to call np.apply_along_axis to apply the 1D xr.concat() function along each axis in turn, something like ```python from numpy import apply_along_axis from xarray import concat

def concat_nd(obj_grid, concat_dims=None): """ Concatenates a structured ndarray of xarray Datasets along multiple dimensions.

Parameters
----------
obj_grid : numpy array of Dataset and DataArray objects
    N-dimensional numpy object array containing xarray objects in the shape they
    are to be concatenated. Each object is expected to
    consist of variables and coordinates with matching shapes except for
    along the concatenated dimension.
concat_dims : list of str or DataArray or pandas.Index
    Names of the dimensions to concatenate along. Each dimension in this argument
    is passed on to :py:func:`xarray.concat` along with the dataset objects.
    Should therefore be a list of valid dimension arguments to xarray.concat().

Returns
-------
combined : xarray.Dataset
"""

# Combine datasets along one dimension at a time
# Start with last axis and finish with axis=0
for axis in reversed(range(obj_grid.ndim)):
    obj_grid = apply_along_axis(concat, axis, arr=obj_grid, dim=concat_dims[axis])

# Grid should now only contain one xarray object
return obj_grid.item

However, testing this code withpython def test_concat_1d(self): data = create_test_data()

split_data = [data.isel(dim1=slice(3)), data.isel(dim1=slice(3, None))]

# Will explain why I'm forced to create ndarray like this shortly
split_data_grid = np.empty(shape=(2), dtype=np.object)
split_data_grid[0] = split_data[0]
split_data_grid[1] = split_data[1]

reconstructed = concat_nd(split_data_grid, ['dim1'])

xrt.assert_identical(data, reconstructed)

throws an error from within `np.apply_along_axis` TypeError: cannot directly convert an xarray.Dataset into a numpy array. Instead, create an xarray.DataArray first, either with indexing on the Dataset or by invoking the to_array() method. ```

I think this is because even just the idea of having a ndarray containing xarray datasets seems to cause problems - if I do it with a single item then xarray thinks I'm trying to convert the Dataset into a numpy array and throws the same error: python data = create_test_data() data_grid = np.array(data, dtype='object')

and if I do it with multiple items then numpy will dive down and extract the variables in the dataset instead of just storing a reference to the dataset: python data = create_test_data() split_data = [data.isel(dim1=slice(3)), data.isel(dim1=slice(3, None))] split_data_grid = np.array(split_data, dtype='object') print(split_data_grid) returns [['time' 'dim2' 'dim3' 'var1' 'var2' 'var3' 'numbers'] ['time' 'dim2' 'dim3' 'var1' 'var2' 'var3' 'numbers']] when I expected something more like numpy.array([<xarray.Dataset object at 0x7f94efa62278>, <xarray.Dataset object at 0x9f23phj10582>]) (This is why I had to create an empty array and then fill it afterwards in my example test further up.)

Is this the intended behaviour of xarray? Does this mean I can't use numpy arrays of xarray objects at all for this problem? If so then what structure do you think I should use instead (list of lists etc.)?