Repository: xCDAT/xcdat

Issue:
[Bug]: `get_bounds()` can't get bounds for non-CF axes without `var_key` argument
### What happened?

Notice in the `get_bounds()` conditional below, we retrieve bounds for the Dataset object using `_get_bounds_keys()` (which does not work with non-CF axes). On the other hand, if we pass `var_key`, we get the bounds using the "bounds" attribute of non-CF axes (which works fine).
https://github.com/xCDAT/xcdat/blob/e3dda64feab046046b6a3ee8047faec1ab01f8ec/xcdat/bounds.py#L233-L250 

https://github.com/xCDAT/xcdat/blob/e3dda64feab046046b6a3ee8047faec1ab01f8ec/xcdat/bounds.py#L477-L508


### What did you expect to happen? Are there are possible answers you came across?

We should be able to get bounds for non-CF compliant axes that xCDAT can map to via dimension name (e.g., "lev", "lat", "lon"). We need to update `_get_bounds_keys()` to consider the "bounds" attribute.

### Minimal Complete Verifiable Example (MVCE)

```python
# %%
import xcdat as xc
import xarray as xr

# Create a dummy dataset
data = xr.DataArray([[1, 2, 3], [4, 5, 6]], dims=("lat", "lon"))
ds = xr.Dataset({"data": data})

# Add a "lev" Z axis with bounds
lev = xr.DataArray([1000, 2000], dims=("lev"))
ds["lev"] = lev
ds = ds.bounds.add_missing_bounds(axes="Z")

# %%
# Notice how we can map to "lev" without CF attributes because "lev" is a standard
# Z axis name.
xc.get_dim_coords(ds, "Z")

# %%
# However, we can't map to "lev_bnds".
ds.bounds.get_bounds(axis="Z")
# KeyError: "No bounds data variables were found for the 'Z' axis. Make sure the dataset has bound data vars and their names match the 'bounds' attributes found on their related time coordinate variables. Alternatively, you can add bounds with `ds.bounds.add_missing_bounds()` or `ds.bounds.add_bounds()`."

# %%
# Workaround 1 - add CF attribute(s)
ds_cf = ds.copy()
ds_cf["lev"].attrs["axis"] = "Z"
ds_cf.bounds.get_bounds(axis="Z")

# %%
# Workaround 2 - pass the var_key arg to get bounds via "bounds" key
ds.bounds.get_bounds(axis="Z", var_key="lev")
```


### Relevant log output

_No response_

### Anything else we need to know?

