Repository: tomvothecoder/xcdat

Issue:
Add fetching of dict of bounds for a DataArray, similar to cf-xarray with Datasets
### Is your feature request related to a problem? If yes, Please describe
<!-- Concise description of the problem -->

`cf-xarray` only supports retrieving bounds for Datasets. `xcdat` adds functionality to store bounds in DataArrays, but needs a way to retrieve bounds dynamically.

### Describe your proposed solution
<!-- Concise description of your proposed solution. -->
Base the DataArray implementation of retrieving bounds similar to `cf-xarray` with `Datasets`.

[bounds](https://github.com/xarray-contrib/cf-xarray/blob/5523a8f09ce00f486d06e51fdda60a303c5b844e/cf_xarray/accessor.py#L1876)
[get_bounds](https://github.com/xarray-contrib/cf-xarray/blob/5523a8f09ce00f486d06e51fdda60a303c5b844e/cf_xarray/accessor.py#L1912)
[get_bounds_dim_name](https://github.com/xarray-contrib/cf-xarray/blob/5523a8f09ce00f486d06e51fdda60a303c5b844e/cf_xarray/accessor.py#L323)

### Describe alternatives solutions you've considered
<!-- Concise description of any alternative solutions or features you've considered. -->

### Additional context
<!-- Other context, examples, or screenshots about the feature request here. -->


Discussion:
codecov[bot]: # [Codecov](https://codecov.io/gh/tomvothecoder/xcdat/pull/92?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Tom+Vo) Report
> Merging [#92](https://codecov.io/gh/tomvothecoder/xcdat/pull/92?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Tom+Vo) (64a044b) into [main](https://codecov.io/gh/tomvothecoder/xcdat/commit/5918313acb11ed3f124711863561fa0da694995a?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Tom+Vo) (5918313) will **not change** coverage.
> The diff coverage is `100.00%`.

[![Impacted file tree graph](https://codecov.io/gh/tomvothecoder/xcdat/pull/92/graphs/tree.svg?width=650&height=150&src=pr&token=T8W41LN26O&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Tom+Vo)](https://codecov.io/gh/tomvothecoder/xcdat/pull/92?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Tom+Vo)

```diff
@@            Coverage Diff            @@
##              main       #92   +/-   ##
=========================================
  Coverage   100.00%   100.00%           
=========================================
  Files            5         5           
  Lines          138       165   +27     
=========================================
+ Hits           138       165   +27     
```


| [Impacted Files](https://codecov.io/gh/tomvothecoder/xcdat/pull/92?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Tom+Vo) | Coverage Δ | |
|---|---|---|
| [xcdat/dataset.py](https://codecov.io/gh/tomvothecoder/xcdat/pull/92/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Tom+Vo#diff-eGNkYXQvZGF0YXNldC5weQ==) | `100.00% <ø> (ø)` | |
| [xcdat/bounds.py](https://codecov.io/gh/tomvothecoder/xcdat/pull/92/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Tom+Vo#diff-eGNkYXQvYm91bmRzLnB5) | `100.00% <100.00%> (ø)` | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/tomvothecoder/xcdat/pull/92?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Tom+Vo).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Tom+Vo)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/tomvothecoder/xcdat/pull/92?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Tom+Vo). Last update [5918313...64a044b](https://codecov.io/gh/tomvothecoder/xcdat/pull/92?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Tom+Vo). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Tom+Vo).


Patch:
diff --git a/xcdat/bounds.py b/xcdat/bounds.py
index 47951ec..197c74f 100644
--- a/xcdat/bounds.py
+++ b/xcdat/bounds.py
@@ -1,9 +1,10 @@
 """Bounds module for functions related to coordinate bounds."""
-from typing import Optional, Tuple, get_args
+from typing import Dict, Hashable, List, Optional, Tuple, get_args
 
 import cf_xarray as cfxr  # noqa: F401
 import numpy as np
 import xarray as xr
+from cf_xarray.accessor import _get_bounds, _single, _variables, apply_mapper
 from typing_extensions import Literal
 
 from xcdat.logger import setup_custom_logger
@@ -61,9 +62,7 @@ class DatasetBoundsAccessor:
         allow_generating: bool = True,
         width: float = 0.5,
     ) -> Optional[xr.DataArray]:
-        """Get bounds for a coordinate.
-
-        If bounds don't exist, they will be generated and set in the dataset.
+        """Get bounds for a coordinate, or generate bounds if it doesn't exist.
 
         Parameters
         ----------
@@ -181,14 +180,15 @@ class DatasetBoundsAccessor:
 
         # Add coordinate bounds to the dataset
         var_name = f"{coord}_bnds"
-        self._dataset[da_coord.name] = da_coord.assign_attrs(bounds=var_name)
         self._dataset[var_name] = xr.DataArray(
             name=var_name,
             data=bounds,
             coords={coord: da_coord},
             dims=[coord, "bnds"],
-            attrs={**da_coord.attrs, "is_generated": "True"},
+            attrs={"is_generated": "True"},
         )
+        self._dataset[da_coord.name].attrs["bounds"] = var_name
+
         return self._dataset[var_name]
 
     def _get_coord(self, coord: Coord) -> xr.DataArray:
@@ -271,6 +271,10 @@ class DataArrayBoundsAccessor:
 
         """
         da = self._dataarray.copy()
+
+        # The bounds dimension must be set before adding bounds to the DataArray
+        # coordinates, otherwise the error below is thrown:
+        # "ValueError: cannot add coordinates with new dimensions to a DataArray"
         da = self._set_bounds_dim(dataset)
 
         coords = [*dataset.coords]
@@ -294,11 +298,12 @@ class DataArrayBoundsAccessor:
         return self._dataarray
 
     def _set_bounds_dim(self, dataset: xr.Dataset) -> xr.DataArray:
-        """Sets the bounds dimension in the DataArray using the parent Dataset.
+        """
+        Sets the bounds dimension(s) in the DataArray based on the dims of the
+        parent Dataset.
 
-        The bounds dimension must be set before adding bounds to the DataArray
-        coordinates, otherwise the error below will be thrown:
-        ``ValueError: cannot add coordinates with new dimensions to a DataArray``.
+        This function uses the "bounds" attribute of each coordinate to map to
+        the bounds, then extracts the dimension from each bounds.
 
         Parameters
         ----------
@@ -316,20 +321,118 @@ class DataArrayBoundsAccessor:
             When no bounds dimension exists in the parent Dataset.
         """
         da = self._dataarray.copy()
-        dims = dataset.dims.keys()
+        coords = dataset.cf.coordinates.keys()
+
+        dims = set()
+        for coord in coords:
+            try:
+                dims.add(dataset.cf.get_bounds_dim_name(coord))
+            except KeyError:
+                logger.warning(
+                    f"{coord} has no bounds, or the `bounds` attribute is missing to "
+                    "link to the bounds."
+                )
 
-        if "bnds" in dims:
-            bounds_dim = "bnds"
-        elif "bounds" in dims:
-            bounds_dim = "bounds"
-        else:
+        if len(dims) == 0:
             raise KeyError(
-                "No bounds dimension in the parent dataset. This indicates that there "
-                "are probably no coordinate bounds in the dataset. Pass the "
-                "dataset to `xcdat.dataset.open_dataset` to auto-generate bounds."
+                "No bounds dimension in the parent dataset, which indicates that there "
+                "are probably no coordinate bounds. Try passing the dataset to "
+                "`xcdat.dataset.open_dataset` to auto-generate them."
             )
 
-        da = da.expand_dims(dim={bounds_dim: np.array([0, 1])})
+        for dim in dims:
+            da = da.expand_dims(dim={dim: np.array([0, 1])})
 
         self._dataarray = da
         return self._dataarray
+
+    @property
+    def bounds(self) -> Dict[str, Optional[xr.DataArray]]:
+        """Returns a mapping of coordinate keys to their bounds.
+
+        Missing coordinates are handled by ``self.copy_from_parent()``.
+
+        Returns
+        -------
+        Dict[str, Optional[xr.DataArray]]
+            Dictionary mapping coordinate keys to their bounds.
+
+        Notes
+        -----
+        Based on ``cf_xarray.accessor.CFDatasetAccessor.bounds``.
+        """
+        da = self._dataarray
+        bounds: Dict[str, Optional[xr.DataArray]] = {}
+        for coord, bounds_name in self.bounds_names.items():
+            bounds[coord] = da.coords.get(bounds_name, None)  # type: ignore
+
+        return bounds
+
+    @property
+    def bounds_names(self) -> Dict[Hashable, List[str]]:
+        """Returns a mapping of coordinate keys to the name of their bounds.
+
+        Missing coordinates are handled by ``self.copy_from_parent()``.
+
+        Returns
+        -------
+        Dict[Hashable, List[str]]
+            Dictionary mapping valid keys to the variable names of their bounds.
+
+        Notes
+        -----
+        Based on ``cf_xarray.accessor.CFDatasetAccessor.bounds_names``.
+        """
+        da = self._dataarray
+        keys = da.coords
+
+        vardict = {key: apply_mapper(_get_bounds, da, key, error=False) for key in keys}
+
+        # Each coord should have only one bound, thus select the first index
+        # TODO: Handle when there is more than one bounds per coordinate.
+        return {k: v[0] for k, v in vardict.items() if v}
+
+    def get_bounds(self, coord: Coord) -> xr.DataArray:
+        """Get bounds corresponding to a coordinate key.
+
+        Parameters
+        ----------
+        coord : Coord
+            Name of the coordinate whose bounds are desired.
+
+        Returns
+        -------
+        DataArray
+            The bounds for a coordinate key.
+
+        Notes
+        -----
+        Based on ``cf_xarray.accessor.CFDatasetAccessor.get_bounds``.
+        """
+        # TODO: Handle when there is more than bounds dimension
+        return apply_mapper(_variables(_single(_get_bounds)), self._dataarray, coord)[0]
+
+    def get_bounds_dim_name(self, coord: Coord) -> Hashable:
+        """Get bounds dimension name corresponding to coordinate key.
+
+        Parameters
+        ----------
+        coord : Coord
+            Name of the coordinate whose bounds dimension is desired.
+
+        Returns
+        -------
+        Hashable
+            The bounds dimension name.
+
+        Notes
+        -----
+        Based on ``cf_xarray.accessor.CFDatasetAccessor.get_bounds_dim_name``.
+        """
+        crd = self._dataarray[coord]
+        bounds = self.get_bounds(coord)
+        bounds_dims = set(bounds.dims) - set(crd.dims)
+        assert len(bounds_dims) == 1
+        bounds_dim = bounds_dims.pop()
+        assert self._dataarray.sizes[bounds_dim] in [2, 4]
+        return bounds_dim
diff --git a/xcdat/dataset.py b/xcdat/dataset.py
index e8c2752..9391ada 100644
--- a/xcdat/dataset.py
+++ b/xcdat/dataset.py
@@ -15,6 +15,9 @@ def open_dataset(path: str, **kwargs: Dict[str, Any]) -> xr.Dataset:
     - Decode all time units, including non-CF compliant units (months and years).
     - Generate bounds for supported coordinates if they don't exist.
 
+    NOTE: Using decode_times=False may add incorrect units for existing time
+    bounds (becomes "days since 1970-01-01  00:00:00").
+
     Parameters
     ----------
     path : str
@@ -65,6 +68,9 @@ def open_mfdataset(paths: Union[str, List[str]], **kwargs) -> xr.Dataset:
     - Decode all time units, including non-CF compliant units (months and years).
     - Generate bounds for supported coordinates if they don't exist.
 
+    NOTE: Using decode_times=False may add incorrect units for existing time
+    bounds (becomes "days since 1970-01-01  00:00:00").
+
     Parameters
     ----------
     path : Union[str, List[str]]
@@ -194,11 +200,10 @@ def decode_time_units(dataset: xr.Dataset):
     cf_compliant = units not in non_cf_units_to_freq.keys()
     if cf_compliant:
         dataset = xr.decode_cf(dataset, decode_times=True)
-        dataset.time.attrs["units"] = units_attr
     else:
-        # NOTE: Calendar type for "months" or years" does not matter because the
-        # number of days in a month does not play a factor when generating date
-        # ranges for these units.
+        # NOTE: The "calendar" attribute for units consisting of "months" or
+        # "years" is not factored when generating date ranges. The number of
+        # days in a month is not factored.
         decoded_time = xr.DataArray(
             data=pd.date_range(
                 start=reference_date,
@@ -214,7 +219,7 @@ def decode_time_units(dataset: xr.Dataset):
             "original_shape": decoded_time.shape,
             "units": units_attr,
             # pandas.date_range() returns "proleptic_gregorian" by default
-            "calendar": time.attrs.get("calendar", "proleptic_gregorian"),
+            "calendar": "proleptic_gregorian",
         }
 
         dataset = dataset.assign_coords({"time": decoded_time})


Test patch:
diff --git a/tests/fixtures.py b/tests/fixtures.py
index c93d59c..838268e 100644
--- a/tests/fixtures.py
+++ b/tests/fixtures.py
@@ -9,8 +9,13 @@ import xarray as xr
 # side-effects caused by reference assignment.
 # https://xarray.pydata.org/en/stable/generated/xarray.DataArray.copy.html
 
-# Dataset coordinates
-time = xr.DataArray(
+# NOTE:
+# - Non-CF time includes "units" attr
+# - Coordinates with bounds includes "bounds" attr and vice versa
+
+# TIME
+# ====
+time_cf = xr.DataArray(
     data=[
         datetime(2000, 1, 1),
         datetime(2000, 2, 1),
@@ -26,25 +31,22 @@ time = xr.DataArray(
         datetime(2000, 12, 1),
     ],
     dims=["time"],
+    attrs={
+        "long_name": "time",
+        "standard_name": "time",
+        "axis": "T",
+    },
 )
-time_non_cf_compliant = xr.DataArray(
+time_non_cf = xr.DataArray(
     data=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11],
     dims=["time"],
-    attrs={"units": "months since 2000-01-01"},
-)
-
-lat = xr.DataArray(
-    data=np.array([-90, -88.75, 88.75, 90]),
-    dims=["lat"],
-    attrs={"units": "degrees_north", "axis": "Y"},
-)
-lon = xr.DataArray(
-    data=np.array([0, 1.875, 356.25, 358.125]),
-    dims=["lon"],
-    attrs={"units": "degrees_east", "axis": "X"},
+    attrs={
+        "long_name": "time",
+        "standard_name": "time",
+        "axis": "T",
+    },
 )
 
-# Dataset data variables (bounds)
 time_bnds = xr.DataArray(
     name="time_bnds",
     data=[
@@ -61,12 +63,13 @@ time_bnds = xr.DataArray(
         [datetime(2000, 10, 16, 12), datetime(2000, 11, 16)],
         [datetime(2000, 11, 16), datetime(2000, 12, 16)],
     ],
-    coords={"time": time},
+    coords={"time": time_cf},
     dims=["time", "bnds"],
-    attrs={"is_generated": "True"},
+    attrs={
+        "is_generated": "True",
+    },
 )
-
-time_bnds_non_cf_compliant = xr.DataArray(
+time_bnds_non_cf = xr.DataArray(
     name="time_bnds",
     data=[
         [datetime(1999, 12, 16, 12), datetime(2000, 1, 16, 12)],
@@ -82,16 +85,33 @@ time_bnds_non_cf_compliant = xr.DataArray(
         [datetime(2000, 10, 16, 12), datetime(2000, 11, 16)],
         [datetime(2000, 11, 16), datetime(2000, 12, 16)],
     ],
-    coords={"time": time.data},
+    coords={"time": time_non_cf},
     dims=["time", "bnds"],
     attrs={"is_generated": "True"},
 )
+
+# LATITUDE
+# ========
+lat = xr.DataArray(
+    data=np.array([-90, -88.75, 88.75, 90]),
+    dims=["lat"],
+    attrs={"units": "degrees_north", "axis": "Y"},
+)
 lat_bnds = xr.DataArray(
     name="lat_bnds",
     data=np.array([[-90, -89.375], [-89.375, 0.0], [0.0, 89.375], [89.375, 90]]),
     coords={"lat": lat.data},
     dims=["lat", "bnds"],
-    attrs={"units": "degrees_north", "axis": "Y", "is_generated": "True"},
+    attrs={"is_generated": "True"},
+)
+
+
+# LONGITUDE
+# =========
+lon = xr.DataArray(
+    data=np.array([0, 1.875, 356.25, 358.125]),
+    dims=["lon"],
+    attrs={"units": "degrees_east", "axis": "X"},
 )
 lon_bnds = xr.DataArray(
     name="lon_bnds",
@@ -105,84 +125,108 @@ lon_bnds = xr.DataArray(
     ),
     coords={"lon": lon.data},
     dims=["lon", "bnds"],
-    attrs={"units": "degrees_east", "axis": "X", "is_generated": "True"},
+    attrs={"is_generated": "True"},
 )
 
-# Dataset data variables (variables)
-ts = xr.DataArray(
+# VARIABLES
+# =========
+ts_cf = xr.DataArray(
     name="ts",
     data=np.ones((12, 4, 4)),
-    coords={"time": time, "lat": lat, "lon": lon},
+    coords={"time": time_cf, "lat": lat, "lon": lon},
     dims=["time", "lat", "lon"],
 )
-ts_non_cf_compliant = xr.DataArray(
+
+ts_non_cf = xr.DataArray(
     name="ts",
     data=np.ones((12, 4, 4)),
-    coords={"time": time_non_cf_compliant, "lat": lat, "lon": lon},
+    coords={"time": time_non_cf, "lat": lat, "lon": lon},
     dims=["time", "lat", "lon"],
 )
 
-ts_with_bnds = xr.DataArray(
+# Special cases, after copying bounds from parent Dataset to data variable
+# using xCDAT.
+ts_with_bnds_from_parent_cf = xr.DataArray(
     name="ts",
     data=np.ones((2, 12, 4, 4)),
     coords={
         "bnds": np.array([0, 1]),
-        "time": time.assign_attrs(bounds="time_bnds"),
+        "time": time_cf.assign_attrs(bounds="time_bnds"),
         "lat": lat.assign_attrs(bounds="lat_bnds"),
         "lon": lon.assign_attrs(bounds="lon_bnds"),
+        "lat_bnds": lat_bnds.copy(),
+        "lon_bnds": lon_bnds.copy(),
+        "time_bnds": time_bnds.copy(),
+    },
+    dims=["bnds", "time", "lat", "lon"],
+)
+ts_with_bnds_from_parent_non_cf = xr.DataArray(
+    name="ts",
+    data=np.ones((2, 12, 4, 4)),
+    coords={
+        "bnds": np.array([0, 1]),
+        "time": time_cf.copy(),
+        "lat": lat.copy(),
+        "lon": lon.copy(),
         "lat_bnds": lat_bnds,
         "lon_bnds": lon_bnds,
-        "time_bnds": time_bnds,
+        "time_bnds": time_bnds_non_cf,
     },
-    dims=[
-        "bnds",
-        "time",
-        "lat",
-        "lon",
-    ],
+    dims=["bnds", "time", "lat", "lon"],
 )
 
 
-def generate_dataset(cf_compliant=True, has_bounds: bool = True) -> xr.Dataset:
+def generate_dataset(cf_compliant: bool, has_bounds: bool) -> xr.Dataset:
     """Generates a dataset using coordinate and data variable fixtures.
 
-    NOTE: Using ``.assign()`` to add data variables to an existing dataset will
-    remove attributes from existing coordinates. The workaround is to update a
-    data_vars dict then create the dataset. https://github.com/pydata/xarray/issues/2245
-
     Parameters
     ----------
     cf_compliant : bool, optional
-        CF compliant time units, by default True
+        CF compliant time units.
     has_bounds : bool, optional
-        Include bounds for coordinates, by default True
+        Include bounds for coordinates. This also adds the "bounds" attribute
+        to existing coordinates to link them to their respective bounds.
 
     Returns
     -------
     xr.Dataset
         Test dataset.
     """
-    data_vars = {}
-    coords = {
-        "lat": lat.copy(),
-        "lon": lon.copy(),
-    }
-
-    if cf_compliant:
-        coords.update({"time": time.copy()})
-        data_vars.update({"ts": ts.copy()})
-    else:
-        coords.update({"time": time_non_cf_compliant.copy()})
-        data_vars.update({"ts": ts_non_cf_compliant.copy()})
-
     if has_bounds:
-        data_vars.update(
-            {
-                "time_bnds": time_bnds.copy(),
+        ds = xr.Dataset(
+            data_vars={
+                "ts": ts_cf.copy(),
                 "lat_bnds": lat_bnds.copy(),
                 "lon_bnds": lon_bnds.copy(),
-            }
+            },
+            coords={"lat": lat.copy(), "lon": lon.copy()},
         )
 
-    ds = xr.Dataset(data_vars=data_vars, coords=coords)
+        if cf_compliant:
+            ds = ds.assign({"time_bnds": time_bnds.copy()})
+            ds = ds.assign_coords({"time": time_cf.copy()})
+        elif not cf_compliant:
+            ds = ds.assign({"time_bnds": time_bnds_non_cf.copy()})
+            ds = ds.assign_coords({"time": time_non_cf.copy()})
+            ds["time"] = ds.time.assign_attrs(units="months since 2000-01-01")
+
+        # If the "bounds" attribute is included in an existing DataArray and
+        # added to a new Dataset, it will get dropped. Therefore, it needs to be
+        # assigned to the DataArrays after they are added to Dataset.
+        ds["lat"] = ds.lat.assign_attrs(bounds="lat_bnds")
+        ds["lon"] = ds.lon.assign_attrs(bounds="lon_bnds")
+        ds["time"] = ds.time.assign_attrs(bounds="time_bnds")
+
+    elif not has_bounds:
+        ds = xr.Dataset(
+            data_vars={"ts": ts_cf.copy()},
+            coords={"lat": lat.copy(), "lon": lon.copy()},
+        )
+
+        if cf_compliant:
+            ds = ds.assign_coords({"time": time_cf.copy()})
+        elif not cf_compliant:
+            ds = ds.assign_coords({"time": time_non_cf.copy()})
+            ds["time"] = ds.time.assign_attrs(units="months since 2000-01-01")
+
     return ds
diff --git a/tests/test_bounds.py b/tests/test_bounds.py
index 2b92fa5..238f76b 100644
--- a/tests/test_bounds.py
+++ b/tests/test_bounds.py
@@ -2,67 +2,14 @@ import numpy as np
 import pytest
 import xarray as xr
 
-from tests.fixtures import generate_dataset, ts_with_bnds
+from tests.fixtures import generate_dataset, ts_cf, ts_with_bnds_from_parent_cf
 from xcdat.bounds import DataArrayBoundsAccessor, DatasetBoundsAccessor
 
 
 class TestDatasetBoundsAccessor:
     @pytest.fixture(autouse=True)
     def setup(self):
-        # Coordinate information
-        lat = xr.DataArray(
-            data=np.array([-90, -88.75, 0, 88.75, 90]),
-            dims=["lat"],
-            attrs={"units": "degrees_north", "axis": "Y"},
-        )
-        lon = xr.DataArray(
-            data=np.array([0, 1.875, 178.125, 356.25, 358.125]),
-            dims=["lon"],
-            attrs={"units": "degrees_east", "axis": "X"},
-        )
-
-        # Generated coordinates bounds
-        self.lat_bnds = xr.DataArray(
-            name="lat_bnds",
-            data=np.array(
-                [
-                    [-90, -89.375],
-                    [-89.375, -44.375],
-                    [-44.375, 44.375],
-                    [44.375, 89.375],
-                    [89.375, 90],
-                ]
-            ),
-            coords={"lat": lat},
-            dims=["lat", "bnds"],
-            attrs={
-                "units": "degrees_north",
-                "axis": "Y",
-                "is_generated": "True",
-            },
-        )
-        self.lon_bnds = xr.DataArray(
-            name="lon_bnds",
-            data=np.array(
-                [
-                    [-0.9375, 0.9375],
-                    [0.9375, 90.0],
-                    [90.0, 267.1875],
-                    [267.1875, 357.1875],
-                    [357.1875, 359.0625],
-                ]
-            ),
-            coords={"lon": lon},
-            dims=["lon", "bnds"],
-            attrs={
-                "units": "degrees_east",
-                "axis": "X",
-                "is_generated": "True",
-            },
-        )
-
-        # Create Dataset using coordinates
-        self.ds = xr.Dataset(coords={"lat": lat, "lon": lon})
+        self.ds = generate_dataset(cf_compliant=True, has_bounds=False)
 
     def test__init__(self):
         obj = DatasetBoundsAccessor(self.ds)
@@ -72,18 +19,12 @@ class TestDatasetBoundsAccessor:
         assert self.ds.bounds._dataset.identical(self.ds)
 
     def test_get_bounds_when_bounds_exist_in_dataset(self):
-        ds = self.ds.copy()
-        ds = ds.assign(
-            lat_bnds=self.lat_bnds,
-            lon_bnds=self.lon_bnds,
-        )
-
+        ds = generate_dataset(cf_compliant=True, has_bounds=True)
         lat_bnds = ds.bounds.get_bounds("lat")
-        assert lat_bnds is not None and lat_bnds.identical(self.lat_bnds)
-        assert lat_bnds.is_generated
+        assert lat_bnds.identical(ds.lat_bnds)
 
         lon_bnds = ds.bounds.get_bounds("lon")
-        assert lon_bnds is not None and lon_bnds.identical(self.lon_bnds)
+        assert lon_bnds.identical(ds.lon_bnds)
         assert lon_bnds.is_generated
 
     def test_get_bounds_when_bounds_do_not_exist_in_dataset(self):
@@ -91,12 +32,12 @@ class TestDatasetBoundsAccessor:
 
         lat_bnds = ds.bounds.get_bounds("lat")
         assert lat_bnds is not None
-        assert lat_bnds.identical(self.lat_bnds)
+        assert lat_bnds.identical(ds.lat_bnds)
         assert lat_bnds.is_generated
 
         lon_bnds = ds.bounds.get_bounds("lon")
         assert lon_bnds is not None
-        assert lon_bnds.identical(self.lon_bnds)
+        assert lon_bnds.identical(ds.lon_bnds)
         assert lon_bnds.is_generated
 
         # Check raises error when bounds do not exist and not allowing generated bounds.
@@ -105,21 +46,21 @@ class TestDatasetBoundsAccessor:
             ds.bounds.get_bounds("lat", allow_generating=False)
 
     def test_get_bounds_raises_error_with_incorrect_axis_argument(self):
-        ds = self.ds.copy()
-
         with pytest.raises(ValueError):
-            ds.bounds.get_bounds("incorrect_axis_argument")
+            self.ds.bounds.get_bounds("incorrect_axis_argument")
 
     def test__get_bounds_does_not_drop_attrs_of_existing_coords_when_generating_bounds(
         self,
     ):
-        ds = self.ds.copy()
+        ds_expected = self.ds.copy()
+        ds_expected["lat"].attrs["bounds"] = "lat_bnds"
 
-        lat_bnds = ds.bounds.get_bounds("lat", allow_generating=True)
-        assert lat_bnds.identical(self.lat_bnds)
+        ds_result = self.ds.copy()
+        lat_bnds = ds_result.bounds.get_bounds("lat", allow_generating=True)
+        assert lat_bnds.identical(ds_result.lat_bnds)
 
-        ds = ds.drop_vars("lat_bnds")
-        assert ds.identical(self.ds)
+        ds_result = ds_result.drop_vars("lat_bnds")
+        assert ds_result.identical(ds_expected)
 
     def test__generate_bounds_raises_errors_for_data_dim_and_length(self):
         # Multidimensional
@@ -147,11 +88,11 @@ class TestDatasetBoundsAccessor:
         ds = self.ds.copy()
 
         lat_bnds = ds.bounds._generate_bounds("lat")
-        assert lat_bnds.equals(self.lat_bnds)
+        assert lat_bnds.equals(ds.lat_bnds)
         assert ds.lat_bnds.is_generated
 
         lon_bnds = ds.bounds._generate_bounds("lon")
-        assert lon_bnds.equals(self.lon_bnds)
+        assert lon_bnds.equals(ds.lon_bnds)
         assert ds.lon_bnds.is_generated
 
     def test__get_coord(self):
@@ -168,18 +109,15 @@ class TestDatasetBoundsAccessor:
     def test__get_coord_raises_error_if_coord_does_not_exist(self):
         ds = self.ds.copy()
 
+        ds = ds.drop_dims("lat")
         with pytest.raises(KeyError):
-            ds = ds.drop_vars("lat")
             ds.bounds._get_coord("lat")
 
 
 class TestDataArrayBoundsAccessor:
     @pytest.fixture(autouse=True)
     def setUp(self):
-        self.ds = generate_dataset(has_bounds=True)
-        self.ds.lat.attrs["bounds"] = "lat_bnds"
-        self.ds.lon.attrs["bounds"] = "lon_bnds"
-        self.ds.time.attrs["bounds"] = "time_bnds"
+        self.ds = generate_dataset(cf_compliant=True, has_bounds=True)
 
     def test__init__(self):
         obj = DataArrayBoundsAccessor(self.ds.ts)
@@ -191,27 +129,64 @@ class TestDataArrayBoundsAccessor:
         assert result.identical(expected)
 
     def test__copy_from_parent_copies_bounds(self):
-        ts_expected = ts_with_bnds.copy()
-
-        ts = self.ds["ts"].copy()
-        ts_result = ts.bounds.copy_from_parent(self.ds)
+        ts_expected = ts_with_bnds_from_parent_cf.copy()
+        ts_result = ts_cf.copy().bounds.copy_from_parent(self.ds)
 
         assert ts_result.identical(ts_expected)
 
     def test__set_bounds_dim_adds_bnds(self):
-        ts = self.ds["ts"].copy()
-        ts_expected = ts.copy()
-        ts_expected = ts_expected.expand_dims(bnds=np.array([0, 1]))
-        ts_result = ts.bounds._set_bounds_dim(self.ds)
+        ts_expected = ts_cf.copy().expand_dims(bnds=np.array([0, 1]))
+        ts_result = ts_cf.copy().bounds._set_bounds_dim(self.ds)
 
         assert ts_result.identical(ts_expected)
 
     def test__set_bounds_dim_adds_bounds(self):
         ds = self.ds.swap_dims({"bnds": "bounds"}).copy()
-        ts = ds["ts"].copy()
 
-        ts_expected = ts.copy()
-        ts_expected = ts_expected.expand_dims(bounds=np.array([0, 1]))
-
-        ts_result = ts.bounds._set_bounds_dim(ds)
+        ts_expected = ts_cf.copy().expand_dims(bounds=np.array([0, 1]))
+        ts_result = ts_cf.copy().bounds._set_bounds_dim(ds)
         assert ts_result.identical(ts_expected)
+
+    def test_bounds_property_returns_mapping_of_coordinate_keys_to_bounds_dataarrays(
+        self,
+    ):
+        ts = ts_with_bnds_from_parent_cf.copy()
+        result = ts.bounds.bounds
+        expected = {"time": ts.time_bnds, "lat": ts.lat_bnds, "lon": ts.lon_bnds}
+
+        for key in expected.keys():
+            assert result[key].identical(expected[key])
+
+    def test_bounds_names_property_returns_mapping_of_coordinate_keys_to_names_of_bounds(
+        self,
+    ):
+        ts = ts_with_bnds_from_parent_cf.copy()
+        result = ts.bounds.bounds_names
+        expected = {"time": "time_bnds", "lat": "lat_bnds", "lon": "lon_bnds"}
+        assert result == expected
+
+    def test_get_bounds_returns_bounds_dataarray_for_coordinate_key(self):
+        ts = ts_with_bnds_from_parent_cf.copy()
+        result = ts.bounds.get_bounds("lat")
+        expected = ts.lat_bnds
+
+        assert result.identical(expected)
+
+    def test_get_bounds_returns_keyerror_with_incorrect_key(self):
+        ts = ts_with_bnds_from_parent_cf.copy()
+
+        with pytest.raises(KeyError):
+            ts.bounds.get_bounds("something")
+
+    def test_get_bounds_dim_name_returns_dim_of_coordinate_bounds(self):
+        ts = ts_with_bnds_from_parent_cf.copy()
+        result = ts.bounds.get_bounds_dim_name("lat")
+        expected = "bnds"
+
+        assert result == expected
+
+    def test_get_bounds_dim_name_raises_keyerror_with_incorrect_key(self):
+        ts = ts_with_bnds_from_parent_cf.copy()
+
+        with pytest.raises(KeyError):
+            ts.bounds.get_bounds_dim_name("something")
diff --git a/tests/test_dataset.py b/tests/test_dataset.py
index 96d32cf..5895bc5 100644
--- a/tests/test_dataset.py
+++ b/tests/test_dataset.py
@@ -17,12 +17,19 @@ class TestOpenDataset:
         self.file_path = f"{self.dir}/file.nc"
 
     def test_non_cf_compliant_time_is_decoded(self):
-        # Generate two dummy datasets with non-CF compliant time units.
+        # Generate dummy datasets with non-CF compliant time units that aren't
+        # encoded yet.
         ds = generate_dataset(cf_compliant=False, has_bounds=False)
         ds.to_netcdf(self.file_path)
 
-        # Generate an expected dataset, which is a combination of both datasets
-        # with decoded time units and coordinate bounds.
+        result_ds = open_dataset(self.file_path)
+        # Replicates decode_times=False, which adds units to "time" coordinate.
+        # Refer to xcdat.bounds.DatasetBoundsAccessor._generate_bounds() for
+        # how attributes propagate from coord to coord bounds.
+        result_ds["time_bnds"].attrs["units"] = "months since 2000-01-01"
+
+        # Generate an expected dataset with non-CF compliant time units that are
+        # manually encoded
         expected_ds = generate_dataset(cf_compliant=True, has_bounds=True)
         expected_ds.time.attrs["units"] = "months since 2000-01-01"
         expected_ds.time_bnds.attrs["units"] = "months since 2000-01-01"
@@ -35,26 +42,20 @@ class TestOpenDataset:
         }
 
         # Check that non-cf compliant time was decoded and bounds were generated.
-        result_ds = open_dataset(self.file_path)
         assert result_ds.identical(expected_ds)
 
     def test_preserves_lat_and_lon_bounds_if_they_exist(self):
         # Create expected dataset which includes bounds.
-        expected_ds = generate_dataset(has_bounds=True)
+        expected_ds = generate_dataset(cf_compliant=True, has_bounds=True)
         expected_ds.to_netcdf(self.file_path)
 
-        # Need to add time units since decode time units is also called,
-        # which adds this attribute.
-        expected_ds.time.attrs["units"] = "days since 2000-01-01 00:00:00"
-        expected_ds.time_bnds.attrs["units"] = "days since 2000-01-01 00:00:00"
-
         # Check resulting dataset and expected are identical
         result_ds = open_dataset(self.file_path)
         assert result_ds.identical(expected_ds)
 
     def test_generates_lat_and_lon_bounds_if_they_dont_exist(self):
         # Create expected dataset without bounds.
-        ds = generate_dataset(has_bounds=False)
+        ds = generate_dataset(cf_compliant=True, has_bounds=False)
         ds.to_netcdf(self.file_path)
 
         # Make sure bounds don't exist
@@ -84,11 +85,16 @@ class TestOpenMfDataset:
         # Generate two dummy datasets with non-CF compliant time units.
         ds1 = generate_dataset(cf_compliant=False, has_bounds=False)
         ds1.to_netcdf(self.file_path1)
-
         ds2 = generate_dataset(cf_compliant=False, has_bounds=False)
         ds2 = ds2.rename_vars({"ts": "tas"})
         ds2.to_netcdf(self.file_path2)
 
+        result_ds = open_mfdataset([self.file_path1, self.file_path2])
+        # Replicates decode_times=False, which adds units to "time" coordinate.
+        # Refer to xcdat.bounds.DatasetBoundsAccessor._generate_bounds() for
+        # how attributes propagate from coord to coord bounds.
+        result_ds.time_bnds.attrs["units"] = "months since 2000-01-01"
+
         # Generate an expected dataset, which is a combination of both datasets
         # with decoded time units and coordinate bounds.
         expected_ds = generate_dataset(cf_compliant=True, has_bounds=True)
@@ -104,7 +110,6 @@ class TestOpenMfDataset:
         }
 
         # Check that non-cf compliant time was decoded and bounds were generated.
-        result_ds = open_mfdataset([self.file_path1, self.file_path2])
         assert result_ds.identical(expected_ds)
 
     def test_preserves_lat_and_lon_bounds_if_they_exist(self):
@@ -117,10 +122,8 @@ class TestOpenMfDataset:
         ds2.to_netcdf(self.file_path2)
 
         # Generate expected dataset, which is a combination of the two datasets.
-        expected_ds = generate_dataset(has_bounds=True)
+        expected_ds = generate_dataset(cf_compliant=True, has_bounds=True)
         expected_ds["tas"] = expected_ds["ts"].copy()
-        expected_ds.time.attrs["units"] = "days since 2000-01-01 00:00:00"
-        expected_ds.time_bnds.attrs["units"] = "days since 2000-01-01 00:00:00"
 
         # Check that the result is identical to the expected.
         result_ds = open_mfdataset([self.file_path1, self.file_path2])
@@ -155,7 +158,6 @@ class TestDecodeTimeUnits:
         # on the unit that needs to be tested (days (CF compliant) or months
         # (non-CF compliant).
         self.time_attrs = {
-            "units": None,
             "bounds": "time_bnds",
             "axis": "T",
             "long_name": "time",
@@ -165,7 +167,7 @@ class TestDecodeTimeUnits:
     def test_throws_error_if_function_is_called_on_already_decoded_cf_compliant_dataset(
         self,
     ):
-        ds = generate_dataset(cf_compliant=True)
+        ds = generate_dataset(cf_compliant=True, has_bounds=True)
 
         with pytest.raises(KeyError):
             decode_time_units(ds)
@@ -173,11 +175,6 @@ class TestDecodeTimeUnits:
     def test_decodes_cf_compliant_time_units(self):
         # Create a dummy dataset with CF compliant time units.
         time_attrs = self.time_attrs
-        time_attrs.update({"units": "days since 2000-01-01"})
-        time_coord = xr.DataArray(
-            name="time", data=[0, 1, 2], dims=["time"], attrs=time_attrs
-        )
-        input_ds = xr.Dataset({"time": time_coord})
 
         # Create an expected dataset with properly decoded time units.
         expected_ds = xr.Dataset(
@@ -195,6 +192,13 @@ class TestDecodeTimeUnits:
             }
         )
 
+        # Update the time attrs to mimic decode_times=False
+        time_attrs.update({"units": "days since 2000-01-01"})
+        time_coord = xr.DataArray(
+            name="time", data=[0, 1, 2], dims=["time"], attrs=time_attrs
+        )
+        input_ds = xr.Dataset({"time": time_coord})
+
         # Check the resulting dataset is identical to the expected.
         result_ds = decode_time_units(input_ds)
         assert result_ds.identical(expected_ds)
diff --git a/tests/test_variable.py b/tests/test_variable.py
index 062a9d1..289705d 100644
--- a/tests/test_variable.py
+++ b/tests/test_variable.py
@@ -1,38 +1,41 @@
 import pytest
 
-from tests.fixtures import generate_dataset, lat_bnds, ts_with_bnds
+from tests.fixtures import generate_dataset, ts_with_bnds_from_parent_cf
 from xcdat.variable import open_variable
 
 
 class TestOpenVariable:
     def test_raises_error_if_variable_does_not_exist(self):
-        ds = generate_dataset(has_bounds=False)
+        ds = generate_dataset(cf_compliant=True, has_bounds=False)
 
         with pytest.raises(KeyError):
             open_variable(ds, "invalid_var")
 
     def test_raises_error_if_bounds_dim_is_missing(self):
-        ds = generate_dataset(has_bounds=False)
+        ds = generate_dataset(cf_compliant=True, has_bounds=False)
 
         with pytest.raises(KeyError):
             open_variable(ds, "ts")
 
     def test_raises_error_if_bounds_are_missing_for_coordinates(self):
-        ds = generate_dataset(has_bounds=False)
+        ds = generate_dataset(cf_compliant=True, has_bounds=True)
+
+        # Need to add CF attributes in order for cf_xarray to properly
+        # identify coordinates and bounds.
+        ds = ds.drop_vars("lat_bnds")
 
         # By adding bounds to the parent dataset, it will initiate copying bounds
         # and find that bounds are missing for the other coords (lon and time).
-        ds["lat_bnds"] = lat_bnds.copy()
         with pytest.raises(ValueError):
             open_variable(ds, "ts")
 
     def test_returns_variable_with_bounds(self):
-        ds = generate_dataset(has_bounds=True)
+        ds = generate_dataset(cf_compliant=True, has_bounds=True)
         ds.lat.attrs["bounds"] = "lat_bnds"
         ds.lon.attrs["bounds"] = "lon_bnds"
         ds.time.attrs["bounds"] = "time_bnds"
 
-        ts_expected = ts_with_bnds.copy()
+        ts_expected = ts_with_bnds_from_parent_cf.copy()
 
         ts_result = open_variable(ds, "ts")
         assert ts_result.identical(ts_expected)