Bugfix for issue #191 part 2 (overlooking attributes in xarray.DataArray.encoding) (#264)

* Bugfix for issue #191 part 2

- CFAccessor no longer overlooks cell_measures, formula_terms and bounds when associated attributes are stored in xarray.DataArray.encoding rather than xarray.DataArray.attrs
- Added test and test data

* Replaced netCDF dataset with dummy dataset.

* Execute drop_bounds only for xarray.Datasets

* Added test incl. dataset for drop_bounds

* Removed unnecessary attributes from the new datasets 'ambig' and 'vert'.

* use chainmap rather than nested get

Co-authored-by: Mattia Almansi <[email protected]>

Author: sol1105

cf_xarray/accessor.py

@@ -309,8 +309,9 @@ def _get_measure(obj: Union[DataArray, Dataset], key: str) -> List[str]:
     results = set()
     for var in obj.variables:
         da = obj[var]
-        if "cell_measures" in da.attrs:
-            attr = da.attrs["cell_measures"]
+        attrs_or_encoding = ChainMap(da.attrs, da.encoding)
+        if "cell_measures" in attrs_or_encoding:
+            attr = attrs_or_encoding["cell_measures"]
             measures = parse_cell_methods_attr(attr)
             if key in measures:
                 results.update([measures[key]])
@@ -339,8 +340,9 @@ def _get_bounds(obj: Union[DataArray, Dataset], key: str) -> List[str]:
 
     results = set()
     for var in apply_mapper(_get_all, obj, key, error=False, default=[key]):
-        if "bounds" in obj[var].attrs:
-            results |= {obj[var].attrs["bounds"]}
+        attrs_or_encoding = ChainMap(obj[var].attrs, obj[var].encoding)
+        if "bounds" in attrs_or_encoding:
+            results |= {attrs_or_encoding["bounds"]}
 
     return list(results)
 
@@ -627,8 +629,10 @@ def drop_bounds(names):
         # actual variable. It seems practical to ignore them when indexing
         # with a scalar key. Hopefully these will soon get decoded to IntervalIndex
         # and we can move on...
-        if scalar_key:
-            bounds = {obj[k].attrs.get("bounds", None) for k in names}
+        if not isinstance(obj, DataArray) and scalar_key:
+            bounds = set()
+            for name in names:
+                bounds.update(obj.cf.bounds.get(name, []))
             names = set(names) - bounds
         return names
 
@@ -1364,12 +1368,16 @@ def cell_measures(self) -> Dict[str, List[str]]:
         """
 
         obj = self._obj
-        all_attrs = [da.attrs.get("cell_measures", "") for da in obj.coords.values()]
+        all_attrs = [
+            ChainMap(da.attrs, da.encoding).get("cell_measures", "")
+            for da in obj.coords.values()
+        ]
         if isinstance(obj, DataArray):
-            all_attrs += [obj.attrs.get("cell_measures", "")]
+            all_attrs += [ChainMap(obj.attrs, obj.encoding).get("cell_measures", "")]
         elif isinstance(obj, Dataset):
             all_attrs += [
-                da.attrs.get("cell_measures", "") for da in obj.data_vars.values()
+                ChainMap(da.attrs, da.encoding).get("cell_measures", "")
+                for da in obj.data_vars.values()
             ]
 
         keys = {}
@@ -2144,12 +2152,13 @@ def formula_terms(self) -> Dict[str, str]:
             {parametric_coord_name: {standard_term_name: variable_name}}
         """
         da = self._obj
-        if "formula_terms" not in da.attrs:
+        if "formula_terms" not in ChainMap(da.attrs, da.encoding):
             var = da[_single(_get_dims)(da, "Z")[0]]
         else:
             var = da
+
         terms = {}
-        formula_terms = var.attrs.get("formula_terms", "")
+        formula_terms = ChainMap(var.attrs, var.encoding).get("formula_terms", "")
         for mapping in re.sub(r"\s*:\s*", ":", formula_terms).split():
             key, value = mapping.split(":")
             terms[key] = value

cf_xarray/datasets.py

@@ -291,3 +291,172 @@
     },
     name="basin",
 )
