I had the chance to play around with stack and unstack, and it appears that these actually do nearly all the work needed here, so you can disregard my last comment. The only logic which is somewhat unwieldy is code which creates a DataArray from the eofs dask array. This is a complete example using the air dataset: ```python air = load_dataset("air_temperature").air

A = air.stack(features=['lat', 'lon']).chunk() A-= A.mean('features')

,,eofs = svd_compressed(A.data, 4)

wrap eofs in dataarray

dims = ['modes', 'features'] coords = {}

for i, dim in enumerate(dims): if dim in A.dims: coords[dim] = A[dim] elif dim in coords: pass else: coords[dim] = np.arange(eofs.shape[i])

eofs = xr.DataArray(eofs, dims=dims, coords=coords).unstack('features') `` This is pretty compact as is, so maybe the ugly final bit could be replaced with a convenience function likeunstack_array(eofs, dims, coords)or a method callA.unstack_array(eofs, dims, new_coords={})`.