Repository: pydata/xarray

Issue:
Can't remove coordinates attribute from DataArrays 
<!-- Please include a self-contained copy-pastable example that generates the issue if possible.

Please be concise with code posted. See guidelines below on how to provide a good bug report:

- Craft Minimal Bug Reports: http://matthewrocklin.com/blog/work/2018/02/28/minimal-bug-reports
- Minimal Complete Verifiable Examples: https://stackoverflow.com/help/mcve

Bug reports that follow these guidelines are easier to diagnose, and so are often handled much more quickly.
-->

**What happened**:
Coordinates added to some variables unexpectedly. 

I noticed this after outputting to netCDF. What I have:
```
variables:
        double time(time) ;
                time:bounds = "time_bnds" ;
                time:axis = "T" ;
                time:long_name = "valid_time" ;
                time:standard_name = "time" ;
                time:units = "days since 1850-01-01" ;
                time:calendar = "gregorian" ;
        double time_bnds(time, bnds) ;
                time_bnds:_FillValue = NaN ;
                time_bnds:coordinates = "reftime leadtime height" ;
        double lat(lat) ;
                lat:bounds = "lat_bnds" ;
                lat:units = "degrees_north" ;
                lat:axis = "Y" ;
                lat:long_name = "latitude" ;
                lat:standard_name = "latitude" ;
        double lat_bnds(lat, bnds) ;
                lat_bnds:_FillValue = NaN ;
                lat_bnds:coordinates = "reftime height" ;
        double lon(lon) ;
                lon:bounds = "lon_bnds" ;
                lon:units = "degrees_east" ;
                lon:axis = "X" ;
                lon:long_name = "Longitude" ;
                lon:standard_name = "longitude" ;
        double lon_bnds(lon, bnds) ;
                lon_bnds:_FillValue = NaN ;
                lon_bnds:coordinates = "reftime height" ;
        double height ;
                height:_FillValue = NaN ;
                height:units = "m" ;
                height:axis = "Z" ;
                height:positive = "up" ;
                height:long_name = "height" ;
                height:standard_name = "height" ;
        float tas(time, lat, lon) ;
                tas:_FillValue = 1.e+20f ;
                tas:standard_name = "air_temperature" ;
                tas:long_name = "Near-Surface Air Temperature" ;
                tas:comment = "near-surface (usually, 2 meter) air temperature" ;
                tas:units = "K" ;
                tas:cell_methods = "area: time: mean" ;
                tas:cell_measures = "area: areacella" ;
                tas:history = "2019-05-11T15:53:32Z altered by CMOR: Treated scalar dimension: \'height\'. 2019-05-11T15:53:32Z altered by CMOR: Reordered dimensions, original order: lat lon time." ;
                tas:coordinates = "height reftime leadtime" ;
                tas:missing_value = 1.e+20f ;
        int reftime ;
                reftime:long_name = "Start date of the forecast" ;
                reftime:standard_name = "forecast_reference_time" ;
                reftime:units = "days since 1850-01-01" ;
                reftime:calendar = "gregorian" ;
        double leadtime(time) ;
                leadtime:_FillValue = NaN ;
                leadtime:long_name = "Time elapsed since the start of the forecast" ;
                leadtime:standard_name = "forecast_period" ;
                leadtime:units = "days" ;
        int realization ;
                realization:long_name = "realization" ;
                realization:comment = "For more information on the ripf, refer to the variant_label, initialization_description, physics_description and forcing_description global attributes" ;
                realization:coordinates = "reftime height" ;
```

On `time_bnds`, `lon_bnds`, `lat_bnds` and `realization` there is `coordinates` that I wouldn't expect to be there.


**What you expected to happen**:
Looking only at the coordinates attribute, I expected my ncdump to show:

