Repository: xCDAT/xcdat

Issue:
[Bug]: Failure to decode non-CF compliant time
### What happened?

I tried opening a dataset with non-CF compliant time units. We don't currently have code to process this particular type of non-CF compliant time axis so I got an error. 

### What did you expect to happen?

Ideally, we could change `open_dataset` (and `open_mfdataset`) to attempt to `decode_non_cf_time` and if it fails, simply do not decode the time axis, e.g., [here](https://github.com/xCDAT/xcdat/blob/main/xcdat/dataset.py#L87):

```
        if cf_compliant_time is False:
            # XCDAT handles decoding time values with non-CF units.
            ds = xr.open_dataset(path, decode_times=False, **kwargs)
            try:
                 ds = decode_non_cf_time(ds)
            except:
                 # some log message to note that decoding failed
```

[Note that there would also need to be some error handling in `_has_cf_compliant_time`]

### Minimal Complete Verifiable Example

```python
fn = '/p/user_pub/climate_work/pochedley1/surface/best.nc'
ds = xcdat.open_dataset(fn)
```


### Relevant log output

```python
~/code/xcdat/xcdat/dataset.py in _split_time_units_attr(units_attr)
    612         raise KeyError("No 'units' attribute found for the dataset's time coordinates.")
    613 
--> 614     units, reference_date = units_attr.split(" since ")
    615 
    616     return units, reference_date

ValueError: not enough values to unpack (expected 2, got 1)
```

This error results because we were expecting units of the form "X since Y", but the units are `year A.D.` (time axis has decimal years, i.e., 1850.041667, 1850.125, ...).

### Anything else we need to know?

Should I attempt a PR?

### Environment

Recent-ish xcdat dev environment.

Discussion:
durack1: @pochedls if you want a quick check of CF validity, point your file at https://pumatest.nerc.ac.uk/cgi-bin/cf-checker.pl - catches all manner of things that aren't obvious through a quick `ncdump` scan
tomvothecoder: Oh yeah, I also rebased this branch on the latest `main` so there are no merge conflicts.
pochedls: > I updated the assertion statement since we expect the result to be the same as the original dataset.

Thank you...this makes a lot of sense and I wish this is what I had done.

> Oh yeah, I also rebased this branch on the latest `main` so there are no merge conflicts.

You just needed to do this because you just merged #262, right? Or did I make a mistake?

I will squash and merge shortly. Thank you for the review, @tomvothecoder.
tomvothecoder: > Just a note that I am getting some warnings, but I don't think I caused them:
> 
> > tests/test_regrid.py::TestRegrid2Regridder::test_output_bounds
> > tests/test_regrid.py::TestAccessor::test_convenience_methods
> > /home/pochedley1/code/xcdat/xcdat/regridder/regrid2.py:233: RuntimeWarning: invalid value encountered in true_divide
> > np.multiply(input_lon_segment, dot_weight).sum(
> > tests/test_regrid.py::TestXESMFRegridder::test_regrid
> > /home/pochedley1/bin/anaconda3/envs/xcdat_dev/lib/python3.9/site-packages/numba/np/arraymath.py:3845: DeprecationWarning: `np.MachAr` is deprecated (NumPy 1.22).
> > @overload(np.MachAr)

This is from `numba`, which was introduced in PR #44. I'll see if there's a way to fix it in a separate PR.

> > Oh yeah, I also rebased this branch on the latest `main` so there are no merge conflicts.
> 
> You just needed to do this because you just merged #262, right? Or did I make a mistake?
> 
> I will squash and merge shortly. Thank you for the review, @tomvothecoder.

No mistakes were made, #262 was merged so `main` was ahead of this branch. Just needed to rebase on top of the latest `main` to avoid conflicts.

Thanks for the great work!

Patch:
diff --git a/xcdat/dataset.py b/xcdat/dataset.py
index 320f583..d480649 100644
--- a/xcdat/dataset.py
+++ b/xcdat/dataset.py
@@ -41,9 +41,10 @@ def open_dataset(
         the Dataset, by default True. Bounds are required for many xCDAT
         features.
     decode_times: bool, optional
-        If True, decode times encoded in the standard NetCDF datetime format
-        into datetime objects. Otherwise, leave them encoded as numbers.
-        This keyword may not be supported by all the backends, by default True.
+        If True, attempt to decode times encoded in the standard NetCDF
+        datetime format into datetime objects. Otherwise, leave them encoded
+        as numbers. This keyword may not be supported by all the backends,
+        by default True.
     center_times: bool, optional
         If True, center time coordinates using the midpoint between its upper
         and lower bounds. Otherwise, use the provided time coordinates, by
@@ -84,6 +85,7 @@ def open_dataset(
         if cf_compliant_time is False:
             # XCDAT handles decoding time values with non-CF units.
             ds = xr.open_dataset(path, decode_times=False, **kwargs)
+            # attempt to decode non-cf-compliant time axis
             ds = decode_non_cf_time(ds)
         else:
             ds = xr.open_dataset(path, decode_times=True, **kwargs)
@@ -225,13 +227,14 @@ def decode_non_cf_time(dataset: xr.Dataset) -> xr.Dataset:
     numerically encoded time values (representing the offset from the reference
     date) to pandas DateOffset objects. These offset values are added to the
     reference date, forming DataArrays of datetime objects that replace the time
-    coordinate and time bounds (if they exist) values in the Dataset.
+    coordinate and time bounds (if they exist) in the Dataset.
 
     Parameters
     ----------
     dataset : xr.Dataset
         Dataset with numerically encoded time coordinates and time bounds (if
-        they exist).
+        they exist). If the time coordinates cannot be decoded then the original
+        dataset is returned.
 
     Returns
     -------
@@ -304,7 +307,14 @@ def decode_non_cf_time(dataset: xr.Dataset) -> xr.Dataset:
     time = ds.cf["T"]
     time_bounds = ds.get(time.attrs.get("bounds"), None)
     units_attr = time.attrs.get("units")
-    units, ref_date = _split_time_units_attr(units_attr)
+
+    # If the time units cannot be split into a unit and reference date, it
+    # cannot be decoded so the original dateset is returned.
+    try:
+        units, ref_date = _split_time_units_attr(units_attr)
+    except ValueError:
+        return ds
+
     ref_date = pd.to_datetime(ref_date)
 
     data = [ref_date + pd.DateOffset(**{units: offset}) for offset in time.data]
@@ -403,7 +413,13 @@ def _has_cf_compliant_time(
         return None
 
     time = ds.cf["T"]
-    units = _split_time_units_attr(time.attrs.get("units"))[0]
+
+    # If the time units attr cannot be split, it is not cf_compliant.
+    try:
+        units = _split_time_units_attr(time.attrs.get("units"))[0]
+    except ValueError:
+        return False
+
     cf_compliant = units not in NON_CF_TIME_UNITS
 
     return cf_compliant
@@ -589,6 +605,7 @@ def _preprocess_non_cf_dataset(
     if callable:
         ds_new = callable(ds)
 
+    # Attempt to decode non-cf-compliant time axis.
     ds_new = decode_non_cf_time(ds_new)
 
     return ds_new
@@ -606,11 +623,23 @@ def _split_time_units_attr(units_attr: str) -> Tuple[str, str]:
     -------
     Tuple[str, str]
         The units (e.g, "months") and the reference date (e.g., "1800-01-01").
-        If the units attribute doesn't exist for the time coordinates.
+
+    Raises
+    ------
+    KeyError
+        If the time units attribute was not found.
+
+    ValueError
+        If the time units attribute is not of the form `X since Y`.
     """
     if units_attr is None:
-        raise KeyError("No 'units' attribute found for the dataset's time coordinates.")
+        raise KeyError("The dataset's time coordinates does not have a 'units' attr.")
 
-    units, reference_date = units_attr.split(" since ")
+    if "since" in units_attr:
+        units, reference_date = units_attr.split(" since ")
+    else:
+        raise ValueError(
+            "This dataset does not have time coordinates of the form 'X since Y'."
+        )
 
     return units, reference_date


Test patch:
diff --git a/tests/fixtures.py b/tests/fixtures.py
index ca68354..d79fb05 100644
--- a/tests/fixtures.py
+++ b/tests/fixtures.py
@@ -52,7 +52,15 @@ time_non_cf = xr.DataArray(
         "standard_name": "time",
     },
 )
-
+time_non_cf_unsupported = xr.DataArray(
+    data=np.arange(1850 + 1 / 24.0, 1851 + 3 / 12.0, 1 / 12.0),
+    dims=["time"],
+    attrs={
+        "units": "year A.D.",
+        "long_name": "time",
+        "standard_name": "time",
+    },
+)
 time_bnds = xr.DataArray(
     name="time_bnds",
     data=np.array(
@@ -104,6 +112,16 @@ time_bnds_non_cf = xr.DataArray(
     dims=["time", "bnds"],
     attrs={"xcdat_bounds": "True"},
 )
+tb = []
+for t in time_non_cf_unsupported:
+    tb.append([t - 1 / 24.0, t + 1 / 24.0])
+time_bnds_non_cf_unsupported = xr.DataArray(
+    name="time_bnds",
+    data=tb,
+    coords={"time": time_non_cf_unsupported},
+    dims=["time", "bnds"],
+    attrs={"is_generated": "True"},
+)
 
 # LATITUDE
 # ========
@@ -159,7 +177,9 @@ ts_non_cf = xr.DataArray(
 )
 
 
-def generate_dataset(cf_compliant: bool, has_bounds: bool) -> xr.Dataset:
+def generate_dataset(
+    cf_compliant: bool, has_bounds: bool, unsupported: bool = False
+) -> xr.Dataset:
     """Generates a dataset using coordinate and data variable fixtures.
 
     Parameters
@@ -169,12 +189,22 @@ def generate_dataset(cf_compliant: bool, has_bounds: bool) -> xr.Dataset:
     has_bounds : bool, optional
         Include bounds for coordinates. This also adds the "bounds" attribute
         to existing coordinates to link them to their respective bounds.
+    unsupported : bool, optional
+        Create time units that are unsupported and cannot be decoded.
+        Note that cf_compliant must be set to False.
 
     Returns
     -------
     xr.Dataset
         Test dataset.
     """
+
+    if unsupported & cf_compliant:
+        raise ValueError(
+            "Cannot set cf_compliant=True and unsupported=True. \n"
+            "Set cf_compliant=False."
+        )
+
     if has_bounds:
         ds = xr.Dataset(
             data_vars={
@@ -189,8 +219,12 @@ def generate_dataset(cf_compliant: bool, has_bounds: bool) -> xr.Dataset:
             ds.coords["time"] = time_cf.copy()
             ds["time_bnds"] = time_bnds.copy()
         elif not cf_compliant:
-            ds.coords["time"] = time_non_cf.copy()
-            ds["time_bnds"] = time_bnds_non_cf.copy()
+            if unsupported:
+                ds.coords["time"] = time_non_cf_unsupported.copy()
+                ds["time_bnds"] = time_bnds_non_cf_unsupported.copy()
+            else:
+                ds.coords["time"] = time_non_cf.copy()
+                ds["time_bnds"] = time_bnds_non_cf.copy()
 
         # 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
diff --git a/tests/test_dataset.py b/tests/test_dataset.py
index ba44740..ed83809 100644
--- a/tests/test_dataset.py
+++ b/tests/test_dataset.py
@@ -37,6 +37,16 @@ class TestOpenDataset:
         expected = generate_dataset(cf_compliant=False, has_bounds=True)
         assert result.identical(expected)
 
+    def test_non_cf_compliant_and_unsupported_time_is_not_decoded(self):
+        ds = generate_dataset(cf_compliant=False, has_bounds=True, unsupported=True)
+        ds.to_netcdf(self.file_path)
+
+        # even though decode_times=True, it should fail to decode unsupported time axis
+        result = open_dataset(self.file_path, decode_times=True)
+        expected = ds
+
+        assert result.identical(expected)
+
     def test_non_cf_compliant_time_is_decoded(self):
         ds = generate_dataset(cf_compliant=False, has_bounds=False)
         ds.to_netcdf(self.file_path)