add cf.grid_mapping_names (#391)

* added grid_mapping attribute

* added get_grid_mapping_name function

* updated rotds with grid_mapping

* renamed dummy variable

* updated grid mapping scalar

* update repr tests and rotds dims

* api docs update

* updated examples in grid_mappings docstring

* added references

* remove zero dim from rotds

* format fix

* added basic grid mappings section in docs

* update section

* update docs section

* added grid_mappgin to associated variables

* added grid_mapping accessor tests

* updated docs section

* syntax fix

* eval-rst

* typo

* fix eval-rst

* polish

* Update cf_xarray/accessor.py

Co-authored-by: Deepak Cherian <[email protected]>

* Update cf_xarray/accessor.py

Co-authored-by: Deepak Cherian <[email protected]>

* Update cf_xarray/accessor.py

Co-authored-by: Deepak Cherian <[email protected]>

* Update cf_xarray/accessor.py

Co-authored-by: Deepak Cherian <[email protected]>

* Update doc/grid_mappings.md

Co-authored-by: Deepak Cherian <[email protected]>

* updated see also sections

* added warning tests

* extract warning tests

* remove skip_grid_mappings

* added dataarray grid_mapping accessor

* added dataarray grid_mapping_name

* docstrings

* docs updates

* fix

* skip grid_mapping coords repr

* cf.grid_mappings switchted logic

* update grid_mapping accessor logic

* test updates

* revert changes

* fixed grid_mapping key critera

* adapted rotds dim names

* small cleanup

* Add "grid_mapping" as special name

* Swtich to .cf.grid_mapping_names

* Avoid DataArray.cf.grid_mapping

* Avoid misleading n/a grid mapping in DataArray repr

* Fix repr

* Fix typing

* Update docs.

* fix cartopy cell

* fix typing

* Add to quickstart too

* docs udpates

---------

Co-authored-by: Deepak Cherian <[email protected]>
Co-authored-by: dcherian <[email protected]>

Author: larsbuntemeyer

cf_xarray/accessor.py

@@ -30,7 +30,12 @@
 from xarray.core.rolling import Coarsen, Rolling
 from xarray.core.weighted import Weighted
 
-from .criteria import cf_role_criteria, coordinate_criteria, regex
+from .criteria import (
+    cf_role_criteria,
+    coordinate_criteria,
+    grid_mapping_var_criteria,
+    regex,
+)
 from .helpers import _guess_bounds_1d, _guess_bounds_2d, bounds_to_vertices
 from .options import OPTIONS
 from .utils import (
@@ -369,6 +374,41 @@ def _get_bounds(obj: DataArray | Dataset, key: Hashable) -> list[Hashable]:
     return list(results)
 
 
+def _get_grid_mapping_name(obj: DataArray | Dataset, key: str) -> list[str]:
+    """
+    Translate from grid mapping name attribute to appropriate variable name.
+    This function interprets the ``grid_mapping`` attribute on DataArrays.
+
+    Parameters
+    ----------
+    obj : DataArray, Dataset
+        DataArray belonging to the coordinate to be checked
+    key : str
+        key to check for.
+
+    Returns
+    -------
+    List[str], Variable name(s) in parent xarray object that matches grid_mapping_name `key`
+    """
+
+    if isinstance(obj, DataArray):
+        obj = obj._to_temp_dataset()
+
+    results = set()
+    for var in obj.variables:
+        da = obj[var]
+        attrs_or_encoding = ChainMap(da.attrs, da.encoding)
+        if "grid_mapping" in attrs_or_encoding:
+            grid_mapping_var_name = attrs_or_encoding["grid_mapping"]
+            if grid_mapping_var_name not in obj.variables:
+                raise ValueError(
+                    f"{var} defines non-existing grid_mapping variable {grid_mapping_var_name}."
+                )
+            if key == obj[grid_mapping_var_name].attrs["grid_mapping_name"]:
+                results.update([grid_mapping_var_name])
+    return list(results)
+
+
 def _get_with_standard_name(
     obj: DataArray | Dataset, name: Hashable | Iterable[Hashable]
 ) -> list[Hashable]:
@@ -395,8 +435,10 @@ def _get_all(obj: DataArray | Dataset, key: Hashable) -> list[Hashable]:
     all_mappers: tuple[Mapper] = (
         _get_custom_criteria,
         functools.partial(_get_custom_criteria, criteria=cf_role_criteria),  # type: ignore
+        functools.partial(_get_custom_criteria, criteria=grid_mapping_var_criteria),
         _get_axis_coord,
         _get_measure,
+        _get_grid_mapping_name,
         _get_with_standard_name,
     )
     results = apply_mapper(all_mappers, obj, key, error=False, default=None)
@@ -706,6 +748,15 @@ def check_results(names, key):
         measures = []
         warnings.warn("Ignoring bad cell_measures attribute.", UserWarning)
 
+    if isinstance(obj, Dataset):
+        grid_mapping_names = list(accessor.grid_mapping_names)
+    else:
+        try:
+            grid_mapping_names = [accessor.grid_mapping_name]
+        except ValueError:
+            grid_mapping_names = []
+    grid_mapping_names.append("grid_mapping")
+
     custom_criteria = ChainMap(*OPTIONS["custom_criteria"])
 
     varnames: list[Hashable] = []
@@ -724,6 +775,12 @@ def check_results(names, key):
             successful[k] = bool(measure)
             if measure:
                 varnames.extend(measure)
+        elif "grid_mapping_names" not in skip and k in grid_mapping_names:
+            grid_mapping = _get_all(obj, k)
+            check_results(grid_mapping, k)
+            successful[k] = bool(grid_mapping)
+            if grid_mapping:
+                varnames.extend(grid_mapping)
         elif k in custom_criteria or k in cf_role_criteria:
             names = _get_all(obj, k)
             check_results(names, k)
@@ -1415,13 +1472,15 @@ def make_text_section(subtitle, attr, valid_values=None, default_keys=None):
         text += make_text_section("Standard Names", "standard_names", coords)
         text += make_text_section("Bounds", "bounds", coords)
         if isinstance(self._obj, Dataset):
+            text += make_text_section("Grid Mappings", "grid_mapping_names", coords)
             data_vars = self._obj.data_vars
             text += "\nData Variables:"
             text += make_text_section(
                 "Cell Measures", "cell_measures", data_vars, _CELL_MEASURES
             )
             text += make_text_section("Standard Names", "standard_names", data_vars)
             text += make_text_section("Bounds", "bounds", data_vars)
+            text += make_text_section("Grid Mappings", "grid_mapping_names", data_vars)
 
         return text
 
@@ -1442,6 +1501,14 @@ def keys(self) -> set[Hashable]:
         varnames.extend(list(self.cell_measures))
         varnames.extend(list(self.standard_names))
         varnames.extend(list(self.cf_roles))
+        if isinstance(self._obj, xr.Dataset):
+            varnames.extend(list(self.grid_mapping_names))
+        else:
+            try:
+                gmname = self.grid_mapping_name
+                varnames.extend(list(gmname))
+            except ValueError:
+                pass
 
         return set(varnames)
 
@@ -1604,6 +1671,7 @@ def get_associated_variable_names(
             2. "bounds"
             3. "cell_measures"
             4. "coordinates"
+            5. "grid_mapping"
         to a list of variable names referred to in the appropriate attribute
 
         Parameters
@@ -1618,7 +1686,13 @@ def get_associated_variable_names(
         names : dict
             Dictionary with keys "ancillary_variables", "cell_measures", "coordinates", "bounds".
         """
-        keys = ["ancillary_variables", "cell_measures", "coordinates", "bounds"]
+        keys = [
+            "ancillary_variables",
+            "cell_measures",
+            "coordinates",
+            "bounds",
+            "grid_mapping",
+        ]
         coords: dict[str, list[Hashable]] = {k: [] for k in keys}
         attrs_or_encoding = ChainMap(self._obj[name].attrs, self._obj[name].encoding)
 
@@ -1660,6 +1734,9 @@ def get_associated_variable_names(
                 if dbounds:
                     coords["bounds"].append(dbounds)
 
+        if "grid_mapping" in attrs_or_encoding:
+            coords["grid_mapping"] = [attrs_or_encoding["grid_mapping"]]
+
         allvars = itertools.chain(*coords.values())
         missing = set(allvars) - set(self._maybe_to_dataset()._variables)
         if missing:
@@ -2048,6 +2125,8 @@ def __getitem__(self, key: Hashable | Iterable[Hashable]) -> DataArray | Dataset
               - cell measures: "area", "volume", or other names present in the \
                              ``cell_measures`` attribute
               - standard names: names present in ``standard_name`` attribute
+              - cf roles: 'timeseries_id', 'profile_id', 'trajectory_id', 'mesh_topology', 'grid_topology'
+              - grid mappings: 'grid_mapping' or a grid_mapping_name like 'rotated_latitude_longitude'
 
         Returns
         -------
@@ -2372,6 +2451,51 @@ def bounds_to_vertices(
                 )
         return obj
 
+    @property
+    def grid_mapping_names(self) -> dict[str, list[str]]:
+        """
+        Property that returns a dictionary mapping the CF grid mapping name
+        to the variable name containing the grid mapping attributes.
+
+        Returns
+        -------
+        dict
+            Dictionary mapping the CF grid mapping name to the grid mapping variable name.
+
+        See Also
+        --------
+        DataArray.cf.grid_mapping
+
+        References
+        ----------
+        Please refer to the CF conventions document : https://cfconventions.org/Data/cf-conventions/cf-conventions-1.10/cf-conventions.html#grid-mappings-and-projections
+
+        For a list of valid grid_mapping names, refer to: https://cfconventions.org/Data/cf-conventions/cf-conventions-1.10/cf-conventions.html#appendix-grid-mappings
+
+        Examples
+        --------
+        >>> from cf_xarray.datasets import rotds
+        >>> rotds.cf.grid_mapping_names
+        {'rotated_latitude_longitude': ['rotated_pole']}
+        """
+
+        obj = self._obj
+        keys = set(obj.variables)
+
+        vardict = {
+            key: obj.variables[key].attrs["grid_mapping_name"]
+            for key in keys
+            if "grid_mapping_name" in obj.variables[key].attrs
+        }
+
+        results = {}
+        for k, v in vardict.items():
+            if v not in results:
+                results[v] = [k]
+            else:
+                results[v].append(k)
+        return results
+
     def decode_vertical_coords(self, *, outnames=None, prefix=None):
         """
         Decode parameterized vertical coordinates in place.
@@ -2547,6 +2671,43 @@ def formula_terms(self) -> dict[str, str]:
             terms[key] = value
         return terms
 
+    @property
+    def grid_mapping_name(self) -> str:
+        """
+        Get CF grid mapping name associated with this variable.
+
+        Parameters
+        ----------
+        key : str
+            Name of variable whose grid_mapping name is desired.
+
+        Returns
+        -------
+        str
+            CF Name of the associated grid mapping.
+
+        See Also
+        --------
+        Dataset.cf.grid_mapping_names
+
+        Examples
+        --------
+        >>> from cf_xarray.datasets import rotds
+        >>> rotds.cf["temp"].cf.grid_mapping_name
+        'rotated_latitude_longitude'
+
+        """
+
+        da = self._obj
+
+        attrs_or_encoding = ChainMap(da.attrs, da.encoding)
+        grid_mapping = attrs_or_encoding.get("grid_mapping", None)
+        if not grid_mapping:
+            raise ValueError("No 'grid_mapping' attribute present.")
+
+        grid_mapping_var = da[grid_mapping]
+        return grid_mapping_var.attrs["grid_mapping_name"]
+
     def __getitem__(self, key: Hashable | Iterable[Hashable]) -> DataArray:
         """
         Index into a DataArray making use of CF attributes.
@@ -2561,6 +2722,8 @@ def __getitem__(self, key: Hashable | Iterable[Hashable]) -> DataArray:
                              ``cell_measures`` attribute
               - standard names: names present in ``standard_name`` attribute of \
                 coordinate variables
+              - cf roles: 'timeseries_id', 'profile_id', 'trajectory_id', 'mesh_topology', 'grid_topology'
+              - grid mappings: 'grid_mapping' or a grid_mapping_name like 'rotated_latitude_longitude'
 
         Returns
         -------

cf_xarray/criteria.py

@@ -4,9 +4,12 @@
 Copyright (c) 2017 MetPy Developers.
 """
 
+try:
+    import regex as re
+except ImportError:
+    import re  # type: ignore
 
-import re
-from typing import Mapping, MutableMapping, Tuple
+from typing import Any, Mapping, MutableMapping, Tuple
 
 cf_role_criteria: Mapping[str, Mapping[str, str]] = {
     k: {"cf_role": k}
@@ -22,6 +25,11 @@
     )
 }
 
+# A grid mapping varibale is anything with a grid_mapping_name attribute
+grid_mapping_var_criteria: Mapping[str, Mapping[str, Any]] = {
+    "grid_mapping": {"grid_mapping_name": re.compile(".")}
+}
+
 coordinate_criteria: MutableMapping[str, MutableMapping[str, Tuple]] = {
     "latitude": {
         "standard_name": ("latitude",),

cf_xarray/datasets.py

@@ -228,6 +228,10 @@ def _create_mollw_ds():
 def _create_inexact_bounds():
     # Dataset that creates rotated pole curvilinear coordinates with CF bounds in
     # counterclockwise order that have precision issues.
+    # dataset created using: https://gist.github.com/larsbuntemeyer/105d83c1eb39b1462150d3fabca0b66b
+    rlon = np.array([17.935, 18.045, 18.155])
+    rlat = np.array([21.615, 21.725, 21.835])
+
     lon = np.array(
         [
             [64.21746363939087, 64.42305921561967, 64.62774455060337],
@@ -296,23 +300,70 @@ def _create_inexact_bounds():
 
     rotated = xr.Dataset(
         coords=dict(
+            rlon=xr.DataArray(
+                rlon,
+                dims="rlon",
+                attrs={
+                    "units": "degrees",
+                    "axis": "X",
+                    "standard_name": "grid_longitude",
+                },
+            ),
+            rlat=xr.DataArray(
+                rlat,
+                dims="rlat",
+                attrs={
+                    "units": "degrees",
+                    "axis": "Y",
+                    "standard_name": "grid_latitude",
+                },
+            ),
             lon=xr.DataArray(
                 lon,
-                dims=("x", "y"),
-                attrs={"units": "degrees_east", "bounds": "lon_bounds"},
+                dims=("rlon", "rlat"),
+                attrs={
+                    "units": "degrees_east",
+                    "bounds": "lon_bounds",
+                    "standard_name": "longitude",
+                },
             ),
             lat=xr.DataArray(
                 lat,
-                dims=("x", "y"),
-                attrs={"units": "degrees_north", "bounds": "lat_bounds"},
+                dims=("rlon", "rlat"),
+                attrs={
+                    "units": "degrees_north",
+                    "bounds": "lat_bounds",
+                    "standard_name": "latitude",
+                },
             ),
         ),
         data_vars=dict(
             lon_bounds=xr.DataArray(
-                lon_bounds, dims=("bounds", "x", "y"), attrs={"units": "degrees_east"}
+                lon_bounds,
+                dims=("bounds", "rlon", "rlat"),
+                attrs={"units": "degrees_east"},
             ),
             lat_bounds=xr.DataArray(
-                lat_bounds, dims=("bounds", "x", "y"), attrs={"units": "degrees_north"}
+                lat_bounds,
+                dims=("bounds", "rlon", "rlat"),
+                attrs={"units": "degrees_north"},
+            ),
+            rotated_pole=xr.DataArray(
+                np.zeros((), dtype=np.int32),
+                dims=None,
+                attrs={
+                    "grid_mapping_name": "rotated_latitude_longitude",
+                    "grid_north_pole_latitude": 39.25,
+                    "grid_north_pole_longitude": -162.0,
+                },
+            ),
+            temp=xr.DataArray(
+                np.random.rand(3, 3),
+                dims=("rlat", "rlon"),
+                attrs={
+                    "standard_name": "air_temperature",
+                    "grid_mapping": "rotated_pole",
+                },
             ),
         ),
     )