```
variables:
        int reftime ;
                reftime:long_name = "Start date of the forecast" ;
                reftime:standard_name = "forecast_reference_time" ;
                reftime:units = "days since 1850-01-01" ;
                reftime:calendar = "gregorian" ;
        double leadtime(time) ;
                leadtime:long_name = "Time elapsed since the start of the forecast" ;
                leadtime:standard_name = "forecast_period" ;
                leadtime:units = "days" ;
        int realization ;
                realization:long_name = "realization" ;
                realization:comment = "For more information on the ripf, refer to the variant_label, initialization_description, physics_description and forcing_description global attributes" ;
        double time(time) ;
                time:bounds = "time_bnds" ;
                time:axis = "T" ;
                time:standard_name = "time" ;
                time:units = "days since 1850-01-01" ;
                time:calendar = "gregorian" ;
                time:long_name = "valid_time" ;
        double time_bnds(time, bnds) ;
                time_bnds:units = "days since 1850-01-01" ;
        double lat(lat) ;
                lat:bounds = "lat_bnds" ;
                lat:units = "degrees_north" ;
                lat:axis = "Y" ;
                lat:long_name = "latitude" ;
                lat:standard_name = "latitude" ;
        double lat_bnds(lat, bnds) ;
        double lon(lon) ;
                lon:bounds = "lon_bnds" ;
                lon:units = "degrees_east" ;
                lon:axis = "X" ;
                lon:long_name = "Longitude" ;
                lon:standard_name = "longitude" ;
        double lon_bnds(lon, bnds) ;
        double height ;
                height:units = "m" ;
                height:axis = "Z" ;
                height:positive = "up" ;
                height:long_name = "height" ;
                height:standard_name = "height" ;
        float tas(time, lat, lon) ;
                tas:standard_name = "air_temperature" ;
                tas:long_name = "Near-Surface Air Temperature" ;
                tas:comment = "near-surface (usually, 2 meter) air temperature" ;
                tas:units = "K" ;
                tas:cell_methods = "area: time: mean" ;
                tas:cell_measures = "area: areacella" ;
                tas:history = "2019-05-11T15:53:32Z altered by CMOR: Treated scalar dimension: \'height\'. 2019-05-11T15:53:32Z altered by CMOR: Reordered dimensions, original order: lat lon time." ;
                tas:missing_value = 1.e+20f ;
                tas:_FillValue = 1.e+20f ;
                tas:coordinates = "height reftime leadtime" ;
```

I tried to remove this in the xarray dataset, but whatever I tried they always ended up back in there:

```
>>> import xarray as xr
>>> ds = xr.open_dataset("file.nc", use_cftime=True)

# show coords on realization
>>> ds.realization
<xarray.DataArray 'realization' ()>
array(1, dtype=int32)
Coordinates:
    height   float64 ...
    reftime  object ...
Attributes:
    long_name:  realization
    comment:    For more information on the ripf, refer to the variant_label,...

# try reset_coords - removes the coords
>>> ds.realization.reset_coords(names=["height", "reftime"], drop=True)
<xarray.DataArray 'realization' ()>
array(1, dtype=int32)
Attributes:
    long_name:  realization
    comment:    For more information on the ripf, refer to the variant_label,...

# set realization with result of reset_coords
>>> ds["realization"] = ds.realization.reset_coords(names=["height", "reftime"], drop=True)

# coords back in
>>> ds.realization
<xarray.DataArray 'realization' ()>
array(1, dtype=int32)
Coordinates:
    height   float64 ...
    reftime  object ...
Attributes:
    long_name:  realization
    comment:    For more information on the ripf, refer to the variant_label,...

# try drop_vars - same thing happens
>>> ds.realization.drop_vars(("height", "reftime"))
<xarray.DataArray 'realization' ()>
array(1, dtype=int32)
Attributes:
    long_name:  realization
    comment:    For more information on the ripf, refer to the variant_label,...
    
>>> ds["realization"] = ds.realization.drop_vars(("height", "reftime"))

>>> ds.realization
<xarray.DataArray 'realization' ()>
array(1, dtype=int32)
Coordinates:
    height   float64 ...
    reftime  object ...
Attributes:
    long_name:  realization
    comment:    For more information on the ripf, refer to the variant_label,...
    
# tried creating a new variable to see if the same thing happens  - it does
>>> ds["test"] = ds.realization.drop_vars(("height", "reftime")) 
>>> ds.test
<xarray.DataArray 'test' ()>
array(1, dtype=int32)
Coordinates:
    height   float64 ...
    reftime  object ...
Attributes:
    long_name:  realization
    comment:    For more information on the ripf, refer to the variant_label,...
```

