Repository: pydata/xarray

Issue:
⚠️ Nightly upstream-dev CI failed ⚠️
[Workflow Run URL](https://github.com/pydata/xarray/actions/runs/13021957757)
<details><summary>Python 3.12 Test Summary</summary>

```
xarray/tests/test_calendar_ops.py::test_convert_calendar_360_days[date-standard-360_day-D]: AssertionError: 
Arrays are not equal

Mismatched elements: 1 / 12 (8.33%)
Max absolute difference among violations: 1
Max relative difference among violations: 0.03448276
 ACTUAL: array(30)
 DESIRED: array([30, 29, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30])
xarray/tests/test_calendar_ops.py::test_convert_calendar_360_days[date-360_day-proleptic_gregorian-D]: AssertionError: 
Arrays are not equal

Mismatched elements: 1 / 12 (8.33%)
Max absolute difference among violations: 1
Max relative difference among violations: 0.03448276
 ACTUAL: array(30)
 DESIRED: array([30, 29, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30])
xarray/tests/test_calendar_ops.py::test_convert_calendar_360_days[date-proleptic_gregorian-360_day-4h]: AssertionError: 
Arrays are not equal

Mismatched elements: 1 / 12 (8.33%)
Max absolute difference among violations: 1
Max relative difference among violations: 0.03448276
 ACTUAL: array(30)
 DESIRED: array([30, 29, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30])
xarray/tests/test_calendar_ops.py::test_convert_calendar_360_days[year-standard-360_day-D]: AssertionError: 
Arrays are not equal

Mismatched elements: 1 / 12 (8.33%)
Max absolute difference among violations: 1
Max relative difference among violations: 0.03448276
 ACTUAL: array(30)
 DESIRED: array([30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 29])
xarray/tests/test_calendar_ops.py::test_convert_calendar_360_days[year-360_day-proleptic_gregorian-D]: AssertionError: 
Arrays are not equal

Mismatched elements: 6 / 12 (50%)
Max absolute difference among violations: 1
Max relative difference among violations: 0.03448276
 ACTUAL: array(30)
 DESIRED: array([30, 29, 30, 30, 31, 30, 30, 31, 30, 31, 29, 31])
xarray/tests/test_calendar_ops.py::test_convert_calendar_360_days[year-proleptic_gregorian-360_day-4h]: AssertionError: 
Arrays are not equal

Mismatched elements: 1 / 12 (8.33%)
Max absolute difference among violations: 1
Max relative difference among violations: 0.03448276
 ACTUAL: array(30)
 DESIRED: array([30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 29])
xarray/tests/test_dask.py::TestDataArrayAndDataset::test_groupby_first_last[first]: AssertionError: assert False
 +  where False = isinstance(array('a', dtype='<U1'), <class 'dask.array.core.Array'>)
 +    where array('a', dtype='<U1') = <xarray.DataArray 'ab' ()> Size: 4B\narray('a', dtype='<U1')\nCoordinates:\n    ab       <U1 4B 'a'.data
 +    and   <class 'dask.array.core.Array'> = da.Array
xarray/tests/test_dask.py::TestDataArrayAndDataset::test_groupby_first_last[last]: AssertionError: assert False
 +  where False = isinstance(array('a', dtype='<U1'), <class 'dask.array.core.Array'>)
 +    where array('a', dtype='<U1') = <xarray.DataArray 'ab' ()> Size: 4B\narray('a', dtype='<U1')\nCoordinates:\n    ab       <U1 4B 'a'.data
 +    and   <class 'dask.array.core.Array'> = da.Array
xarray/tests/test_groupby.py::TestDataArrayGroupBy::test_stack_groupby_unsorted_coord: AssertionError: Left and right DataArray objects are not equal
Differing dimensions:
    () != (z: 4)
Differing values:
L
    array(0)
R
    array([0, 1, 2, 3])
Differing coordinates:
L   z        object 8B (0, 3)
R * z        (z) object 32B MultiIndex
Coordinates only on the right object:
  * x        (z) int64 32B 0 0 1 1
  * y        (z) int64 32B 3 2 3 2
xarray/tests/test_groupby.py::TestDataArrayGroupBy::test_groupby_first_and_last: AssertionError: Left and right DataArray objects are not identical
Differing dimensions:
    (ab: 2) != ()
Differing values:
L
    array([1, 3])
R
    array(1)
Differing coordinates:
L * ab       (ab) <U1 8B 'a' 'b'
R   ab       <U1 4B 'a'
xarray/tests/test_groupby.py::TestDataArrayResample::test_resample_first: AssertionError: Left and right DataArray objects are not identical
Differing dimensions:
    (time: 10) != ()
Differing values:
L
    array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
R
    array(0)
Differing coordinates:
L * time     (time) datetime64[ns] 80B 2000-01-01 ... 2000-01-03T06:00:00
R   time     datetime64[ns] 8B 2000-01-01
xarray/tests/test_groupby.py::TestDatasetResample::test_resample_and_first: AssertionError
xarray/tests/test_groupby.py::TestDatasetResample::test_resample_by_last_discarding_attrs: TypeError: string indices must be integers, not 'str'
```

</details>


Discussion:

Patch:
diff --git a/doc/internals/time-coding.rst b/doc/internals/time-coding.rst
index a7e0d5de..442b749a 100644
--- a/doc/internals/time-coding.rst
+++ b/doc/internals/time-coding.rst
@@ -473,3 +473,65 @@ on-disk resolution, if possible.
 
     coder = xr.coders.CFDatetimeCoder(time_unit="s")
     xr.open_dataset("test-datetimes2.nc", decode_times=coder)
+
+Similar logic applies for decoding timedelta values. The default resolution is
+``"ns"``:
+
+.. ipython:: python
+
+    attrs = {"units": "hours"}
+    ds = xr.Dataset({"time": ("time", [0, 1, 2, 3], attrs)})
+    ds.to_netcdf("test-timedeltas1.nc")
+
+.. ipython:: python
+    :okwarning:
+
+    xr.open_dataset("test-timedeltas1.nc")
+
+By default, timedeltas will be decoded to the same resolution as datetimes:
+
+.. ipython:: python
+    :okwarning:
+
+    coder = xr.coders.CFDatetimeCoder(time_unit="s")
+    xr.open_dataset("test-timedeltas1.nc", decode_times=coder)
+
+but if one would like to decode timedeltas to a different resolution, one can
+provide a coder specifically for timedeltas to ``decode_timedelta``:
+
+.. ipython:: python
+
+    timedelta_coder = xr.coders.CFTimedeltaCoder(time_unit="ms")
+    xr.open_dataset(
+        "test-timedeltas1.nc", decode_times=coder, decode_timedelta=timedelta_coder
+    )
+
+As with datetimes, if a coarser unit is requested the timedeltas are decoded
+into their native on-disk resolution, if possible:
+
+.. ipython:: python
+
+    attrs = {"units": "milliseconds"}
+    ds = xr.Dataset({"time": ("time", [0, 1, 2, 3], attrs)})
+    ds.to_netcdf("test-timedeltas2.nc")
+
+.. ipython:: python
+    :okwarning:
+
+    xr.open_dataset("test-timedeltas2.nc")
+
+.. ipython:: python
+    :okwarning:
+
+    coder = xr.coders.CFDatetimeCoder(time_unit="s")
+    xr.open_dataset("test-timedeltas2.nc", decode_times=coder)
+
+To opt-out of timedelta decoding (see issue `Undesired decoding to timedelta64 <https://github.com/pydata/xarray/issues/1621>`_) pass ``False`` to ``decode_timedelta``:
+
+.. ipython:: python
+
+    xr.open_dataset("test-timedeltas2.nc", decode_timedelta=False)
+
+.. note::
+    Note that in the future the default value of ``decode_timedelta`` will be
+    ``False`` rather than ``None``.
diff --git a/doc/whats-new.rst b/doc/whats-new.rst
index 75ec2273..70b5bf2c 100644
--- a/doc/whats-new.rst
+++ b/doc/whats-new.rst
@@ -19,23 +19,36 @@ What's New
 v2025.01.2 (unreleased)
 -----------------------
 
-This release brings non-nanosecond datetime resolution to xarray. In the
-last couple of releases xarray has been prepared for that change. The code had
-to be changed and adapted in numerous places, affecting especially the test suite.
-The documentation has been updated accordingly and a new internal chapter
-on :ref:`internals.timecoding` has been added.
-
-To make the transition as smooth as possible this is designed to be fully backwards
-compatible, keeping the current default of ``'ns'`` resolution on decoding.
-To opt-in decoding into other resolutions (``'us'``, ``'ms'`` or ``'s'``) the
-new :py:class:`coders.CFDatetimeCoder` is used as parameter to ``decode_times``
-kwarg (see also :ref:`internals.default_timeunit`):
+This release brings non-nanosecond datetime and timedelta resolution to xarray.
+In the last couple of releases xarray has been prepared for that change. The
+code had to be changed and adapted in numerous places, affecting especially the
+test suite. The documentation has been updated accordingly and a new internal
+chapter on :ref:`internals.timecoding` has been added.
+
+To make the transition as smooth as possible this is designed to be fully
+backwards compatible, keeping the current default of ``'ns'`` resolution on
+decoding. To opt-into decoding to other resolutions (``'us'``, ``'ms'`` or
+``'s'``) an instance of the newly public :py:class:`coders.CFDatetimeCoder`
+class can be passed through the ``decode_times`` keyword argument (see also
+:ref:`internals.default_timeunit`):
 
 .. code-block:: python
 
     coder = xr.coders.CFDatetimeCoder(time_unit="s")
     ds = xr.open_dataset(filename, decode_times=coder)
 
+Similar control of the resoution of decoded timedeltas can be achieved through
+passing a :py:class:`coders.CFTimedeltaCoder` instance to the
+``decode_timedelta`` keyword argument:
+
+.. code-block:: python
+
+    coder = xr.coders.CFTimedeltaCoder(time_unit="s")
+    ds = xr.open_dataset(filename, decode_timedelta=coder)
+
+though by default timedeltas will be decoded to the same ``time_unit`` as
+datetimes.
+
 There might slight changes when encoding/decoding times as some warning and
 error messages have been removed or rewritten. Xarray will now also allow
 non-nanosecond datetimes (with ``'us'``, ``'ms'`` or ``'s'`` resolution) when
@@ -50,7 +63,7 @@ eventually be deprecated.
 
 New Features
 ~~~~~~~~~~~~
-- Relax nanosecond datetime restriction in CF time decoding (:issue:`7493`, :pull:`9618`, :pull:`9977`).
+- Relax nanosecond datetime / timedelta restriction in CF time decoding (:issue:`7493`, :pull:`9618`, :pull:`9966`, :pull:`9977`).
   By `Kai Mühlbauer <https://github.com/kmuehlbauer>`_ and `Spencer Clark <https://github.com/spencerkclark>`_.
 - Enable the ``compute=False`` option in :py:meth:`DataTree.to_zarr`. (:pull:`9958`).
   By `Sam Levang <https://github.com/slevang>`_.
@@ -72,6 +85,12 @@ Breaking changes
 
 Deprecations
 ~~~~~~~~~~~~
+- In a future version of xarray decoding of variables into
+  :py:class:`numpy.timedelta64` values will be disabled by default. To silence
+  warnings associated with this, set ``decode_timedelta`` to ``True``,
+  ``False``, or a :py:class:`coders.CFTimedeltaCoder` instance when opening
+  data (:issue:`1621`, :pull:`9966`). By `Spencer Clark
+  <https://github.com/spencerkclark>`_.
 
 
 Bug fixes
diff --git a/xarray/backends/api.py b/xarray/backends/api.py
index ef97bfb8..db76ac98 100644
--- a/xarray/backends/api.py
+++ b/xarray/backends/api.py
@@ -33,7 +33,7 @@ from xarray.backends.common import (
     _normalize_path,
 )
 from xarray.backends.locks import _get_scheduler
-from xarray.coders import CFDatetimeCoder
+from xarray.coders import CFDatetimeCoder, CFTimedeltaCoder
 from xarray.core import indexing
 from xarray.core.combine import (
     _infer_concat_order_from_positions,
@@ -487,7 +487,10 @@ def open_dataset(
     | CFDatetimeCoder
     | Mapping[str, bool | CFDatetimeCoder]
     | None = None,
-    decode_timedelta: bool | Mapping[str, bool] | None = None,
+    decode_timedelta: bool
+    | CFTimedeltaCoder
+    | Mapping[str, bool | CFTimedeltaCoder]
+    | None = None,
     use_cftime: bool | Mapping[str, bool] | None = None,
     concat_characters: bool | Mapping[str, bool] | None = None,
     decode_coords: Literal["coordinates", "all"] | bool | None = None,
@@ -555,11 +558,14 @@ def open_dataset(
         Pass a mapping, e.g. ``{"my_variable": False}``,
         to toggle this feature per-variable individually.
         This keyword may not be supported by all the backends.
-    decode_timedelta : bool or dict-like, optional
+    decode_timedelta : bool, CFTimedeltaCoder, or dict-like, 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.
+        If None (default), assume the same value of ``decode_times``; if
+        ``decode_times`` is a :py:class:`coders.CFDatetimeCoder` instance, this
+        takes the form of a :py:class:`coders.CFTimedeltaCoder` instance with a
+        matching ``time_unit``.
         Pass a mapping, e.g. ``{"my_variable": False}``,
         to toggle this feature per-variable individually.
         This keyword may not be supported by all the backends.
@@ -712,7 +718,7 @@ def open_dataarray(
     | CFDatetimeCoder
     | Mapping[str, bool | CFDatetimeCoder]
     | None = None,
-    decode_timedelta: bool | None = None,
+    decode_timedelta: bool | CFTimedeltaCoder | None = None,
     use_cftime: bool | None = None,
     concat_characters: bool | None = None,
     decode_coords: Literal["coordinates", "all"] | bool | None = None,
@@ -785,7 +791,10 @@ def open_dataarray(
         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.
+        If None (default), assume the same value of ``decode_times``; if
+        ``decode_times`` is a :py:class:`coders.CFDatetimeCoder` instance, this
+        takes the form of a :py:class:`coders.CFTimedeltaCoder` instance with a
+        matching ``time_unit``.
         This keyword may not be supported by all the backends.
     use_cftime: bool, optional
         Only relevant if encoded dates come from a standard calendar
@@ -927,7 +936,10 @@ def open_datatree(
     | CFDatetimeCoder
     | Mapping[str, bool | CFDatetimeCoder]
     | None = None,
-    decode_timedelta: bool | Mapping[str, bool] | None = None,
+    decode_timedelta: bool
+    | CFTimedeltaCoder
+    | Mapping[str, bool | CFTimedeltaCoder]
+    | None = None,
     use_cftime: bool | Mapping[str, bool] | None = None,
     concat_characters: bool | Mapping[str, bool] | None = None,
     decode_coords: Literal["coordinates", "all"] | bool | None = None,
@@ -995,7 +1007,10 @@ def open_datatree(
         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.
+        If None (default), assume the same value of ``decode_times``; if
+        ``decode_times`` is a :py:class:`coders.CFDatetimeCoder` instance, this
+        takes the form of a :py:class:`coders.CFTimedeltaCoder` instance with a
+        matching ``time_unit``.
         Pass a mapping, e.g. ``{"my_variable": False}``,
         to toggle this feature per-variable individually.
         This keyword may not be supported by all the backends.
@@ -1150,7 +1165,10 @@ def open_groups(
     | CFDatetimeCoder
     | Mapping[str, bool | CFDatetimeCoder]
     | None = None,
-    decode_timedelta: bool | Mapping[str, bool] | None = None,
+    decode_timedelta: bool
+    | CFTimedeltaCoder
+    | Mapping[str, bool | CFTimedeltaCoder]
+    | None = None,
     use_cftime: bool | Mapping[str, bool] | None = None,
     concat_characters: bool | Mapping[str, bool] | None = None,
     decode_coords: Literal["coordinates", "all"] | bool | None = None,
@@ -1222,9 +1240,10 @@ def open_groups(
         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.
-        Pass a mapping, e.g. ``{"my_variable": False}``,
-        to toggle this feature per-variable individually.
+        If None (default), assume the same value of ``decode_times``; if
+        ``decode_times`` is a :py:class:`coders.CFDatetimeCoder` instance, this
+        takes the form of a :py:class:`coders.CFTimedeltaCoder` instance with a
+        matching ``time_unit``.
         This keyword may not be supported by all the backends.
     use_cftime: bool or dict-like, optional
         Only relevant if encoded dates come from a standard calendar
diff --git a/xarray/coders.py b/xarray/coders.py
index 238ac714..4f0a32dc 100644
--- a/xarray/coders.py
+++ b/xarray/coders.py
@@ -3,8 +3,6 @@ This module provides coder objects that encapsulate the
 "encoding/decoding" process.
 """
 
-from xarray.coding.times import CFDatetimeCoder
+from xarray.coding.times import CFDatetimeCoder, CFTimedeltaCoder
 
-__all__ = [
-    "CFDatetimeCoder",
-]
+__all__ = ["CFDatetimeCoder", "CFTimedeltaCoder"]
diff --git a/xarray/coding/times.py b/xarray/coding/times.py
index fd99a55a..ad5e8653 100644
--- a/xarray/coding/times.py
+++ b/xarray/coding/times.py
@@ -1343,6 +1343,21 @@ class CFDatetimeCoder(VariableCoder):
 
 
 class CFTimedeltaCoder(VariableCoder):
+    """Coder for CF Timedelta coding.
+
+    Parameters
+    ----------
+    time_unit : PDDatetimeUnitOptions
+          Target resolution when decoding timedeltas. Defaults to "ns".
+    """
+
+    def __init__(
+        self,
+        time_unit: PDDatetimeUnitOptions = "ns",
+    ) -> None:
+        self.time_unit = time_unit
+        self._emit_decode_timedelta_future_warning = False
+
     def encode(self, variable: Variable, name: T_Name = None) -> Variable:
         if np.issubdtype(variable.data.dtype, np.timedelta64):
             dims, data, attrs, encoding = unpack_for_encoding(variable)
@@ -1359,12 +1374,21 @@ class CFTimedeltaCoder(VariableCoder):
     def decode(self, variable: Variable, name: T_Name = None) -> Variable:
         units = variable.attrs.get("units", None)
         if isinstance(units, str) and units in TIME_UNITS:
+            if self._emit_decode_timedelta_future_warning:
+                emit_user_level_warning(
+                    "In a future version of xarray decode_timedelta will "
+                    "default to False rather than None. To silence this "
+                    "warning, set decode_timedelta to True, False, or a "
+                    "'CFTimedeltaCoder' instance.",
+                    FutureWarning,
+                )
             dims, data, attrs, encoding = unpack_for_decoding(variable)
 
             units = pop_to(attrs, encoding, "units")
-            transform = partial(decode_cf_timedelta, units=units)
-            # todo: check, if we can relax this one here, too
-            dtype = np.dtype("timedelta64[ns]")
+            dtype = np.dtype(f"timedelta64[{self.time_unit}]")
+            transform = partial(
+                decode_cf_timedelta, units=units, time_unit=self.time_unit
+            )
             data = lazy_elemwise_func(data, transform, dtype=dtype)
 
             return Variable(dims, data, attrs, encoding, fastpath=True)
diff --git a/xarray/conventions.py b/xarray/conventions.py
index 485c9ac0..f67af95b 100644
--- a/xarray/conventions.py
+++ b/xarray/conventions.py
@@ -1,14 +1,15 @@
 from __future__ import annotations
 
 import itertools
+import warnings
 from collections import defaultdict
 from collections.abc import Hashable, Iterable, Mapping, MutableMapping
 from typing import TYPE_CHECKING, Any, Literal, TypeVar, Union
 
 import numpy as np
 
-from xarray.coders import CFDatetimeCoder
-from xarray.coding import strings, times, variables
+from xarray.coders import CFDatetimeCoder, CFTimedeltaCoder
+from xarray.coding import strings, variables
 from xarray.coding.variables import SerializationWarning, pop_to
 from xarray.core import indexing
 from xarray.core.common import (
@@ -90,7 +91,7 @@ def encode_cf_variable(
 
     for coder in [
         CFDatetimeCoder(),
-        times.CFTimedeltaCoder(),
+        CFTimedeltaCoder(),
         variables.CFScaleOffsetCoder(),
         variables.CFMaskCoder(),
         variables.NativeEnumCoder(),
@@ -114,7 +115,7 @@ def decode_cf_variable(
     decode_endianness: bool = True,
     stack_char_dim: bool = True,
     use_cftime: bool | None = None,
-    decode_timedelta: bool | None = None,
+    decode_timedelta: bool | CFTimedeltaCoder | None = None,
 ) -> Variable:
     """
     Decodes a variable which may hold CF encoded information.
@@ -158,6 +159,8 @@ def decode_cf_variable(
 
         .. deprecated:: 2025.01.1
            Please pass a :py:class:`coders.CFDatetimeCoder` instance initialized with ``use_cftime`` to the ``decode_times`` kwarg instead.
+    decode_timedelta : None, bool, or CFTimedeltaCoder
+        Decode cf timedeltas ("hours") to np.timedelta64.
 
     Returns
     -------
@@ -170,8 +173,12 @@ def decode_cf_variable(
 
     original_dtype = var.dtype
 
+    decode_timedelta_was_none = decode_timedelta is None
     if decode_timedelta is None:
-        decode_timedelta = True if decode_times else False
+        if isinstance(decode_times, CFDatetimeCoder):
+            decode_timedelta = CFTimedeltaCoder(time_unit=decode_times.time_unit)
+        else:
+            decode_timedelta = True if decode_times else False
 
     if concat_characters:
         if stack_char_dim:
@@ -193,7 +200,12 @@ def decode_cf_variable(
             var = coder.decode(var, name=name)
 
     if decode_timedelta:
-        var = times.CFTimedeltaCoder().decode(var, name=name)
+        if not isinstance(decode_timedelta, CFTimedeltaCoder):
+            decode_timedelta = CFTimedeltaCoder()
+        decode_timedelta._emit_decode_timedelta_future_warning = (
+            decode_timedelta_was_none
+        )
+        var = decode_timedelta.decode(var, name=name)
     if decode_times:
         # remove checks after end of deprecation cycle
         if not isinstance(decode_times, CFDatetimeCoder):
@@ -335,13 +347,20 @@ def decode_cf_variables(
     decode_coords: bool | Literal["coordinates", "all"] = True,
     drop_variables: T_DropVariables = None,
     use_cftime: bool | Mapping[str, bool] | None = None,
-    decode_timedelta: bool | Mapping[str, bool] | None = None,
+    decode_timedelta: bool
+    | CFTimedeltaCoder
+    | Mapping[str, bool | CFTimedeltaCoder]
+    | None = None,
 ) -> tuple[T_Variables, T_Attrs, set[Hashable]]:
     """
     Decode several CF encoded variables.
 
     See: decode_cf_variable
     """
+    # Only emit one instance of the decode_timedelta default change
+    # FutureWarning. This can be removed once this change is made.
+    warnings.filterwarnings("once", "decode_timedelta", FutureWarning)
+
     dimensions_used_by = defaultdict(list)
     for v in variables.values():
         for d in v.dims:
@@ -472,7 +491,10 @@ def decode_cf(
     decode_coords: bool | Literal["coordinates", "all"] = True,
     drop_variables: T_DropVariables = None,
     use_cftime: bool | None = None,
-    decode_timedelta: bool | None = None,
+    decode_timedelta: bool
+    | CFTimedeltaCoder
+    | Mapping[str, bool | CFTimedeltaCoder]
+    | None = None,
 ) -> Dataset:
     """Decode the given Dataset or Datastore according to CF conventions into
     a new Dataset.
@@ -516,11 +538,14 @@ def decode_cf(
         .. deprecated:: 2025.01.1
            Please pass a :py:class:`coders.CFDatetimeCoder` instance initialized with ``use_cftime`` to the ``decode_times`` kwarg instead.
 
-    decode_timedelta : bool, optional
-        If True, decode variables and coordinates with time units in
+    decode_timedelta : bool | CFTimedeltaCoder | Mapping[str, bool | CFTimedeltaCoder], optional
+        If True or :py:class:`CFTimedeltaCoder`, 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.
+        If None (default), assume the same behavior as decode_times. The
+        resolution of the decoded timedeltas can be configured with the
+        ``time_unit`` argument in the :py:class:`CFTimedeltaCoder` passed.
 
     Returns
     -------
diff --git a/xarray/core/groupby.py b/xarray/core/groupby.py
index 291632f9..ceae7903 100644
--- a/xarray/core/groupby.py
+++ b/xarray/core/groupby.py
@@ -1357,12 +1357,7 @@ class GroupBy(Generic[T_Xarray]):
         """
         return ops.where_method(self, cond, other)
 
-    def _first_or_last(
-        self,
-        op: Literal["first" | "last"],
-        skipna: bool | None,
-        keep_attrs: bool | None,
-    ):
+    def _first_or_last(self, op, skipna, keep_attrs):
         if all(
             isinstance(maybe_slice, slice)
             and (maybe_slice.stop == maybe_slice.start + 1)
@@ -1373,30 +1368,17 @@ class GroupBy(Generic[T_Xarray]):
             return self._obj
         if keep_attrs is None:
             keep_attrs = _get_keep_attrs(default=True)
-        if (
-            module_available("flox", minversion="0.9.16")
-            and OPTIONS["use_flox"]
-            and contains_only_chunked_or_numpy(self._obj)
-        ):
-            result, *_ = self._flox_reduce(
-                dim=None, func=op, skipna=skipna, keep_attrs=keep_attrs
-            )
-        else:
-            result = self.reduce(
-                getattr(duck_array_ops, op),
-                dim=[self._group_dim],
-                skipna=skipna,
-                keep_attrs=keep_attrs,
-            )
-        return result
+        return self.reduce(
+            op, dim=[self._group_dim], skipna=skipna, keep_attrs=keep_attrs
+        )
 
     def first(self, skipna: bool | None = None, keep_attrs: bool | None = None):
         """Return the first element of each group along the group dimension"""
-        return self._first_or_last("first", skipna, keep_attrs)
+        return self._first_or_last(duck_array_ops.first, skipna, keep_attrs)
 
     def last(self, skipna: bool | None = None, keep_attrs: bool | None = None):
         """Return the last element of each group along the group dimension"""
-        return self._first_or_last("last", skipna, keep_attrs)
+        return self._first_or_last(duck_array_ops.last, skipna, keep_attrs)
 
     def assign_coords(self, coords=None, **coords_kwargs):
         """Assign coordinates by group.
diff --git a/xarray/core/resample.py b/xarray/core/resample.py
index ebd3d46e..5cc98c96 100644
--- a/xarray/core/resample.py
+++ b/xarray/core/resample.py
@@ -2,7 +2,7 @@ from __future__ import annotations
 
 import warnings
 from collections.abc import Callable, Hashable, Iterable, Sequence
-from typing import TYPE_CHECKING, Any, Literal
+from typing import TYPE_CHECKING, Any
 
 from xarray.core._aggregations import (
     DataArrayResampleAggregations,
@@ -103,21 +103,6 @@ class Resample(GroupBy[T_Xarray]):
         (grouper,) = self.groupers
         return self._shuffle_obj(chunks).drop_vars(RESAMPLE_DIM)
 
-    def _first_or_last(
-        self, op: Literal["first", "last"], skipna: bool | None, keep_attrs: bool | None
-    ) -> T_Xarray:
-        from xarray.core.dataset import Dataset
-
-        result = super()._first_or_last(op=op, skipna=skipna, keep_attrs=keep_attrs)
-        if isinstance(result, Dataset):
-            # Can't do this in the base class because group_dim is RESAMPLE_DIM
-            # which is not present in the original object
-            for var in result.data_vars:
-                result._variables[var] = result._variables[var].transpose(
-                    *self._obj._variables[var].dims
-                )
-        return result
-
     def _drop_coords(self) -> T_Xarray:
         """Drop non-dimension coordinates along the resampled dimension."""
         obj = self._obj


Test patch:
diff --git a/xarray/tests/test_backends.py b/xarray/tests/test_backends.py
index 420e30b8..bc6e16ff 100644
--- a/xarray/tests/test_backends.py
+++ b/xarray/tests/test_backends.py
@@ -49,7 +49,7 @@ from xarray.backends.netCDF4_ import (
 from xarray.backends.pydap_ import PydapDataStore
 from xarray.backends.scipy_ import ScipyBackendEntrypoint
 from xarray.backends.zarr import ZarrStore
-from xarray.coders import CFDatetimeCoder
+from xarray.coders import CFDatetimeCoder, CFTimedeltaCoder
 from xarray.coding.cftime_offsets import cftime_range
 from xarray.coding.strings import check_vlen_dtype, create_vlen_dtype
 from xarray.coding.variables import SerializationWarning
@@ -636,7 +636,9 @@ class DatasetIOBase:
         #  to support large ranges
         time_deltas = pd.to_timedelta(["1h", "2h", "NaT"]).as_unit("s")  # type: ignore[arg-type, unused-ignore]
         expected = Dataset({"td": ("td", time_deltas), "td0": time_deltas[0]})
-        with self.roundtrip(expected) as actual:
+        with self.roundtrip(
+            expected, open_kwargs={"decode_timedelta": CFTimedeltaCoder(time_unit="ns")}
+        ) as actual:
             assert_identical(expected, actual)
 
     def test_roundtrip_float64_data(self) -> None:
@@ -3282,7 +3284,13 @@ class ZarrBase(CFEncodedBase):
     def test_chunked_datetime64_or_timedelta64(self, dtype) -> None:
         # Generalized from @malmans2's test in PR #8253
         original = create_test_data().astype(dtype).chunk(1)
-        with self.roundtrip(original, open_kwargs={"chunks": {}}) as actual:
+        with self.roundtrip(
+            original,
+            open_kwargs={
+                "chunks": {},
+                "decode_timedelta": CFTimedeltaCoder(time_unit="ns"),
+            },
+        ) as actual:
             for name, actual_var in actual.variables.items():
                 assert original[name].chunks == actual_var.chunks
             assert original.chunks == actual.chunks
diff --git a/xarray/tests/test_coding_times.py b/xarray/tests/test_coding_times.py
index cdf97f08..6e1eda28 100644
--- a/xarray/tests/test_coding_times.py
+++ b/xarray/tests/test_coding_times.py
@@ -19,7 +19,7 @@ from xarray import (
     date_range,
     decode_cf,
 )
-from xarray.coders import CFDatetimeCoder
+from xarray.coders import CFDatetimeCoder, CFTimedeltaCoder
 from xarray.coding.times import _STANDARD_CALENDARS as _STANDARD_CALENDARS_UNSORTED
 from xarray.coding.times import (
     _encode_datetime_with_cftime,
@@ -1513,7 +1513,10 @@ def test_roundtrip_timedelta64_nanosecond_precision(
 
     encoded_var = conventions.encode_cf_variable(var)
     decoded_var = conventions.decode_cf_variable(
-        "foo", encoded_var, decode_times=CFDatetimeCoder(time_unit=time_unit)
+        "foo",
+        encoded_var,
+        decode_times=CFDatetimeCoder(time_unit=time_unit),
+        decode_timedelta=CFTimedeltaCoder(time_unit=time_unit),
     )
 
     assert_identical(var, decoded_var)
@@ -1540,7 +1543,9 @@ def test_roundtrip_timedelta64_nanosecond_precision_warning() -> None:
     assert encoded_var.dtype == np.int64
     assert encoded_var.attrs["units"] == needed_units
     assert encoded_var.attrs["_FillValue"] == 20
-    decoded_var = conventions.decode_cf_variable("foo", encoded_var)
+    decoded_var = conventions.decode_cf_variable(
+        "foo", encoded_var, decode_timedelta=CFTimedeltaCoder(time_unit="ns")
+    )
     assert_identical(var, decoded_var)
     assert decoded_var.encoding["dtype"] == np.int64
 
@@ -1588,7 +1593,9 @@ def test_roundtrip_float_times(fill_value, times, units, encoded_values) -> None
     assert encoded_var.attrs["units"] == units
     assert encoded_var.attrs["_FillValue"] == fill_value
 
-    decoded_var = conventions.decode_cf_variable("foo", encoded_var)
+    decoded_var = conventions.decode_cf_variable(
+        "foo", encoded_var, decode_timedelta=CFTimedeltaCoder(time_unit="ns")
+    )
     assert_identical(var, decoded_var)
     assert decoded_var.encoding["units"] == units
     assert decoded_var.encoding["_FillValue"] == fill_value
@@ -1779,7 +1786,9 @@ def test_encode_cf_timedelta_casting_value_error(use_dask) -> None:
         with pytest.warns(UserWarning, match="Timedeltas can't be serialized"):
             encoded = conventions.encode_cf_variable(variable)
         assert encoded.attrs["units"] == "hours"
-        decoded = conventions.decode_cf_variable("name", encoded)
+        decoded = conventions.decode_cf_variable(
+            "name", encoded, decode_timedelta=CFTimedeltaCoder(time_unit="ns")
+        )
         assert_equal(variable, decoded)
     else:
         with pytest.raises(ValueError, match="Not possible"):
@@ -1800,3 +1809,88 @@ def test_encode_cf_timedelta_casting_overflow_error(use_dask, dtype) -> None:
     with pytest.raises(OverflowError, match="Not possible"):
         encoded = conventions.encode_cf_variable(variable)
         encoded.compute()
+
+
+_DECODE_TIMEDELTA_TESTS = {
+    "default": (True, None, np.dtype("timedelta64[ns]"), True),
+    "decode_timdelta=False": (True, False, np.dtype("int64"), False),
+    "inherit-time_unit-from-decode_times": (
+        CFDatetimeCoder(time_unit="s"),
+        None,
+        np.dtype("timedelta64[s]"),
+        True,
+    ),
+    "set-time_unit-via-CFTimedeltaCoder-decode_times=True": (
+        True,
+        CFTimedeltaCoder(time_unit="s"),
+        np.dtype("timedelta64[s]"),
+        False,
+    ),
+    "set-time_unit-via-CFTimedeltaCoder-decode_times=False": (
+        False,
+        CFTimedeltaCoder(time_unit="s"),
+        np.dtype("timedelta64[s]"),
+        False,
+    ),
+    "override-time_unit-from-decode_times": (
+        CFDatetimeCoder(time_unit="ns"),
+        CFTimedeltaCoder(time_unit="s"),
+        np.dtype("timedelta64[s]"),
+        False,
+    ),
+}
+
+
+@pytest.mark.parametrize(
+    ("decode_times", "decode_timedelta", "expected_dtype", "warns"),
+    list(_DECODE_TIMEDELTA_TESTS.values()),
+    ids=list(_DECODE_TIMEDELTA_TESTS.keys()),
+)
+def test_decode_timedelta(
+    decode_times, decode_timedelta, expected_dtype, warns
+) -> None:
+    timedeltas = pd.timedelta_range(0, freq="d", periods=3)
+    var = Variable(["time"], timedeltas)
+    encoded = conventions.encode_cf_variable(var)
+    if warns:
+        with pytest.warns(FutureWarning, match="decode_timedelta"):
+            decoded = conventions.decode_cf_variable(
+                "foo",
+                encoded,
+                decode_times=decode_times,
+                decode_timedelta=decode_timedelta,
+            )
+    else:
+        decoded = conventions.decode_cf_variable(
+            "foo", encoded, decode_times=decode_times, decode_timedelta=decode_timedelta
+        )
+    if decode_timedelta is False:
+        assert_equal(encoded, decoded)
+    else:
+        assert_equal(var, decoded)
+    assert decoded.dtype == expected_dtype
+
+
+def test_lazy_decode_timedelta_unexpected_dtype() -> None:
+    attrs = {"units": "seconds"}
+    encoded = Variable(["time"], [0, 0.5, 1], attrs=attrs)
+    decoded = conventions.decode_cf_variable(
+        "foo", encoded, decode_timedelta=CFTimedeltaCoder(time_unit="s")
+    )
+
+    expected_dtype_upon_lazy_decoding = np.dtype("timedelta64[s]")
+    assert decoded.dtype == expected_dtype_upon_lazy_decoding
+
+    expected_dtype_upon_loading = np.dtype("timedelta64[ms]")
+    with pytest.warns(SerializationWarning, match="Can't decode floating"):
+        assert decoded.load().dtype == expected_dtype_upon_loading
+
+
+def test_lazy_decode_timedelta_error() -> None:
+    attrs = {"units": "seconds"}
+    encoded = Variable(["time"], [0, np.iinfo(np.int64).max, 1], attrs=attrs)
+    decoded = conventions.decode_cf_variable(
+        "foo", encoded, decode_timedelta=CFTimedeltaCoder(time_unit="ms")
+    )
+    with pytest.raises(OutOfBoundsTimedelta, match="overflow"):
+        decoded.load()
diff --git a/xarray/tests/test_conventions.py b/xarray/tests/test_conventions.py
index 346ad1c9..8d3827fa 100644
--- a/xarray/tests/test_conventions.py
+++ b/xarray/tests/test_conventions.py
@@ -18,7 +18,7 @@ from xarray import (
 )
 from xarray.backends.common import WritableCFDataStore
 from xarray.backends.memory import InMemoryDataStore
-from xarray.coders import CFDatetimeCoder
+from xarray.coders import CFDatetimeCoder, CFTimedeltaCoder
 from xarray.conventions import decode_cf
 from xarray.testing import assert_identical
 from xarray.tests import (
@@ -536,9 +536,11 @@ class TestDecodeCF:
         )
 
         dsc = conventions.decode_cf(
-            ds, decode_times=CFDatetimeCoder(time_unit=time_unit)
+            ds,
+            decode_times=CFDatetimeCoder(time_unit=time_unit),
+            decode_timedelta=CFTimedeltaCoder(time_unit=time_unit),
         )
-        assert dsc.timedelta.dtype == np.dtype("m8[ns]")
+        assert dsc.timedelta.dtype == np.dtype(f"m8[{time_unit}]")
         assert dsc.time.dtype == np.dtype(f"M8[{time_unit}]")
         dsc = conventions.decode_cf(ds, decode_times=False)
         assert dsc.timedelta.dtype == np.dtype("int64")
@@ -655,3 +657,15 @@ def test_encode_cf_variable_with_vlen_dtype() -> None:
     encoded_v = conventions.encode_cf_variable(v)
     assert encoded_v.data.dtype.kind == "O"
     assert coding.strings.check_vlen_dtype(encoded_v.data.dtype) is str
+
+
+def test_decode_cf_variables_decode_timedelta_warning() -> None:
+    v = Variable(["time"], [1, 2], attrs={"units": "seconds"})
+    variables = {"a": v}
+
+    with warnings.catch_warnings():
+        warnings.filterwarnings("error", "decode_timedelta", FutureWarning)
+        conventions.decode_cf_variables(variables, {}, decode_timedelta=True)
+
+    with pytest.warns(FutureWarning, match="decode_timedelta"):
+        conventions.decode_cf_variables(variables, {})
diff --git a/xarray/tests/test_groupby.py b/xarray/tests/test_groupby.py
index 2287130e..d42f86f5 100644
--- a/xarray/tests/test_groupby.py
+++ b/xarray/tests/test_groupby.py
@@ -1618,8 +1618,6 @@ class TestDataArrayGroupBy:
         expected = array  # should be a no-op
         assert_identical(expected, actual)
 
-        # TODO: groupby_bins too
-
     def make_groupby_multidim_example_array(self) -> DataArray:
         return DataArray(
             [[[0, 1], [2, 3]], [[5, 10], [15, 20]]],