Repository: corteva/rioxarray

Issue:
open_rasterio fails on NetCDF files with non-standard time units
<!-- Please search existing issues to avoid creating duplicates. -->


#### Code Sample

NOTE: The dataset is available from the [LUH2 site](http://gsweb1vh2.umd.edu/LUH2/LUH2_v2f/MESSAGE/multiple-states_input4MIPs_landState_ScenarioMIP_UofMD-MESSAGE-ssp245-2-1-f_gn_2015-2100.nc).

```python
import rioxarray as rxr
rxr.open_rasterio('multiple-states_input4MIPs_landState_ScenarioMIP_UofMD-MESSAGE-ssp245-2-1-f_gn_2015-2100.nc')
```

#### Problem description

The [LUH2](https://luh.umd.edu) datasets have a non-standard time units (`years since 2015-01-01 0:0:0`) which causes `open_rasterio()` to fail with an exception. The error message suggest passing the flag `decode_times=False` or installing `cftime`. Neither approach solved the problem for me.

It would be great if **rioxarray** had the equivalent of `decode_times` for datasets like this one.

Caveat: I don't get an exception when opening the file remotely via http://  You will need to download the file and run it locally. 

#### Expected Output

`open_rasterio()` succeeds.

#### Environment Information

```
rioxarray (0.3.1) deps:
  rasterio: 1.2.2
    xarray: 0.17.0
      GDAL: 3.2.2

Other python deps:
     scipy: 1.6.2
    pyproj: 3.0.1
```


#### Installation method
 - pypi

#### Packages
```
Package           Version
----------------- ---------
affine            2.3.0
attrs             20.3.0
blosc             1.9.2
bokeh             2.3.1
certifi           2020.12.5
cftime            1.4.1
click             7.1.2
click-plugins     1.1.1
cligj             0.7.1
cloudpickle       1.6.0
cycler            0.10.0
dask              2021.4.0
distributed       2021.4.0
fsspec            2021.4.0
HeapDict          1.0.1
Jinja2            2.11.3
kiwisolver        1.3.1
locket            0.2.1
lz4               3.1.1
MarkupSafe        1.1.1
matplotlib        3.4.1
msgpack           1.0.0
netCDF4           1.5.6
numpy             1.18.1
packaging         20.9
pandas            1.2.4
partd             1.2.0
Pillow            8.2.0
pip               21.0.1
psutil            5.8.0
pyflakes          2.3.1
pyparsing         2.4.7
pyproj            3.0.1
python-dateutil   2.8.1
pytz              2021.1
PyYAML            5.4.1
rasterio          1.2.2
rioxarray         0.3.1
scipy             1.6.2
setuptools        56.0.0
six               1.15.0
snuggs            1.4.7
sortedcontainers  2.3.0
tblib             1.7.0
toolz             0.11.1
tornado           6.1
typing-extensions 3.7.4.3
wheel             0.36.2
xarray            0.17.0
zict              2.0.0
```

Discussion:

Patch:
diff --git a/docs/history.rst b/docs/history.rst
index dd91a886..b9b0820c 100644
--- a/docs/history.rst
+++ b/docs/history.rst
@@ -9,6 +9,7 @@ Latest
 - ENH: enable `engine="rasterio"` via xarray backend API (issue #197 pull #281)
 - ENH: Generate 2D coordinates for non-rectilinear sources (issue #290)
 - ENH: Add `encoded` kwarg to `rio.write_nodata` (discussions #313)
+- ENH: Added `decode_times` and `decode_timedelta` kwargs to `rioxarray.open_rasterio` (issue #316)
 - BUG: Use float32 for smaller dtypes when masking (discussions #302)
 - BUG: Return correct transform in `rio.transform` with non-rectilinear transform (discussions #280)
 - BUG: Update to handle WindowError in rasterio 1.2.2 (issue #286)
diff --git a/rioxarray/_io.py b/rioxarray/_io.py
index 0f72ffbd..9812279b 100644
--- a/rioxarray/_io.py
+++ b/rioxarray/_io.py
@@ -26,6 +26,7 @@ from xarray.core import indexing
 from xarray.core.dataarray import DataArray
 from xarray.core.dtypes import maybe_promote
 from xarray.core.utils import is_scalar
+from xarray.core.variable import as_variable
 
 from rioxarray.exceptions import RioXarrayError
 from rioxarray.rioxarray import _generate_spatial_coords
@@ -452,60 +453,37 @@ def _get_rasterio_attrs(riods):
     return attrs
 
 
-def _decode_datetime_cf(data_array):
+def _decode_datetime_cf(data_array, decode_times, decode_timedelta):
     """
     Decide the datetime based on CF conventions
     """
+    if decode_timedelta is None:
+        decode_timedelta = decode_times
+
     for coord in data_array.coords:
-        # stage 1: timedelta
-        if (
-            "units" in data_array[coord].attrs
-            and data_array[coord].attrs["units"] in times.TIME_UNITS
-        ):
-            units = times.pop_to(
-                data_array[coord].attrs, data_array[coord].encoding, "units"
-            )
-            new_values = times.decode_cf_timedelta(
-                data_array[coord].values, units=units
-            )
-            data_array = data_array.assign_coords(
-                {
-                    coord: IndexVariable(
-                        dims=data_array[coord].dims,
-                        data=new_values.astype(np.dtype("timedelta64[ns]")),
-                        attrs=data_array[coord].attrs,
-                        encoding=data_array[coord].encoding,
-                    )
-                }
+        time_var = None
+        if decode_times and "since" in data_array[coord].attrs.get("units", ""):
+            time_var = times.CFDatetimeCoder(use_cftime=True).decode(
+                as_variable(data_array[coord]), name=coord
             )
-
-        # stage 2: datetime
-        if (
-            "units" in data_array[coord].attrs
-            and "since" in data_array[coord].attrs["units"]
+        elif (
+            decode_timedelta
+            and data_array[coord].attrs.get("units") in times.TIME_UNITS
         ):
-            units = times.pop_to(
-                data_array[coord].attrs, data_array[coord].encoding, "units"
-            )
-            calendar = times.pop_to(
-                data_array[coord].attrs, data_array[coord].encoding, "calendar"
+            time_var = times.CFTimedeltaCoder().decode(
+                as_variable(data_array[coord]), name=coord
             )
-            dtype = times._decode_cf_datetime_dtype(
-                data_array[coord].values, units, calendar, True
-            )
-            new_values = times.decode_cf_datetime(
-                data_array[coord].values,
-                units=units,
-                calendar=calendar,
-                use_cftime=True,
+        if time_var is not None:
+            dimensions, data, attributes, encoding = variables.unpack_for_decoding(
+                time_var
             )
             data_array = data_array.assign_coords(
                 {
                     coord: IndexVariable(
-                        dims=data_array[coord].dims,
-                        data=new_values.astype(dtype),
-                        attrs=data_array[coord].attrs,
-                        encoding=data_array[coord].encoding,
+                        dims=dimensions,
+                        data=data,
+                        attrs=attributes,
+                        encoding=encoding,
                     )
                 }
             )
@@ -539,6 +517,9 @@ def _load_subdatasets(
     lock,
     masked,
     mask_and_scale,
+    decode_times,
+    decode_timedelta,
+    **open_kwargs,
 ):
     """
     Load in rasterio subdatasets
@@ -562,6 +543,9 @@ def _load_subdatasets(
             masked=masked,
             mask_and_scale=mask_and_scale,
             default_name=subdataset.split(":")[-1].lstrip("/").replace("/", "_"),
+            decode_times=decode_times,
+            decode_timedelta=decode_timedelta,
+            **open_kwargs,
         )
         if shape not in dim_groups:
             dim_groups[shape] = {rioda.name: rioda}
@@ -648,6 +632,8 @@ def open_rasterio(
     variable=None,
     group=None,
     default_name=None,
+    decode_times=True,
+    decode_timedelta=None,
     **open_kwargs,
 ):
     # pylint: disable=too-many-statements,too-many-locals,too-many-branches
@@ -706,6 +692,14 @@ def open_rasterio(
         Group name or names to use to filter loading.
     default_name: str, optional
         The name of the data array if none exists. Default is None.
+    decode_times: bool, optional
+        If True, decode times encoded in the standard NetCDF datetime format
+        into datetime objects. Otherwise, leave them encoded as numbers.
+    decode_timedelta: bool, optional
+        If True, decode variables and coordinates with time units in
+        {“days”, “hours”, “minutes”, “seconds”, “milliseconds”, “microseconds”}
+        into timedelta objects. If False, leave them encoded as numbers.
+        If None (default), assume the same value of decode_time.
     **open_kwargs: kwargs, optional
         Optional keyword arguments to pass into rasterio.open().
 
@@ -777,6 +771,9 @@ def open_rasterio(
             lock=lock,
             masked=masked,
             mask_and_scale=mask_and_scale,
+            decode_times=decode_times,
+            decode_timedelta=decode_timedelta,
+            **open_kwargs,
         )
 
     if vrt_params is not None:
@@ -840,7 +837,9 @@ def open_rasterio(
 
     # update attributes from NetCDF attributess
     _load_netcdf_attrs(riods.tags(), result)
-    result = _decode_datetime_cf(result)
+    result = _decode_datetime_cf(
+        result, decode_times=decode_times, decode_timedelta=decode_timedelta
+    )
 
     # make sure the _FillValue is correct dtype
     if "_FillValue" in attrs:
diff --git a/rioxarray/xarray_plugin.py b/rioxarray/xarray_plugin.py
index 449fbc90..87c27eaa 100644
--- a/rioxarray/xarray_plugin.py
+++ b/rioxarray/xarray_plugin.py
@@ -41,6 +41,8 @@ class RasterioBackend(xr.backends.common.BackendEntrypoint):
         variable=None,
         group=None,
         default_name="band_data",
+        decode_times=True,
+        decode_timedelta=None,
         open_kwargs=None,
     ):
         if open_kwargs is None:
@@ -56,6 +58,8 @@ class RasterioBackend(xr.backends.common.BackendEntrypoint):
             variable=variable,
             group=group,
             default_name=default_name,
+            decode_times=decode_times,
+            decode_timedelta=decode_timedelta,
             **open_kwargs,
         )
         if isinstance(ds, xr.DataArray):


Test patch:
diff --git a/test/integration/test_integration__io.py b/test/integration/test_integration__io.py
index dbe7c214..19edf3ea 100644
--- a/test/integration/test_integration__io.py
+++ b/test/integration/test_integration__io.py
@@ -981,6 +981,27 @@ def test_nc_attr_loading(open_rasterio):
         assert str(rds.time.values[1]) == "2016-12-29 12:52:42.347451"
 
 
+@pytest.mark.xfail(
+    LooseVersion(rasterio.__gdal_version__) < LooseVersion("3.0.4"),
+    reason="This was fixed in GDAL 3.0.4",
+)
+def test_nc_attr_loading__disable_decode_times(open_rasterio):
+    with open_rasterio(
+        os.path.join(TEST_INPUT_DATA_DIR, "PLANET_SCOPE_3D.nc"), decode_times=False
+    ) as rds:
+        assert rds.dims == {"y": 10, "x": 10, "time": 2}
+        assert rds.attrs == {"coordinates": "spatial_ref"}
+        assert rds.y.attrs["units"] == "metre"
+        assert rds.x.attrs["units"] == "metre"
+        assert rds.time.encoding == {}
+        assert np.isnan(rds.time.attrs.pop("_FillValue"))
+        assert rds.time.attrs == {
+            "units": "seconds since 2016-12-19T10:27:29.687763",
+            "calendar": "proleptic_gregorian",
+        }
+        assert_array_equal(rds.time.values, [0, 872712.659688])
+
+
 def test_lockless():
     with rioxarray.open_rasterio(
         os.path.join(TEST_INPUT_DATA_DIR, "PLANET_SCOPE_3D.nc"), lock=False, chunks=True