Found this bug through e3sm_diags GH issue [863](https://github.com/E3SM-Project/e3sm_diags/issues/863)

### Environment

Latest xCDAT and `main`

Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/xCDAT/xcdat/pull/708?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Project coverage is 100.00%. Comparing base [(`e3dda64`)](https://app.codecov.io/gh/xCDAT/xcdat/commit/e3dda64feab046046b6a3ee8047faec1ab01f8ec?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT) to head [(`80795b2`)](https://app.codecov.io/gh/xCDAT/xcdat/commit/80795b2c8c1f9e25759137b0a7d7c3d495c6d40f?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT).

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff            @@
##              main      #708   +/-   ##
=========================================
  Coverage   100.00%   100.00%           
=========================================
  Files           15        15           
  Lines         1546      1555    +9     
=========================================
+ Hits          1546      1555    +9     
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/xCDAT/xcdat/pull/708?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT).


Patch:
diff --git a/Makefile b/Makefile
index 9135ed7..68e1e7d 100644
--- a/Makefile
+++ b/Makefile
@@ -83,4 +83,4 @@ docs: ## generate Sphinx HTML documentation, including API docs
 # Build
 # ----------------------
 install: clean ## install the package to the active Python's site-packages
-	python setup.py install
+	python -m pip install .
diff --git a/docs/api.rst b/docs/api.rst
index 1ced17a..1e10764 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -27,6 +27,7 @@ Below is a list of top-level API functions that are available in ``xcdat``.
     compare_datasets
     get_dim_coords
     get_dim_keys
+    create_bounds
     create_axis
     create_gaussian_grid
     create_global_mean_grid
@@ -83,13 +84,13 @@ Classes
 .. autosummary::
    :toctree: generated/
 
-    xcdat.bounds.BoundsAccessor
-    xcdat.spatial.SpatialAccessor
-    xcdat.temporal.TemporalAccessor
-    xcdat.regridder.accessor.RegridderAccessor
-    xcdat.regridder.regrid2.Regrid2Regridder
-    xcdat.regridder.xesmf.XESMFRegridder
-    xcdat.regridder.xgcm.XGCMRegridder
+    bounds.BoundsAccessor
+    spatial.SpatialAccessor
+    temporal.TemporalAccessor
+    regridder.accessor.RegridderAccessor
+    regridder.regrid2.Regrid2Regridder
+    regridder.xesmf.XESMFRegridder
+    regridder.xgcm.XGCMRegridder
 
 .. currentmodule:: xarray
 
diff --git a/xcdat/__init__.py b/xcdat/__init__.py
index cae9909..3881359 100644
--- a/xcdat/__init__.py
+++ b/xcdat/__init__.py
@@ -6,7 +6,7 @@ from xcdat.axis import (  # noqa: F401
     get_dim_keys,
     swap_lon_axis,
 )
-from xcdat.bounds import BoundsAccessor  # noqa: F401
+from xcdat.bounds import BoundsAccessor, create_bounds  # noqa: F401
 from xcdat.dataset import decode_time, open_dataset, open_mfdataset  # noqa: F401
 from xcdat.regridder.accessor import RegridderAccessor  # noqa: F401
 from xcdat.regridder.grid import (  # noqa: F401
diff --git a/xcdat/bounds.py b/xcdat/bounds.py
index dfbc0af..4ee799c 100644
--- a/xcdat/bounds.py
+++ b/xcdat/bounds.py
@@ -236,18 +236,7 @@ class BoundsAccessor:
         else:
             # Get the obj in the Dataset using the key.
             obj = _get_data_var(self._dataset, key=var_key)
-
-            # Check if the object is a data variable or a coordinate variable.
-            # If it is a data variable, derive the axis coordinate variable.
-            if obj.name in list(self._dataset.data_vars):
-                coord = get_dim_coords(obj, axis)
-            elif obj.name in list(self._dataset.coords):
-                coord = obj
-
-            try:
-                bounds_keys = [coord.attrs["bounds"]]
-            except KeyError:
-                bounds_keys = []
+            bounds_keys = self._get_bounds_from_attr(obj, axis)
 
         if len(bounds_keys) == 0:
             raise KeyError(
@@ -505,8 +494,53 @@ class BoundsAccessor:
         except KeyError:
             pass
 
+        keys_from_attr = self._get_bounds_from_attr(self._dataset, axis)
+        keys = keys + keys_from_attr
+
         return list(set(keys))
 
+    def _get_bounds_from_attr(
+        self, obj: xr.DataArray | xr.Dataset, axis: CFAxisKey
+    ) -> List[str]:
+        """Retrieve bounds attribute keys from the given xarray object.
+
+        This method extracts the "bounds" attribute keys from the coordinates
+        of the specified axis in the provided xarray DataArray or Dataset.
+
+        Parameters:
+        -----------
+        obj : xr.DataArray | xr.Dataset
+            The xarray object from which to retrieve the bounds attribute keys.
+        axis : CFAxisKey
+            The CF axis key ("X", "Y", "T", or "Z").
+
+        Returns:
+        --------
+        List[str]
+            A list of bounds attribute keys found in the coordinates of the
+            specified axis. Otherwise, an empty list is returned.
+        """
+        coords_obj = get_dim_coords(obj, axis)
+        bounds_keys: List[str] = []
+
+        if isinstance(coords_obj, xr.DataArray):
+            bounds_keys = self._extract_bounds_key(coords_obj, bounds_keys)
+        elif isinstance(coords_obj, xr.Dataset):
+            for coord in coords_obj.coords.values():
+                bounds_keys = self._extract_bounds_key(coord, bounds_keys)
+
+        return bounds_keys
+
+    def _extract_bounds_key(
+        self, coords_obj: xr.DataArray, bounds_keys: List[str]
+    ) -> List[str]:
+        bnds_key = coords_obj.attrs.get("bounds")
+
+        if bnds_key is not None:
+            bounds_keys.append(bnds_key)
+
+        return bounds_keys
+
     def _create_time_bounds(  # noqa: C901
         self,
         time: xr.DataArray,


Test patch:
diff --git a/tests/test_bounds.py b/tests/test_bounds.py
index 2c519fc..00838e3 100644
--- a/tests/test_bounds.py
+++ b/tests/test_bounds.py
@@ -287,6 +287,27 @@ class TestGetBounds:
 
         assert result.identical(expected)
 
+    def test_returns_single_dataset_axis_bounds_as_a_dataarray_object_for_non_cf_axis(
+        self,
+    ):
+        ds = xr.Dataset(
+            coords={
+                "lat": xr.DataArray(
+                    data=np.ones(3),
+                    dims="lat",
+                    attrs={"bounds": "lat_bnds"},
+                )
+            },
+            data_vars={
+                "lat_bnds": xr.DataArray(data=np.ones((3, 3)), dims=["lat", "bnds"])
+            },
+        )
+
+        result = ds.bounds.get_bounds("Y")
+        expected = ds.lat_bnds
+
+        assert result.identical(expected)
+
     def test_returns_multiple_dataset_axis_bounds_as_a_dataset_object(self):
         ds = xr.Dataset(
             coords={
@@ -321,6 +342,40 @@ class TestGetBounds:
 
         assert result.identical(expected)
 
+    def test_returns_multiple_dataset_axis_bounds_as_a_dataset_object_for_non_cf_axis(
+        self,
+    ):
+        ds = xr.Dataset(
+            coords={
+                "lat": xr.DataArray(
+                    data=np.ones(3),
+                    dims="lat",
+                    attrs={
+                        "bounds": "lat_bnds",
+                    },
+                ),
+                "latitude": xr.DataArray(
+                    data=np.ones(3),
+                    dims="latitude",
+                    attrs={
+                        "bounds": "latitude_bnds",
+                    },
+                ),
+            },
+            data_vars={
+                "var": xr.DataArray(data=np.ones(3), dims=["lat"]),
+                "lat_bnds": xr.DataArray(data=np.ones((3, 3)), dims=["lat", "bnds"]),
+                "latitude_bnds": xr.DataArray(
+                    data=np.ones((3, 3)), dims=["latitude", "bnds"]
+                ),
+            },
+        )
+
+        result = ds.bounds.get_bounds("Y")
+        expected = ds.drop_vars("var")
+
+        assert result.identical(expected)
+
 
 class TestAddBounds:
     @pytest.fixture(autouse=True)