Cannot recover original Dataset structure with to_unstacked_dataset after to_stacked_array

MCVE Code Sample

temp = 15 + 8 * np.random.randn(2, 2, 3)
precip = 10 * np.random.rand(2, 2, 3)
lon = [[-99.83, -99.32], [-99.79, -99.23]]
lat = [[42.25, 42.21], [42.63, 42.59]]

ds = xr.Dataset({'temperature': (['x', 'y', 'time'],  temp),
                'precipitation': (['x', 'y', 'time'], precip)},
                coords={'lon': (['x', 'y'], lon),
                'lat': (['x', 'y'], lat),
                'time': pd.date_range('2014-09-06', periods=3),
                'reference_time': pd.Timestamp('2014-09-05')})

ds.to_stacked_array('features', sample_dims=['time']).to_unstacked_dataset('features')

<xarray.Dataset>
Dimensions:         (features: 4, time: 3)
Coordinates:
    reference_time  datetime64[ns] 2014-09-05
  * time            (time) datetime64[ns] 2014-09-06 2014-09-07 2014-09-08
    lon             (features) float64 -99.83 -99.32 -99.79 -99.23
    lat             (features) float64 42.25 42.21 42.63 42.59
  * features        (features) MultiIndex
  - x               (features) int64 0 0 1 1
  - y               (features) int64 0 1 0 1
Data variables:
    precipitation   (time, features) float64 7.405 9.145 6.56 ... 5.423 2.48
    temperature     (time, features) float64 18.55 17.48 35.49 ... 13.13 20.8

Expected Output

The original array:

<xarray.Dataset>
Dimensions:         (time: 3, x: 2, y: 2)
Coordinates:
    lon             (x, y) float64 -99.83 -99.32 -99.79 -99.23
    lat             (x, y) float64 42.25 42.21 42.63 42.59
  * time            (time) datetime64[ns] 2014-09-06 2014-09-07 2014-09-08
    reference_time  datetime64[ns] 2014-09-05
Dimensions without coordinates: x, y
Data variables:
    temperature     (x, y, time) float64 13.11 2.632 14.3 ... 15.46 24.5 5.13
    precipitation   (x, y, time) float64 0.9878 4.014 1.916 ... 2.716 0.6272

Problem Description

After stacking a Dataset to array, the original Dataset is not recovered with to_unstacked_dataset -- this is not always the case E.g.

ds = xr.Dataset({'X1': ('Y', [1,2,3]), 'X2': ('Y', [4,5,6])})
ds = ds.to_stacked_array('features', sample_dims=[])
ds.to_unstacked_dataset('features')

can recover the original Dataset.

Output of xr.show_versions()

INSTALLED VERSIONS
------------------
commit: None
python: 3.7.3 (default, Apr 24 2019, 15:29:51) [MSC v.1915 64 bit (AMD64)]
python-bits: 64
OS: Windows
OS-release: 10
machine: AMD64
processor: Intel64 Family 6 Model 69 Stepping 1, GenuineIntel
byteorder: little
LC_ALL: None
LANG: None
LOCALE: None.None
libhdf5: 1.10.4
libnetcdf: 4.6.1

xarray: 0.14.0
pandas: 0.24.2
numpy: 1.16.4
scipy: 1.2.2
netCDF4: 1.4.2
pydap: None
h5netcdf: None
h5py: None
Nio: None
zarr: None
cftime: 1.0.3.4
nc_time_axis: None
PseudoNetCDF: None
rasterio: None
cfgrib: None
iris: None
bottleneck: None
dask: 2.6.0
distributed: 2.6.0
matplotlib: 3.1.1
cartopy: None
seaborn: 0.9.0
numbagg: None
setuptools: 41.0.1
pip: 19.1.1
conda: 4.7.12
pytest: 5.2.1
IPython: 7.9.0
sphinx: None

Start by reading the implementations of Dataset.to_stacked_array and DataArray.to_unstacked_dataset, then reproduce the reported example with the stated xarray versions. Add a regression test covering multidimensional sample dimensions, and consider the work complete when stacking and unstacking recover the original dimensions, coordinates, and data variables.