This seems like incorrect behaviour, but perhaps it is expected?

**Minimal Complete Verifiable Example**:

```python
>>> data = xr.DataArray(np.random.randn(2, 3), dims=("x", "y"), coords={"x": [10, 20]})
>>> ds = xr.Dataset({"foo": data, "bar": ("x", [1, 2]), "fake": 10})
>>> ds = ds.assign_coords({"reftime":np.array("2004-11-01T00:00:00", dtype=np.datetime64)}) 
>>> ds = ds.assign({"test": 1})

>>> ds.test
<xarray.DataArray 'test' ()>
array(1)
Coordinates:
    reftime  datetime64[ns] 2004-11-01
    
>>> ds.test.reset_coords(names=["reftime"], drop=True)
<xarray.DataArray 'test' ()>
array(1)

>>> ds["test"] = ds.test.reset_coords(names=["reftime"], drop=True)
>>> ds.test
<xarray.DataArray 'test' ()>
array(1)
Coordinates:
    reftime  datetime64[ns] 2004-11-01
    
 ds.to_netcdf("file.nc")
```

```
ncdump -h file.nc
netcdf file {
dimensions:
	x = 2 ;
	y = 3 ;
variables:
	int64 x(x) ;
	double foo(x, y) ;
		foo:_FillValue = NaN ;
		foo:coordinates = "reftime" ;
	int64 bar(x) ;
		bar:coordinates = "reftime" ;
	int64 fake ;
		fake:coordinates = "reftime" ;
	int64 reftime ;
		reftime:units = "days since 2004-11-01 00:00:00" ;
		reftime:calendar = "proleptic_gregorian" ;
	int64 test ;
		test:coordinates = "reftime" ;
}
```

**Environment**:

<details><summary>Output of <tt>xr.show_versions()</tt></summary>

INSTALLED VERSIONS
------------------
commit: None
python: 3.7.3 (default, Mar 27 2019, 16:54:48) 
[Clang 4.0.1 (tags/RELEASE_401/final)]
python-bits: 64
OS: Darwin
OS-release: 18.7.0
machine: x86_64
processor: i386
byteorder: little
LC_ALL: None
LANG: en_GB.UTF-8
LOCALE: ('en_GB', 'UTF-8')
libhdf5: 1.10.5
libnetcdf: 4.6.3

xarray: 0.18.2
pandas: 1.1.3
numpy: 1.19.2
scipy: None
netCDF4: 1.5.4
pydap: None
h5netcdf: None
h5py: None
Nio: None
zarr: None
cftime: 1.4.1
nc_time_axis: None
PseudoNetCDF: None
rasterio: None
cfgrib: None
iris: None
bottleneck: 1.3.2
dask: 2.30.0
distributed: 2.30.0
matplotlib: None
cartopy: None
seaborn: None
numbagg: None
pint: None
setuptools: 54.1.1
pip: 21.0.1
conda: None
pytest: 6.2.2
IPython: 7.21.0
sphinx: 1.8.1


</details>


Discussion:
github-actions[bot]: ## Unit Test Results
         6 files           6 suites   48m 46s :stopwatch:
16 157 tests 14 441 :heavy_check_mark: 1 716 :zzz: 0 :x:
90 138 runs  82 014 :heavy_check_mark: 8 124 :zzz: 0 :x:

Results for commit f56b37a6.


