diff --git a/cf_xarray/accessor.py b/cf_xarray/accessor.py
index cc8617d9..2ddbd90c 100644
--- a/cf_xarray/accessor.py
+++ b/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,6 +1472,7 @@ 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(
@@ -1422,6 +1480,7 @@ def make_text_section(subtitle, attr, valid_values=None, default_keys=None):
             )
             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
         -------
diff --git a/cf_xarray/criteria.py b/cf_xarray/criteria.py
index 302a2619..f2928cd4 100644
--- a/cf_xarray/criteria.py
+++ b/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",),
diff --git a/cf_xarray/datasets.py b/cf_xarray/datasets.py
index aaadfa32..d2d69e28 100644
--- a/cf_xarray/datasets.py
+++ b/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",
+                },
             ),
         ),
     )
diff --git a/cf_xarray/tests/test_accessor.py b/cf_xarray/tests/test_accessor.py
index fb94dbfa..1a04e6ab 100644
--- a/cf_xarray/tests/test_accessor.py
+++ b/cf_xarray/tests/test_accessor.py
@@ -84,12 +84,16 @@ def test_repr() -> None:
 
     - Bounds:   n/a
 
+    - Grid Mappings:   n/a
+
     Data Variables:
     - Cell Measures:   area, volume: n/a
 
     - Standard Names:   air_temperature: ['air']
 
     - Bounds:   n/a
+
+    - Grid Mappings:   n/a
     """
     assert actual == dedent(expected)
 
@@ -137,6 +141,8 @@ def test_repr() -> None:
 
     - Bounds:   n/a
 
+    - Grid Mappings:   n/a
+
     Data Variables:
     - Cell Measures:   area, volume: n/a
 
@@ -144,6 +150,8 @@ def test_repr() -> None:
                         sea_water_x_velocity: ['UVEL']
 
     - Bounds:   n/a
+
+    - Grid Mappings:   n/a
     """
     assert actual == dedent(expected)
 
@@ -173,12 +181,16 @@ def test_repr() -> None:
 
     - Bounds:   n/a
 
+    - Grid Mappings:   n/a
+
     Data Variables:
     - Cell Measures:   area, volume: n/a
 
     - Standard Names:   air_temperature: [<this-array>]
 
     - Bounds:   n/a
+
+    - Grid Mappings:   n/a
     """
     assert actual == dedent(expected)
 
@@ -199,12 +211,16 @@ def test_repr() -> None:
 
     - Bounds:   n/a
 
+    - Grid Mappings:   n/a
+
     Data Variables:
     - Cell Measures:   area, volume: n/a
 
     - Standard Names:   n/a
 
     - Bounds:   n/a
+
+    - Grid Mappings:   n/a
     """
     assert actual == dedent(expected)
 
@@ -270,6 +286,8 @@ def test_cell_measures() -> None:
                         foo_std_name: ['foo']
 
     - Bounds:   n/a
