Scalars coordinates have no memory on their DataArray

What happened?

I am noticing an issue which might be a bug or a feature, but leaves me quite puzzled. I understand that scalar coordinates are transported all together when slicing/selection is operated to make smooth operations across datasets.

However, when I have a dataset with two variables depending on two different coordinates but sharing a third one, if I operate selection on both the first two coordinates and then select one of the two variable available, I obtain a DataArray where both my coordinates are scalar and there is no way - apparently - to know which coordinate is associated to my DataArray. I tried to provide an example below, and the output looks agnostic to my coordinates y and z. In the example below var1 and var2 have the same dimensions x and the same (scalar) coordinates y and z. Both DataArray could lie on both of coordinate and I have no way to know it without knowing its history.

Thanks for any hint!

What did you expect to happen?

To have an attribute, index or anything else into the final DataArray that helps me to guess what is the coordinate on which my variable depends on. In the example below var1 and var2 have the same dimensions x and the same (scalar) coordinates y and z.

Minimal Complete Verifiable Example

import xarray as xr
import numpy as np

arr = xr.Dataset(
    {'var1': (('y', 'x'), np.zeros((1,5))),
    'var2': (('z', 'x'), np.zeros((1,5)))},
     coords={'y': ('y', [10]), 'z': ('z', [20])})
arr.isel(y=0, z=0)['var1']

Steps to reproduce

No response

MVCE confirmation

• Minimal example — the example is as focused as reasonably possible to demonstrate the underlying issue in xarray.
• Complete example — the example is self-contained, including all data and the text of any traceback.
• Verifiable example — the example copy & pastes into an IPython prompt or Binder notebook, returning the result.
• New issue — a search of GitHub Issues suggests this is not a duplicate.
• Recent environment — the issue occurs with the latest version of xarray and its dependencies.

Relevant log output

Anything else we need to know?

I checked other issues, the only one that seems relevant is https://github.com/pydata/xarray/issues/4501 but I am not sure I fully understand the implications

Environment

Start by running the minimal xarray and NumPy example in the issue, then read the related discussion in issue #4501. Define how a selected DataArray should preserve or expose the coordinate association, and add regression coverage demonstrating the expected result.