issues

What happened?

When a netcdf file containing a mapping variable and a height coordinates is opened, the height coordinates is added to the mapping variable.

The minimal example generates a netcdf similar to CMIP6 data remapped on a curvilinear grid (for this example in fact it is also a lat / lon grid), where the curvilinear coordinates x / y are connected to lat / lon and tas with the PlateCarree variable. The PlateCarree variable doesn't contains any coordinates when the netcdf is created, and it is also the case when the netcdf is read with netCDF4. But when it is read with xr.open_dataset, the height coordinates is added to the PlateCarree variable.

What did you expect to happen?

Variables without coordinates should not have any after reading with xarray

Minimal Complete Verifiable Example

```Python import sys,os

import numpy as np import pandas as pd import xarray as xr import netCDF4

Functions

def build_data(): ## Start by build the coordinates x = np.arange( -5 , 10 + 0.1 , 1 ) y = np.arange( 41 , 52 + 0.1 , 1 ) lon,lat = np.meshgrid(x,y) time = pd.date_range( "2000-01-01" , "2005-12-31" ).values

lon = xr.DataArray( lon , dims = ["y","x"] , coords = [y,x] )
lat = xr.DataArray( lat , dims = ["y","x"] , coords = [y,x] )


## Variables
T = xr.DataArray( np.random.normal( size = (time.size,y.size,x.size) ) , dims = ["time","y","x"]  , coords = { "time" : time , "y" : y , "x" : x , "lat" : lat , "lon" : lon } )

## Dataset
odata = xr.Dataset( { "tas" : T , "height" : 2. , "PlateCarree" : 1 } )

## Add attributes
odata.lon.attrs["axis"]          = "x"
odata.lon.attrs["long_name"]     = "longitude coordinate"
odata.lon.attrs["standard_name"] = "longitude"
odata.lon.attrs["units"]         = "degrees_east"

odata.lat.attrs["axis"]          = "y"
odata.lat.attrs["long_name"]     = "latitude coordinate"
odata.lat.attrs["standard_name"] = "latitude"
odata.lat.attrs["units"]         = "degrees_north"

odata.x.attrs["units"]         = "degrees_north"
odata.x.attrs["long_name"]     = "x coordinate of projection"
odata.x.attrs["standard_name"] = "projection_x_coordinate"

odata.y.attrs["units"]         = "degrees_east"
odata.y.attrs["long_name"]     = "y coordinate of projection"
odata.y.attrs["standard_name"] = "projection_y_coordinate"

odata.time.attrs["axis"]          = "T"
odata.time.attrs["standard_name"] = "time"
odata.time.attrs["long_name"]     = "Time axis"

odata.PlateCarree.attrs["grid_mapping_name"] = "plate_carree"
odata.PlateCarree.attrs["epsg"]              = "4326"

odata.tas.attrs["grid_mapping"]  = "PlateCarree"
odata.tas.attrs["standard_name"] = "air_temperature"
odata.tas.attrs["long_name"]     = "Daily Mean Near-Surface Air Temperature"
odata.tas.attrs["units"]         = "Kelvin"
odata.tas.attrs["coordinates"]   = "height lat lon"

odata.height.attrs["standard_name"] = "height"
odata.height.attrs["long_name"] = "height"
odata.height.attrs["units"] = "m"
odata.height.attrs["positive"] = "up"
odata.height.attrs["axis"] = "Z"
odata.height.attrs["name"] = "height"

## Encoding
ny = y.size
nx = x.size
encoding = { "time" : { "dtype" : "double"  , "zlib" : True , "complevel" : 5 , "chunksizes" : (1,) , "units" : "days since 1950-01-01" , "calendar" : "standard" } , 
             "y"    : { "dtype" : "double"  , "zlib" : True , "complevel" : 5 , "chunksizes" : (ny,) } , 
             "x"    : { "dtype" : "double"  , "zlib" : True , "complevel" : 5 , "chunksizes" : (nx,) } , 
             "lat"  : { "dtype" : "double"  , "zlib" : True , "complevel" : 5 , "chunksizes" : (ny,nx) } , 
             "lon"  : { "dtype" : "double"  , "zlib" : True , "complevel" : 5 , "chunksizes" : (ny,nx) } , 
             "tas"  : { "dtype" : "float32" , "zlib" : True , "complevel" : 5 , "chunksizes" : (1,ny,nx) } , 
          "height"  : { "dtype" : "double" } , 
     "PlateCarree"  : { "dtype" : "int32" } }

odata.to_netcdf( "test.nc" , encoding = encoding )

main

if name == "main":

## Start by generate data
build_data()

## Open with xarray
data = xr.open_dataset("test.nc").load()

print("With xarray:")
print(data.PlateCarree) ## height is a coordinate of 'PlateCarree'

## And open with netCDF4
print("="*80)
print("With netCDF4")
with netCDF4.Dataset( "test.nc" , mode = "r" ) as ncfile:
    print(ncfile.variables["PlateCarree"].get_dims()) ## height is not a coordinate of 'PlateCarree'

```

MVCE confirmation

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

Relevant log output

No response

Anything else we need to know?

No response

Environment