+
+    - Grid Mappings:   n/a
     """
     assert actual.endswith(dedent(expected_repr))
 
@@ -308,6 +326,7 @@ def test_accessor_getattr_and_describe() -> None:
     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 ds_verta.cf.grid_mapping_names == ds_vertb.cf.grid_mapping_names
     assert str(ds_verta.cf) == str(ds_vertb.cf)
 
 
@@ -836,7 +855,7 @@ def test_add_bounds_nd_variable() -> None:
     # 2D rotated ds
     lon_bounds = (
         rotds.drop_vars(["lon_bounds"])
-        .assign(x=rotds["x"], y=rotds["y"])
+        .assign(x=rotds.cf["X"], y=rotds.cf["Y"])
         .cf.add_bounds(["lon"])
         .lon_bounds
     )
@@ -950,6 +969,95 @@ def test_get_bounds_dim_name() -> None:
     assert mollwds.cf.get_bounds_dim_name("lon") == "bounds"
 
 
+def test_grid_mappings():
+    ds = rotds.copy(deep=False)
+
+    actual = ds.cf.grid_mapping_names
+    expected = {"rotated_latitude_longitude": ["rotated_pole"]}
+    assert actual == expected
+
+    expected = ds.rotated_pole
+    actual = ds.cf["rotated_latitude_longitude"]
+    assert_identical(actual, expected)
+
+    expected = ds.rotated_pole
+    actual = ds.cf["grid_mapping"]
+    assert_identical(actual, expected)
+
+    # Propagation
+    actual = ds.cf[["temp"]]
+    assert "rotated_pole" in actual.coords
+
+    actual = ds.cf["temp"]
+    assert "rotated_pole" in actual.coords
+
+    actual = ds.cf["temp"].cf["grid_mapping"]
+    assert_identical(actual, expected)
+
+    actual = ds.cf["temp"].coords["rotated_pole"].drop_vars("rotated_pole")
+    assert_identical(actual, expected)
+
+    actual = ds.cf["temp"].cf["rotated_latitude_longitude"]
+    assert_identical(actual, expected)
+
+    # Test repr
+    expected = """\
+    - Grid Mappings:   rotated_latitude_longitude: ['rotated_pole']
+    """
+    assert dedent(expected) in ds.cf.__repr__()
+    # assert dedent(expected) in ds.cf["temp"].cf.__repr__()
+
+    # grid_mapping_name
+    assert ds.cf["temp"].cf.grid_mapping_name == "rotated_latitude_longitude"
+
+    # what if there are really 2 grid mappins?
+    ds["temp2"] = ds.temp
+    ds["temp2"].attrs["grid_mapping"] = "rotated_pole2"
+    ds["rotated_pole2"] = ds.rotated_pole
+    expected = """\
+    - Grid Mappings:   rotated_latitude_longitude: ['rotated_pole', 'rotated_pole2']
+    """
+    assert dedent(expected) in ds.cf.__repr__()
+
+    with pytest.raises(KeyError):
+        ds.cf["grid_mapping"]
+
+    assert "rotated_latitude_longitude" in ds.cf.keys()
+
+    with pytest.raises(KeyError):
+        ds.cf["grid_mapping"]
+    actual = ds.cf[["grid_mapping"]]
+    expected = ds[["rotated_pole", "rotated_pole2"]].reset_coords()
+    assert_identical(expected, actual)
+
+    expected = ds.rotated_pole
+    actual = ds.cf["temp"].cf["grid_mapping"]
+    assert_identical(actual, expected)
+    expected = ds.rotated_pole2
+    actual = ds.cf["temp2"].cf["grid_mapping"]
+    assert_identical(actual, expected)
+
+    # test _get_all with grid_mapping_var mapper
+    ds = ds.cf.set_coords("grid_mapping")
+    assert "rotated_pole" in ds.coords
+
+
+def test_bad_grid_mapping_attribute():
+    ds = rotds.copy(deep=False)
+    ds.temp.attrs["grid_mapping"] = "foo"
+    # warning when extracting a Datarray and grid_mapping does not exist
+    with pytest.warns(UserWarning):
+        ds.cf["temp"]
+    # warning when extracting a Dataset and grid_mapping does not exist
+    with pytest.warns(UserWarning):
+        ds.cf[["temp"]]
+    # this should probably also raise a warning (like cell_measures)
+    # with pytest.warns(UserWarning):
+    #    assert ds.cf.grid_mappings == {}
+    with pytest.warns(UserWarning):
+        ds.cf.get_associated_variable_names("temp", error=False)
+
+
 def test_docstring() -> None:
     assert "One of ('X'" in airds.cf.groupby.__doc__
     assert "Time variable accessor e.g. 'T.month'" in airds.cf.groupby.__doc__
diff --git a/doc/api.rst b/doc/api.rst
index 3f73c94e..8a9060dc 100644
--- a/doc/api.rst
+++ b/doc/api.rst
@@ -36,6 +36,7 @@ Attributes
     DataArray.cf.cf_roles
     DataArray.cf.coordinates
     DataArray.cf.formula_terms
+    DataArray.cf.grid_mapping_name
     DataArray.cf.is_flag_variable
     DataArray.cf.standard_names
     DataArray.cf.plot
@@ -94,6 +95,7 @@ Attributes
     Dataset.cf.cf_roles
     Dataset.cf.coordinates
     Dataset.cf.formula_terms
+    Dataset.cf.grid_mapping_names
     Dataset.cf.standard_names
 
 .. _dsmeth:
diff --git a/doc/cartopy_rotated_pole.png b/doc/cartopy_rotated_pole.png
new file mode 100644
index 00000000..41472d34
Binary files /dev/null and b/doc/cartopy_rotated_pole.png differ
diff --git a/doc/grid_mappings.md b/doc/grid_mappings.md
new file mode 100644
index 00000000..a653f6cc
--- /dev/null
+++ b/doc/grid_mappings.md
@@ -0,0 +1,115 @@
+---
+jupytext:
+  text_representation:
+    format_name: myst
+kernelspec:
+  display_name: Python 3
+  name: python3
+---
+
+```{eval-rst}
+.. currentmodule:: xarray
+```
+
+# Grid Mappings
+
+See
+
+1. {py:attr}`Dataset.cf.grid_mapping_names`,
+1. {py:attr}`DataArray.cf.grid_mapping_name`
+
+`cf_xarray` understands the concept of coordinate projections using the [grid_mapping](https://cfconventions.org/Data/cf-conventions/cf-conventions-1.10/cf-conventions.html#grid-mappings-and-projections) attribute convention. For example, the dataset might contain two sets of coordinates:
+
+- native coordinates in which the data is defined, e.g., regular 1D coordinates
+- projected coordinates which probably denote some "real" coordinates in [latitude and longitude](https://en.wikipedia.org/wiki/Geographic_coordinate_system#Latitude_and_longitude)
+
+Due to the projection, those real coordinates are probably 2D data variables. The `grid_mapping` attribute of a data variable makes a connection to another data variable defining the coordinate reference system (CRS) of those native coordinates. It should enable you to project the native coordinates into any other CRS, including the real 2D latitude and longitude coordinates. This is often useful for plotting, e.g., you can [tell cartopy how to correctly plot coastlines](https://scitools.org.uk/cartopy/docs/latest/tutorials/understanding_transform.html) for the CRS your data is defined in.
+
+## Extracting grid mapping info
+
+### Dataset
+
+To access `grid_mapping` attributes, consider this example:
+
+```{code-cell}
+from cf_xarray.datasets import rotds
+
+rotds
+```
+
+The related grid mappings can be discovered using `Dataset.cf.grid_mapping_names` which maps a
+["grid mapping name"](http://cfconventions.org/cf-conventions/cf-conventions.html#appendix-grid-mappings) to the
+appropriate variable name:
+
+```{code-cell}
+rotds.cf.grid_mapping_names
+```
+
+Access the `grid_mapping` variable as
+
+```{code-cell}
+rotds.cf["grid_mapping"]
+```
+
+### DataArrays
+
+Grid mapping variables are propagated when extracting DataArrays:
+
+```{code-cell}
+da = rotds.cf["temp"]
+da
+```
+
+To find the grid mapping name use the singular {py:attr}`DataArray.cf.grid_mapping_name`
+
+```{code-cell}
+da.cf.grid_mapping_name
+```
+
+And to get the grid mapping variable
+
+```{code-cell}
+da.cf["grid_mapping"]
+```
+
+## Use `grid_mapping` in projections
+
+The grid mapping information use very useful in projections, e.g., for plotting. [pyproj](https://pyproj4.github.io/pyproj/stable/api/crs/crs.html#pyproj.crs.CRS.from_cf) understands CF conventions right away, e.g.
+
+```python
+from pyproj import CRS
+
+CRS.from_cf(rotds.cf["grid_mapping"].attrs)
+```
+
+gives you more details on the projection:
+
+```
+<Derived Geographic 2D CRS: {"$schema": "https://proj.org/schemas/v0.2/projjso ...>
+Name: undefined
+Axis Info [ellipsoidal]:
+- lon[east]: Longitude (degree)
+- lat[north]: Latitude (degree)
+Area of Use:
+- undefined
+Coordinate Operation:
+- name: Pole rotation (netCDF CF convention)
+- method: Pole rotation (netCDF CF convention)
+Datum: World Geodetic System 1984
+- Ellipsoid: WGS 84
+- Prime Meridian: Greenwich
+```
+
+For use in cartopy, there is some more overhead due to [this issue](https://github.com/SciTools/cartopy/issues/2099). So you should select the right cartopy CRS and just feed in the grid mapping info:
+
+```python
+from cartopy import crs as ccrs
+
+grid_mapping = rotds.cf["grid_mapping"]
+pole_latitude = grid_mapping.grid_north_pole_latitude
+pole_longitude = grid_mapping.grid_north_pole_longitude
+
+ccrs.RotatedPole(pole_longitude, pole_latitude)
+```
+
+![cartopy rotated pole projection](cartopy_rotated_pole.png)
diff --git a/doc/index.rst b/doc/index.rst
index b5c14bac..1e26ea52 100644
--- a/doc/index.rst
+++ b/doc/index.rst
@@ -51,6 +51,7 @@ or using ``conda``
    units
    parametricz
    bounds
+   grid_mappings
    coding
    dsg
    sgrid_ugrid
diff --git a/doc/quickstart.md b/doc/quickstart.md
index eb8ebdc0..b86546ce 100644
--- a/doc/quickstart.md
+++ b/doc/quickstart.md
@@ -91,7 +91,8 @@ Sometimes it is more useful to extract the actual variable names associated with
 - {py:attr}`Dataset.cf.cell_measures`,
 - {py:attr}`Dataset.cf.cf_roles`,
 - {py:attr}`Dataset.cf.coordinates`,
-- {py:attr}`Dataset.cf.formula_terms`, and
+- {py:attr}`Dataset.cf.formula_terms`,
+- {py:attr}`Dataset.cf.grid_mapping_names`, and
 - {py:attr}`Dataset.cf.standard_names`.
 
 These properties all return dictionaries mapping a standard key name to a list of matching variable names in the Dataset or DataArray.
diff --git a/doc/whats-new.rst b/doc/whats-new.rst
index 629d2f55..4d753f23 100644
--- a/doc/whats-new.rst
+++ b/doc/whats-new.rst
@@ -3,8 +3,12 @@
 What's New
 ----------
 
-v (unreleased)
-==============
+v0.8.0 (unreleased)
+===================
+
+- Support interpreting the ``grid_mapping`` attribute: (:pr:`391`).
+  See :py:meth:`Dataset.cf.grid_mapping_names` and ``Dataset.cf["grid_mapping"]``,
+  ``DataArray.cf["grid_mapping"]``. By `Lars Buntemeyer`_ and `Deepak Cherian`_.
 
 
 v0.7.9 (Jan 31, 2023)