Patch:
diff --git a/doc/whats-new.rst b/doc/whats-new.rst
index 57ec8539..daddc91b 100644
--- a/doc/whats-new.rst
+++ b/doc/whats-new.rst
@@ -48,6 +48,9 @@ New Features
   By `Justus Magin <https://github.com/keewis>`_.
 - Allow plotting categorical data (:pull:`5464`).
   By `Jimmy Westling <https://github.com/illviljan>`_.
+- Allow removal of the coordinate attribute ``coordinates`` on variables by setting ``.attrs['coordinates']= None``
+  (:issue:`5510`).
+  By `Elle Smith <https://github.com/ellesmith88>`_.
 
 Breaking changes
 ~~~~~~~~~~~~~~~~
diff --git a/xarray/conventions.py b/xarray/conventions.py
index 901d19bd..c3a05e42 100644
--- a/xarray/conventions.py
+++ b/xarray/conventions.py
@@ -751,6 +751,18 @@ def _encode_coordinates(variables, attributes, non_dim_coord_names):
                 f"'coordinates' found in both attrs and encoding for variable {name!r}."
             )
 
+        # if coordinates set to None, don't write coordinates attribute
+        if (
+            "coordinates" in attrs
+            and attrs.get("coordinates") is None
+            or "coordinates" in encoding
+            and encoding.get("coordinates") is None
+        ):
+            # make sure "coordinates" is removed from attrs/encoding
+            attrs.pop("coordinates", None)
+            encoding.pop("coordinates", None)
+            continue
+
         # this will copy coordinates from encoding to attrs if "coordinates" in attrs
         # after the next line, "coordinates" is never in encoding
         # we get support for attrs["coordinates"] for free.
diff --git a/xarray/core/dataarray.py b/xarray/core/dataarray.py
index 7392f7e1..cfb7230d 100644
--- a/xarray/core/dataarray.py
+++ b/xarray/core/dataarray.py
@@ -1902,7 +1902,7 @@ class DataArray(AbstractArray, DataWithCoords, DataArrayArithmetic):
         indexes: Mapping[Hashable, Union[Hashable, Sequence[Hashable]]] = None,
         append: bool = False,
         **indexes_kwargs: Union[Hashable, Sequence[Hashable]],
-    ) -> Optional["DataArray"]:
+    ) -> "DataArray":
         """Set DataArray (multi-)indexes using one or more existing
         coordinates.
 
@@ -1958,7 +1958,7 @@ class DataArray(AbstractArray, DataWithCoords, DataArrayArithmetic):
         self,
         dims_or_levels: Union[Hashable, Sequence[Hashable]],
         drop: bool = False,
-    ) -> Optional["DataArray"]:
+    ) -> "DataArray":
         """Reset the specified index(es) or multi-index level(s).
 
         Parameters


Test patch:
diff --git a/xarray/tests/__init__.py b/xarray/tests/__init__.py
index 64fdbca2..9029dc1c 100644
--- a/xarray/tests/__init__.py
+++ b/xarray/tests/__init__.py
@@ -229,3 +229,15 @@ def create_test_data(seed=None, add_attrs=True):
     obj.encoding = {"foo": "bar"}
     assert all(obj.data.flags.writeable for obj in obj.variables.values())
     return obj
+
+
+_CFTIME_CALENDARS = [
+    "365_day",
+    "360_day",
+    "julian",
+    "all_leap",
+    "366_day",
+    "gregorian",
+    "proleptic_gregorian",
+    "standard",
+]
diff --git a/xarray/tests/test_cftime_offsets.py b/xarray/tests/test_cftime_offsets.py
index cdac9a2d..6d2d9907 100644
--- a/xarray/tests/test_cftime_offsets.py
+++ b/xarray/tests/test_cftime_offsets.py
@@ -26,22 +26,11 @@ from xarray.coding.cftime_offsets import (
     to_cftime_datetime,
     to_offset,
 )
+from xarray.tests import _CFTIME_CALENDARS
 
 cftime = pytest.importorskip("cftime")
 
 
