Unstack non-index coordinate

What happened:

Non-indexed coordinates are not unstacked correctly. The example below should be self-explanatory.

In the example below, I create a surface of the torus which should be used for some calculation for each point - therefore I have to use multi-index. The surface is given by the set of (R, Z) -- indexed by pol_idx -- which is the same for each toroidal angle phi. At the end of the day I would like to have the mesh of ('phi', 'pol_idx').

What you expected to happen:

ds.stack(multiindex=("a", "b")).unstack("multindex")should returnds`.

Minimal Complete Verifiable Example:

theta = np.linspace(0, np.pi)
R0 = 1

grid = xr.Dataset(coords={"R": ("pol_idx", R0 * np.cos(theta)), 
                          "Z": ("pol_idx", R0 * np.sin(theta)), 
                          "phi": np.linspace(0, 2*np.pi, N_phi_points, endpoint=False)})
> <xarray.Dataset>
> Dimensions:  (phi: 256, pol_idx: 50)
> Coordinates:
>    R        (pol_idx) float64 1.0 0.9979 0.9918 0.9816 ... -0.9918 -0.9979 -1.0
>    Z        (pol_idx) float64 0.0 0.06407 0.1279 ... 0.1279 0.06407 1.225e-16
>  * phi      (phi) float64 0.0 0.02454 0.04909 0.07363 ... 6.21 6.234 6.259
> Dimensions without coordinates: pol_idx
> Data variables:
>   *empty*

grid = grid.stack(index=("pol_idx", "phi"))

> <xarray.Dataset>
> Dimensions:  (index: 12800)
> Coordinates:
>     R        (index) float64 1.0 1.0 1.0 1.0 1.0 ... -1.0 -1.0 -1.0 -1.0 -1.0
>     Z        (index) float64 0.0 0.0 0.0 0.0 ... 1.225e-16 1.225e-16 1.225e-16
>   * index    (index) MultiIndex
>   - pol_idx  (index) int64 0 0 0 0 0 0 0 0 0 0 ... 49 49 49 49 49 49 49 49 49 49
>   - phi      (index) float64 0.0 0.02454 0.04909 0.07363 ... 6.21 6.234 6.259
> Data variables:
>     *empty*

# Here I would like to perform some calculation on the prepared grid over the multiindex.

grid = grid.unstack("index")

> <xarray.Dataset>
> Dimensions:  (phi: 256, pol_idx: 50)
> Coordinates:
>   R        (pol_idx, phi) float64 1.0 1.0 1.0 1.0 1.0 ... -1.0 -1.0 -1.0 -1.0
>   Z        (pol_idx, phi) float64 0.0 0.0 0.0 ... 1.225e-16 1.225e-16
>  * pol_idx  (pol_idx) int64 0 1 2 3 4 5 6 7 8 9 ... 41 42 43 44 45 46 47 48 49
>   * phi      (phi) float64 0.0 0.02454 0.04909 0.07363 ... 6.21 6.234 6.259
> Data variables:
>     *empty*

Comment:

I am not sure whether this is a bug. However, I find this behaviour counterintuitive.

Environment:

No files or tests are named. Start by running the provided xarray stack/unstack example and compare the original coordinates with the unstacked result; done means unstack restores the non-indexed R and Z coordinates without changing their intended dimensions.