+
+
+ambig = xr.Dataset(
+    data_vars={},
+    coords={
+        "lat": ("lat", np.zeros(5)),
+        "lon": ("lon", np.zeros(5)),
+        "vertices_latitude": (["lat", "bnds"], np.zeros((5, 2))),
+        "vertices_longitude": (["lon", "bnds"], np.zeros((5, 2))),
+    },
+)
+ambig["lat"].attrs = {
+    "bounds": "vertices_latitude",
+    "units": "degrees_north",
+    "standard_name": "latitude",
+    "axis": "Y",
+}
+ambig["lon"].attrs = {
+    "bounds": "vertices_longitude",
+    "units": "degrees_east",
+    "standard_name": "longitude",
+    "axis": "X",
+}
+ambig["vertices_latitude"].attrs = {
+    "units": "degrees_north",
+}
+ambig["vertices_longitude"].attrs = {
+    "units": "degrees_east",
+}
+
+
+vert = xr.Dataset.from_dict(
+    {
+        "coords": {
+            "lat": {
+                "dims": ("lat",),
+                "attrs": {
+                    "standard_name": "latitude",
+                    "axis": "Y",
+                    "bounds": "lat_bnds",
+                    "units": "degrees_north",
+                },
+                "data": [0.0, 1.0],
+            },
+            "lon": {
+                "dims": ("lon",),
+                "attrs": {
+                    "standard_name": "longitude",
+                    "axis": "X",
+                    "bounds": "lon_bnds",
+                    "units": "degrees_east",
+                },
+                "data": [0.0, 1.0],
+            },
+            "lev": {
+                "dims": ("lev",),
+                "attrs": {
+                    "standard_name": "atmosphere_hybrid_sigma_pressure_coordinate",
+                    "formula": "p = ap + b*ps",
+                    "formula_terms": "ap: ap b: b ps: ps",
+                    "postitive": "down",
+                    "axis": "Z",
+                    "bounds": "lev_bnds",
+                },
+                "data": [0.0, 1.0],
+            },
+            "time": {
+                "dims": ("time",),
+                "attrs": {
+                    "standard_name": "time",
+                    "axis:": "T",
+                    "bounds": "time_bnds",
+                    "units": "days since 1850-01-01",
+                    "calendar": "proleptic_gregorian",
+                },
+                "data": [0.5],
+            },
+            "lat_bnds": {
+                "dims": (
+                    "lat",
+                    "bnds",
+                ),
+                "attrs": {
+                    "units": "degrees_north",
+                },
+                "data": [[0.0, 0.5], [0.5, 1.0]],
+            },
+            "lon_bnds": {
+                "dims": (
+                    "lon",
+                    "bnds",
+                ),
+                "attrs": {
+                    "units": "degrees_east",
+                },
+                "data": [[0.0, 0.5], [0.5, 1.0]],
+            },
+            "lev_bnds": {
+                "dims": (
+                    "lev",
+                    "bnds",
+                ),
+                "attrs": {
+                    "standard_name": "atmosphere_hybrid_sigma_pressure_coordinate",
+                    "formula": "p = ap + b*ps",
+                    "formula_terms": "ap: ap b: b ps: ps",
+                },
+                "data": [[0.0, 0.5], [0.5, 1.0]],
+            },
+            "time_bnds": {
+                "dims": ("time", "bnds"),
+                "attrs": {
+                    "units": "days since 1850-01-01",
+                    "calendar": "proleptic_gregorian",
+                },
+                "data": [[0.0, 1.0]],
+            },
+            "ap": {
+                "dims": ("lev",),
+                "data": [0.0, 0.0],
+            },
+            "b": {
+                "dims": ("lev",),
+                "data": [1.0, 0.9],
+            },
+            "ap_bnds": {
+                "dims": (
+                    "lev",
+                    "bnds",
+                ),
+                "data": [[0.0, 0.0], [0.0, 0.0]],
+            },
+            "b_bnds": {
+                "dims": (
+                    "lev",
+                    "bnds",
+                ),
+                "data": [[1.0, 0.95], [0.95, 0.9]],
+            },
+        },
+        "dims": {"time": 1, "lev": 2, "lat": 2, "lon": 2, "bnds": 2},
+        "data_vars": {
+            "o3": {
+                "dims": ("time", "lev", "lat", "lon"),
+                "attrs": {
+                    "cell_methods": "area: time: mean",
+                    "cell_measures": "area: areacella",
+                    "missing_value": 1e20,
+                    "_FillValue": 1e20,
+                },
+                "data": np.ones(8, dtype=np.float32).reshape((1, 2, 2, 2)),
+            },
+            "areacella": {
+                "dims": ("lat", "lon"),
+                "attrs": {
+                    "standard_name": "cell_area",
+                    "cell_methods": "area: sum",
+                    "missing_value": 1e20,
+                    "_FillValue": 1e20,
+                },
+                "data": np.ones(4, dtype=np.float32).reshape((2, 2)),
+            },
+            "ps": {
+                "dims": ("time", "lat", "lon"),
+                "data": np.ones(4, dtype=np.float32).reshape((1, 2, 2)),
+            },
+        },
+    }
+)

cf_xarray/tests/test_accessor.py

@@ -16,6 +16,7 @@
 
 from ..datasets import (
     airds,
+    ambig,
     anc,
     basin,
     ds_no_attrs,
@@ -24,6 +25,7 @@
     multiple,
     popds,
     romsds,
+    vert,
 )
 from . import raise_if_dask_computes, requires_pint
 
@@ -211,6 +213,30 @@ def test_standard_names():
     assert dsnew.cf.standard_names == dict(a=["a", "b"])
 
 
+def test_drop_bounds():
+    assert ambig.cf["latitude"].name == "lat"
+    assert ambig.cf["longitude"].name == "lon"
+    assert ambig.cf.bounds["latitude"] == ["vertices_latitude"]
+    assert ambig.cf.bounds["longitude"] == ["vertices_longitude"]
+
+
+def test_accessor_getattr_and_describe():
+    ds_verta = vert.set_coords(
+        (
+            "ps",
+            "areacella",
+        )
+    )
+    ds_vertb = xr.decode_cf(vert, decode_coords="all")
+
+    assert ds_verta.cf.cell_measures == ds_vertb.cf.cell_measures
+    assert ds_verta.o3.cf.cell_measures == ds_vertb.o3.cf.cell_measures
+    assert ds_verta.cf.formula_terms == ds_vertb.cf.formula_terms
+    assert ds_verta.o3.cf.formula_terms == ds_vertb.o3.cf.formula_terms
+    assert ds_verta.cf.bounds == ds_vertb.cf.bounds
+    assert str(ds_verta.cf) == str(ds_vertb.cf)
+
+
 def test_getitem_standard_name():
     actual = airds.cf["air_temperature"]
     expected = airds["air"]