-_CFTIME_CALENDARS = [
-    "365_day",
-    "360_day",
-    "julian",
-    "all_leap",
-    "366_day",
-    "gregorian",
-    "proleptic_gregorian",
-    "standard",
-]
-
-
 def _id_func(param):
     """Called on each parameter passed to pytest.mark.parametrize"""
     return str(param)
diff --git a/xarray/tests/test_conventions.py b/xarray/tests/test_conventions.py
index cd8e3419..ceea1677 100644
--- a/xarray/tests/test_conventions.py
+++ b/xarray/tests/test_conventions.py
@@ -142,6 +142,40 @@ class TestEncodeCFVariable:
         with pytest.raises(ValueError, match=r"'coordinates' found in both attrs"):
             conventions.encode_dataset_coordinates(orig)
 
+    def test_emit_coordinates_attribute_in_attrs(self):
+        orig = Dataset(
+            {"a": 1, "b": 1},
+            coords={"t": np.array("2004-11-01T00:00:00", dtype=np.datetime64)},
+        )
+
+        orig["a"].attrs["coordinates"] = None
+        enc, _ = conventions.encode_dataset_coordinates(orig)
+
+        # check coordinate attribute emitted for 'a'
+        assert "coordinates" not in enc["a"].attrs
+        assert "coordinates" not in enc["a"].encoding
+
+        # check coordinate attribute not emitted for 'b'
+        assert enc["b"].attrs.get("coordinates") == "t"
+        assert "coordinates" not in enc["b"].encoding
+
+    def test_emit_coordinates_attribute_in_encoding(self):
+        orig = Dataset(
+            {"a": 1, "b": 1},
+            coords={"t": np.array("2004-11-01T00:00:00", dtype=np.datetime64)},
+        )
+
+        orig["a"].encoding["coordinates"] = None
+        enc, _ = conventions.encode_dataset_coordinates(orig)
+
+        # check coordinate attribute emitted for 'a'
+        assert "coordinates" not in enc["a"].attrs
+        assert "coordinates" not in enc["a"].encoding
+
+        # check coordinate attribute not emitted for 'b'
+        assert enc["b"].attrs.get("coordinates") == "t"
+        assert "coordinates" not in enc["b"].encoding
+
     @requires_dask
     def test_string_object_warning(self):
         original = Variable(("x",), np.array(["foo", "bar"], dtype=object)).chunk()
diff --git a/xarray/tests/test_missing.py b/xarray/tests/test_missing.py
index e4c74b40..e2dfac04 100644
--- a/xarray/tests/test_missing.py
+++ b/xarray/tests/test_missing.py
@@ -14,6 +14,7 @@ from xarray.core.missing import (
 )
 from xarray.core.pycompat import dask_array_type
 from xarray.tests import (
+    _CFTIME_CALENDARS,
     assert_allclose,
     assert_array_equal,
     assert_equal,
@@ -23,7 +24,6 @@ from xarray.tests import (
     requires_dask,
     requires_scipy,
 )
-from xarray.tests.test_cftime_offsets import _CFTIME_CALENDARS
 
 
 @pytest.fixture
@@ -542,6 +542,7 @@ def test_get_clean_interp_index_dt(cf_da, calendar, freq):
     np.testing.assert_array_equal(gi, si)
 
 
+@requires_cftime
 def test_get_clean_interp_index_potential_overflow():
     da = xr.DataArray(
         [0, 1, 2],
@@ -592,7 +593,10 @@ def test_interpolate_na_max_gap_errors(da_time):
 
 
 @requires_bottleneck
-@pytest.mark.parametrize("time_range_func", [pd.date_range, xr.cftime_range])
+@pytest.mark.parametrize(
+    "time_range_func",
+    [pd.date_range, pytest.param(xr.cftime_range, marks=requires_cftime)],
+)
 @pytest.mark.parametrize("transform", [lambda x: x, lambda x: x.to_dataset(name="a")])
 @pytest.mark.parametrize(
     "max_gap", ["3H", np.timedelta64(3, "h"), pd.to_timedelta("3H")]