Repository: pydata/xarray

Issue:
Averaging timestamps with non-nanosecond precision
### What happened?

The latest xarray version from the main branch returns unexpected values when averaging timestamps with non-nanosecond precision.

### What did you expect to happen?

The average of [01:00, 02:00] to be 01:30, instead it is 01:00:00.000001800

### Minimal Complete Verifiable Example

```Python
import xarray as xr
import numpy as np

timestamps = xr.DataArray(
    [np.datetime64("1970-01-01 01:00"),
     np.datetime64("1970-01-01 02:00")]
)
timestamps.mean()
```

### MVCE confirmation

- [x] Minimal example — the example is as focused as reasonably possible to demonstrate the underlying issue in xarray.
- [x] Complete example — the example is self-contained, including all data and the text of any traceback.
- [x] Verifiable example — the example copy & pastes into an IPython prompt or [Binder notebook](https://mybinder.org/v2/gh/pydata/xarray/main?urlpath=lab/tree/doc/examples/blank_template.ipynb), returning the result.
- [x] New issue — a search of GitHub Issues suggests this is not a duplicate.
- [x] Recent environment — the issue occurs with the latest version of xarray and its dependencies.

### Relevant log output

```Python
<xarray.DataArray ()> Size: 8B
array('1970-01-01T01:00:00.000001800', dtype='datetime64[ns]')
```

### Anything else we need to know?

_No response_

### Environment

<details>

INSTALLED VERSIONS
------------------
commit: None
python: 3.12.8 | packaged by conda-forge | (main, Dec  5 2024, 14:24:40) [GCC 13.3.0]
python-bits: 64
OS: Linux
OS-release: 4.18.0-513.5.1.el8_9.x86_64
machine: x86_64
processor: x86_64
byteorder: little
LC_ALL: en_US.utf-8
LANG: de_DE.UTF-8
LOCALE: ('en_US', 'UTF-8')
libhdf5: 1.14.4
libnetcdf: 4.9.2
xarray: 2025.1.2.dev9+g609412d8
pandas: 2.2.3
numpy: 2.2.2
scipy: 1.15.1
netCDF4: 1.7.2
pydap: None
h5netcdf: 1.4.1
h5py: 3.12.1
zarr: 3.0.1
cftime: 1.6.4
nc_time_axis: None
iris: None
bottleneck: None
dask: 2024.12.1
distributed: None
matplotlib: None
cartopy: None
seaborn: None
numbagg: None
fsspec: 2024.12.0
cupy: None
pint: None
sparse: None
flox: None
numpy_groupies: None
setuptools: 75.8.0
pip: 24.3.1
conda: None
pytest: 8.3.4
mypy: None
IPython: None
sphinx: None

</details>


Discussion:
kmuehlbauer: > Works fine, thanks a lot!

Great, let's wait for some more feedback. @dcherian and @spencerkclark would you mind taking a look here?
kmuehlbauer: Apart from the loss of resolution when the calculated mean contains fractional part (but is in alignment with pandas and numpy), this is working for all available resolutions. In case of `ns` resolution this replicates the old behaviour.  
kmuehlbauer: > Thanks @kmuehlbauer—this looks good to me now!

Thanks @spencerkclark, good to have a decent solution here for the time being.

Patch:
diff --git a/.github/workflows/ci-additional.yaml b/.github/workflows/ci-additional.yaml
index 96eca004..e2fa49f2 100644
--- a/.github/workflows/ci-additional.yaml
+++ b/.github/workflows/ci-additional.yaml
@@ -123,7 +123,7 @@ jobs:
           python -m mypy --install-types --non-interactive --cobertura-xml-report mypy_report
 
       - name: Upload mypy coverage to Codecov
-        uses: codecov/codecov-action@v5.1.2
+        uses: codecov/codecov-action@v5.3.1
         with:
           file: mypy_report/cobertura.xml
           flags: mypy
@@ -174,7 +174,7 @@ jobs:
           python -m mypy --install-types --non-interactive --cobertura-xml-report mypy_report
 
       - name: Upload mypy coverage to Codecov
-        uses: codecov/codecov-action@v5.1.2
+        uses: codecov/codecov-action@v5.3.1
         with:
           file: mypy_report/cobertura.xml
           flags: mypy-min
@@ -230,7 +230,7 @@ jobs:
           python -m pyright xarray/
 
       - name: Upload pyright coverage to Codecov
-        uses: codecov/codecov-action@v5.1.2
+        uses: codecov/codecov-action@v5.3.1
         with:
           file: pyright_report/cobertura.xml
           flags: pyright
@@ -286,7 +286,7 @@ jobs:
           python -m pyright xarray/
 
       - name: Upload pyright coverage to Codecov
-        uses: codecov/codecov-action@v5.1.2
+        uses: codecov/codecov-action@v5.3.1
         with:
           file: pyright_report/cobertura.xml
           flags: pyright39
diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml
index 6d3f7d63..e296951c 100644
--- a/.github/workflows/ci.yaml
+++ b/.github/workflows/ci.yaml
@@ -159,7 +159,7 @@ jobs:
           path: pytest.xml
 
       - name: Upload code coverage to Codecov
-        uses: codecov/codecov-action@v5.1.2
+        uses: codecov/codecov-action@v5.3.1
         env:
           CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }}
         with:
diff --git a/.github/workflows/pypi-release.yaml b/.github/workflows/pypi-release.yaml
index 6377f59a..add9261f 100644
--- a/.github/workflows/pypi-release.yaml
+++ b/.github/workflows/pypi-release.yaml
@@ -88,7 +88,7 @@ jobs:
           path: dist
       - name: Publish package to TestPyPI
         if: github.event_name == 'push'
-        uses: pypa/gh-action-pypi-publish@v1.12.3
+        uses: pypa/gh-action-pypi-publish@v1.12.4
         with:
           repository_url: https://test.pypi.org/legacy/
           verbose: true
@@ -110,6 +110,6 @@ jobs:
           name: releases
           path: dist
       - name: Publish package to PyPI
-        uses: pypa/gh-action-pypi-publish@v1.12.3
+        uses: pypa/gh-action-pypi-publish@v1.12.4
         with:
           verbose: true
diff --git a/.github/workflows/upstream-dev-ci.yaml b/.github/workflows/upstream-dev-ci.yaml
index f338a222..18a0aca0 100644
--- a/.github/workflows/upstream-dev-ci.yaml
+++ b/.github/workflows/upstream-dev-ci.yaml
@@ -140,7 +140,7 @@ jobs:
         run: |
           python -m mypy --install-types --non-interactive --cobertura-xml-report mypy_report
       - name: Upload mypy coverage to Codecov
-        uses: codecov/codecov-action@v5.1.2
+        uses: codecov/codecov-action@v5.3.1
         with:
           file: mypy_report/cobertura.xml
           flags: mypy
diff --git a/doc/ecosystem.rst b/doc/ecosystem.rst
index d5123669..1f1ca04b 100644
--- a/doc/ecosystem.rst
+++ b/doc/ecosystem.rst
@@ -38,6 +38,7 @@ Geosciences
 - `salem <https://salem.readthedocs.io>`_: Adds geolocalised subsetting, masking, and plotting operations to xarray's data structures via accessors.
 - `SatPy <https://satpy.readthedocs.io/>`_ : Library for reading and manipulating meteorological remote sensing data and writing it to various image and data file formats.
 - `SARXarray <https://tudelftgeodesy.github.io/sarxarray/>`_: xarray extension for reading and processing large Synthetic Aperture Radar (SAR) data stacks.
+- `shxarray <https://shxarray.wobbly.earth/>`_: Convert, filter,and map geodesy related spherical harmonic representations of gravity and terrestrial water storage through an xarray extension.
 - `Spyfit <https://spyfit.readthedocs.io/en/master/>`_: FTIR spectroscopy of the atmosphere
 - `windspharm <https://ajdawson.github.io/windspharm/index.html>`_: Spherical
   harmonic wind analysis in Python.
diff --git a/doc/user-guide/pandas.rst b/doc/user-guide/pandas.rst
index 3d4a6b42..9e070ae6 100644
--- a/doc/user-guide/pandas.rst
+++ b/doc/user-guide/pandas.rst
@@ -78,7 +78,7 @@ To create a ``Dataset`` from a ``DataFrame``, use the
 
     xr.Dataset.from_dataframe(df)
 
-Notice that that dimensions of variables in the ``Dataset`` have now
+Notice that the dimensions of variables in the ``Dataset`` have now
 expanded after the round-trip conversion to a ``DataFrame``. This is because
 every object in a ``DataFrame`` must have the same indices, so we need to
 broadcast the data of each array to the full size of the new ``MultiIndex``.
diff --git a/doc/whats-new.rst b/doc/whats-new.rst
index 9b40a323..75ec2273 100644
--- a/doc/whats-new.rst
+++ b/doc/whats-new.rst
@@ -50,10 +50,21 @@ eventually be deprecated.
 
 New Features
 ~~~~~~~~~~~~
-- Relax nanosecond datetime restriction in CF time decoding (:issue:`7493`, :pull:`9618`).
+- Relax nanosecond datetime restriction in CF time decoding (:issue:`7493`, :pull:`9618`, :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>`_.
 - Improve the error message raised when no key is matching the available variables in a dataset.  (:pull:`9943`)
   By `Jimmy Westling <https://github.com/illviljan>`_.
+- Added a ``time_unit`` argument to :py:meth:`CFTimeIndex.to_datetimeindex`.
+  Note that in a future version of xarray,
+  :py:meth:`CFTimeIndex.to_datetimeindex` will return a microsecond-resolution
+  :py:class:`pandas.DatetimeIndex` instead of a nanosecond-resolution
+  :py:class:`pandas.DatetimeIndex` (:pull:`9965`). By `Spencer Clark
+  <https://github.com/spencerkclark>`_ and `Kai Mühlbauer
+  <https://github.com/kmuehlbauer>`_.
+- :py:meth:`DatasetGroupBy.first` and :py:meth:`DatasetGroupBy.last` can now use ``flox`` if available. (:issue:`9647`)
+  By `Deepak Cherian <https://github.com/dcherian>`_.
 
 Breaking changes
 ~~~~~~~~~~~~~~~~
diff --git a/xarray/backends/zarr.py b/xarray/backends/zarr.py
index a584f78e..e83f5556 100644
--- a/xarray/backends/zarr.py
+++ b/xarray/backends/zarr.py
@@ -444,6 +444,7 @@ def extract_zarr_variable_encoding(
     safe_to_drop = {"source", "original_shape", "preferred_chunks"}
     valid_encodings = {
         "chunks",
+        "shards",
         "compressor",  # TODO: delete when min zarr >=3
         "compressors",
         "filters",
@@ -825,6 +826,7 @@ class ZarrStore(AbstractWritableDataStore):
                 {
                     "compressors": zarr_array.compressors,
                     "filters": zarr_array.filters,
+                    "shards": zarr_array.shards,
                 }
             )
             if self.zarr_group.metadata.zarr_format == 3:
diff --git a/xarray/coding/cftimeindex.py b/xarray/coding/cftimeindex.py
index bd5f5155..34b266b9 100644
--- a/xarray/coding/cftimeindex.py
+++ b/xarray/coding/cftimeindex.py
@@ -42,9 +42,8 @@
 from __future__ import annotations
 
 import math
-import warnings
 from datetime import timedelta
-from typing import TYPE_CHECKING, Any
+from typing import TYPE_CHECKING, Any, Optional
 
 import numpy as np
 import pandas as pd
@@ -58,7 +57,8 @@ from xarray.coding.times import (
 )
 from xarray.core.common import _contains_cftime_datetimes
 from xarray.core.options import OPTIONS
-from xarray.core.utils import attempt_import, is_scalar
+from xarray.core.types import PDDatetimeUnitOptions
+from xarray.core.utils import attempt_import, emit_user_level_warning, is_scalar
 
 if TYPE_CHECKING:
     from xarray.coding.cftime_offsets import BaseCFTimeOffset
@@ -239,6 +239,8 @@ class CFTimeIndex(pd.Index):
     cftime_range
     """
 
+    _data: np.ndarray
+
     year = _field_accessor("year", "The year of the datetime")
     month = _field_accessor("month", "The month of the datetime")
     day = _field_accessor("day", "The days of the datetime")
@@ -544,14 +546,18 @@ class CFTimeIndex(pd.Index):
                 "that can be expressed at the nanosecond resolution."
             ) from err
 
-    def to_datetimeindex(self, unsafe=False):
+    def to_datetimeindex(
+        self, unsafe: bool = False, time_unit: Optional[PDDatetimeUnitOptions] = None
+    ) -> pd.DatetimeIndex:
         """If possible, convert this index to a pandas.DatetimeIndex.
 
         Parameters
         ----------
         unsafe : bool
-            Flag to turn off warning when converting from a CFTimeIndex with
-            a non-standard calendar to a DatetimeIndex (default ``False``).
+            Flag to turn off calendar mismatch warnings (default ``False``).
+        time_unit : str
+            Time resolution of resulting DatetimeIndex. Can be one of `"s"`,
+            ``"ms"``, ``"us"``, or ``"ns"`` (default ``"ns"``).
 
         Returns
         -------
@@ -561,18 +567,20 @@ class CFTimeIndex(pd.Index):
         ------
         ValueError
             If the CFTimeIndex contains dates that are not possible in the
-            standard calendar or outside the nanosecond-precision range.
+            standard calendar or outside the range representable by the
+            specified ``time_unit``.
 
         Warns
         -----
         RuntimeWarning
-            If converting from a non-standard calendar to a DatetimeIndex.
+            If converting from a non-standard calendar, or a Gregorian
+            calendar with dates prior to the reform (1582-10-15).
 
         Warnings
         --------
-        Note that for non-standard calendars, this will change the calendar
-        type of the index.  In that case the result of this method should be
-        used with caution.
+        Note that for non-proleptic Gregorian calendars, this will change the
+        calendar type of the index. In that case the result of this method
+        should be used with caution.
 
         Examples
         --------
@@ -580,26 +588,47 @@ class CFTimeIndex(pd.Index):
         >>> times
         CFTimeIndex([2000-01-01 00:00:00, 2000-01-02 00:00:00],
                     dtype='object', length=2, calendar='standard', freq=None)
-        >>> times.to_datetimeindex()
-        DatetimeIndex(['2000-01-01', '2000-01-02'], dtype='datetime64[us]', freq=None)
+        >>> times.to_datetimeindex(time_unit="ns")
+        DatetimeIndex(['2000-01-01', '2000-01-02'], dtype='datetime64[ns]', freq=None)
         """
 
         if not self._data.size:
             return pd.DatetimeIndex([])
 
-        # transform to us-resolution is needed for DatetimeIndex
-        nptimes = cftime_to_nptime(self, time_unit="us")
+        if time_unit is None:
+            emit_user_level_warning(
+                "In a future version of xarray to_datetimeindex will default "
+                "to returning a 'us'-resolution DatetimeIndex instead of a "
+                "'ns'-resolution DatetimeIndex. This warning can be silenced "
+                "by explicitly passing the `time_unit` keyword argument.",
+                FutureWarning,
+            )
+            time_unit = "ns"
+
+        nptimes = cftime_to_nptime(self, time_unit=time_unit)
         calendar = infer_calendar_name(self)
         if calendar not in _STANDARD_CALENDARS and not unsafe:
-            warnings.warn(
+            emit_user_level_warning(
                 "Converting a CFTimeIndex with dates from a non-standard "
-                f"calendar, {calendar!r}, to a pandas.DatetimeIndex, which uses dates "
-                "from the standard calendar.  This may lead to subtle errors "
-                "in operations that depend on the length of time between "
-                "dates.",
+                f"calendar, {calendar!r}, to a pandas.DatetimeIndex, which "
+                "uses dates from the standard calendar.  This may lead to "
+                "subtle errors in operations that depend on the length of "
+                "time between dates.",
                 RuntimeWarning,
-                stacklevel=2,
             )
+        if calendar == "standard" and not unsafe:
+            reform_date = self.date_type(1582, 10, 15)
+            if self.min() < reform_date:
+                emit_user_level_warning(
+                    "Converting a CFTimeIndex with dates from a Gregorian "
+                    "calendar that fall before the reform date of 1582-10-15 "
+                    "to a pandas.DatetimeIndex. During this time period the "
+                    "Gregorian calendar and the proleptic Gregorian calendar "
+                    "of the DatetimeIndex do not exactly align. This warning "
+                    "can be silenced by setting unsafe=True.",
+                    RuntimeWarning,
+                )
+
         return pd.DatetimeIndex(nptimes)
 
     def strftime(self, date_format):
diff --git a/xarray/core/datatree_io.py b/xarray/core/datatree_io.py
index 3a57ebc7..31f5a93d 100644
--- a/xarray/core/datatree_io.py
+++ b/xarray/core/datatree_io.py
@@ -87,7 +87,7 @@ def _datatree_to_zarr(
     consolidated: bool = True,
     group: str | None = None,
     write_inherited_coords: bool = False,
-    compute: Literal[True] = True,
+    compute: bool = True,
     **kwargs,
 ):
     """This function creates an appropriate datastore for writing a datatree
@@ -103,9 +103,6 @@ def _datatree_to_zarr(
             "specifying a root group for the tree has not been implemented"
         )
 
-    if not compute:
-        raise NotImplementedError("compute=False has not been implemented yet")
-
     if encoding is None:
         encoding = {}
 
@@ -127,6 +124,7 @@ def _datatree_to_zarr(
             mode=mode,
             encoding=encoding.get(node.path),
             consolidated=False,
+            compute=compute,
             **kwargs,
         )
         if "w" in mode:
diff --git a/xarray/core/duck_array_ops.py b/xarray/core/duck_array_ops.py
index 7e7333fd..c3f15980 100644
--- a/xarray/core/duck_array_ops.py
+++ b/xarray/core/duck_array_ops.py
@@ -662,16 +662,10 @@ def _to_pytimedelta(array, unit="us"):
 
 
 def np_timedelta64_to_float(array, datetime_unit):
-    """Convert numpy.timedelta64 to float.
-
-    Notes
-    -----
-    The array is first converted to microseconds, which is less likely to
-    cause overflow errors.
-    """
-    array = array.astype("timedelta64[ns]").astype(np.float64)
-    conversion_factor = np.timedelta64(1, "ns") / np.timedelta64(1, datetime_unit)
-    return conversion_factor * array
+    """Convert numpy.timedelta64 to float, possibly at a loss of resolution."""
+    unit, _ = np.datetime_data(array.dtype)
+    conversion_factor = np.timedelta64(1, unit) / np.timedelta64(1, datetime_unit)
+    return conversion_factor * array.astype(np.float64)
 
 
 def pd_timedelta_to_float(value, datetime_unit):
@@ -715,12 +709,15 @@ def mean(array, axis=None, skipna=None, **kwargs):
     if dtypes.is_datetime_like(array.dtype):
         offset = _datetime_nanmin(array)
 
-        # xarray always uses np.datetime64[ns] for np.datetime64 data
-        dtype = "timedelta64[ns]"
+        # From version 2025.01.2 xarray uses np.datetime64[unit], where unit
+        # is one of "s", "ms", "us", "ns".
+        # To not have to worry about the resolution, we just convert the output
+        # to "timedelta64" (without unit) and let the dtype of offset take precedence.
+        # This is fully backwards compatible with datetime64[ns].
         return (
             _mean(
                 datetime_to_numeric(array, offset), axis=axis, skipna=skipna, **kwargs
-            ).astype(dtype)
+            ).astype("timedelta64")
             + offset
         )
     elif _contains_cftime_datetimes(array):
diff --git a/xarray/core/groupby.py b/xarray/core/groupby.py
index ceae7903..291632f9 100644
--- a/xarray/core/groupby.py
+++ b/xarray/core/groupby.py
@@ -1357,7 +1357,12 @@ class GroupBy(Generic[T_Xarray]):
         """
         return ops.where_method(self, cond, other)
 
-    def _first_or_last(self, op, skipna, keep_attrs):
+    def _first_or_last(
+        self,
+        op: Literal["first" | "last"],
+        skipna: bool | None,
+        keep_attrs: bool | None,
+    ):
         if all(
             isinstance(maybe_slice, slice)
             and (maybe_slice.stop == maybe_slice.start + 1)
@@ -1368,17 +1373,30 @@ class GroupBy(Generic[T_Xarray]):
             return self._obj
         if keep_attrs is None:
             keep_attrs = _get_keep_attrs(default=True)
-        return self.reduce(
-            op, dim=[self._group_dim], skipna=skipna, keep_attrs=keep_attrs
-        )
+        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
 
     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(duck_array_ops.first, skipna, keep_attrs)
+        return self._first_or_last("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(duck_array_ops.last, skipna, keep_attrs)
+        return self._first_or_last("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 5cc98c96..ebd3d46e 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
+from typing import TYPE_CHECKING, Any, Literal
 
 from xarray.core._aggregations import (
     DataArrayResampleAggregations,
@@ -103,6 +103,21 @@ 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_accessor_dt.py b/xarray/tests/test_accessor_dt.py
index 64309966..652b1752 100644
--- a/xarray/tests/test_accessor_dt.py
+++ b/xarray/tests/test_accessor_dt.py
@@ -519,7 +519,9 @@ def test_cftime_strftime_access(data) -> None:
     date_format = "%Y%m%d%H"
     result = data.time.dt.strftime(date_format)
     datetime_array = xr.DataArray(
-        xr.coding.cftimeindex.CFTimeIndex(data.time.values).to_datetimeindex(),
+        xr.coding.cftimeindex.CFTimeIndex(data.time.values).to_datetimeindex(
+            time_unit="ns"
+        ),
         name="stftime",
         coords=data.time.coords,
         dims=data.time.dims,
diff --git a/xarray/tests/test_backends.py b/xarray/tests/test_backends.py
index dc333719..420e30b8 100644
--- a/xarray/tests/test_backends.py
+++ b/xarray/tests/test_backends.py
@@ -2496,6 +2496,24 @@ class ZarrBase(CFEncodedBase):
             with self.roundtrip(data) as actual:
                 pass
 
+    def test_shard_encoding(self) -> None:
+        # These datasets have no dask chunks. All chunking/sharding specified in
+        # encoding
+        if has_zarr_v3 and zarr.config.config["default_zarr_format"] == 3:
+            data = create_test_data()
+            chunks = (1, 1)
+            shards = (5, 5)
+            data["var2"].encoding.update({"chunks": chunks})
+            data["var2"].encoding.update({"shards": shards})
+            with self.roundtrip(data) as actual:
+                assert shards == actual["var2"].encoding["shards"]
+
+            # expect an error with shards not divisible by chunks
+            data["var2"].encoding.update({"chunks": (2, 2)})
+            with pytest.raises(ValueError):
+                with self.roundtrip(data) as actual:
+                    pass
+
     @requires_dask
     @pytest.mark.skipif(
         ON_WINDOWS,
diff --git a/xarray/tests/test_backends_datatree.py b/xarray/tests/test_backends_datatree.py
index 9cbe3f22..bca528ca 100644
--- a/xarray/tests/test_backends_datatree.py
+++ b/xarray/tests/test_backends_datatree.py
@@ -440,6 +440,25 @@ class TestZarrDatatreeIO:
         with pytest.raises(zarr.errors.ContainsGroupError):
             simple_datatree.to_zarr(tmpdir)
 
+    @requires_dask
+    def test_to_zarr_compute_false(self, tmpdir, simple_datatree):
+        import dask.array as da
+
+        filepath = tmpdir / "test.zarr"
+        original_dt = simple_datatree.chunk()
+        original_dt.to_zarr(filepath, compute=False)
+
+        for node in original_dt.subtree:
+            for name, variable in node.dataset.variables.items():
+                var_dir = filepath / node.path / name
+                var_files = var_dir.listdir()
+                assert var_dir / ".zarray" in var_files
+                assert var_dir / ".zattrs" in var_files
+                if isinstance(variable.data, da.Array):
+                    assert var_dir / "0" not in var_files
+                else:
+                    assert var_dir / "0" in var_files
+
     def test_to_zarr_inherited_coords(self, tmpdir):
         original_dt = DataTree.from_dict(
             {
diff --git a/xarray/tests/test_cftime_offsets.py b/xarray/tests/test_cftime_offsets.py
index 4db73048..086e187d 100644
--- a/xarray/tests/test_cftime_offsets.py
+++ b/xarray/tests/test_cftime_offsets.py
@@ -1670,9 +1670,9 @@ def test_new_to_legacy_freq_pd_freq_passthrough(freq, expected):
         pytest.param("-1YE", marks=requires_pandas_3),
     ),
 )
-def test_cftime_range_same_as_pandas(start, end, freq):
+def test_cftime_range_same_as_pandas(start, end, freq) -> None:
     result = date_range(start, end, freq=freq, calendar="standard", use_cftime=True)
-    result = result.to_datetimeindex()
+    result = result.to_datetimeindex(time_unit="ns")
     expected = date_range(start, end, freq=freq, use_cftime=False)
 
     np.testing.assert_array_equal(result, expected)
@@ -1694,8 +1694,8 @@ def test_cftime_range_no_freq(start, end, periods):
     when freq is not provided, but start, end and periods are.
     """
     # Generate date ranges using cftime_range
-    result = cftime_range(start=start, end=end, periods=periods)
-    result = result.to_datetimeindex()
+    cftimeindex = cftime_range(start=start, end=end, periods=periods)
+    result = cftimeindex.to_datetimeindex(time_unit="ns")
     expected = pd.date_range(start=start, end=end, periods=periods)
 
     np.testing.assert_array_equal(result, expected)
diff --git a/xarray/tests/test_cftimeindex.py b/xarray/tests/test_cftimeindex.py
index 8fc79a0c..85622d54 100644
--- a/xarray/tests/test_cftimeindex.py
+++ b/xarray/tests/test_cftimeindex.py
@@ -1219,43 +1219,59 @@ def test_strftime_of_cftime_array(calendar):
 @requires_cftime
 @pytest.mark.parametrize("calendar", _ALL_CALENDARS)
 @pytest.mark.parametrize("unsafe", [False, True])
-def test_to_datetimeindex(calendar, unsafe):
+def test_to_datetimeindex(calendar, unsafe) -> None:
     index = xr.cftime_range("2000", periods=5, calendar=calendar)
-    expected = pd.date_range("2000", periods=5, unit="us")
+    expected = pd.date_range("2000", periods=5, unit="ns")
 
     if calendar in _NON_STANDARD_CALENDARS and not unsafe:
         with pytest.warns(RuntimeWarning, match="non-standard"):
-            result = index.to_datetimeindex()
+            result = index.to_datetimeindex(time_unit="ns")
     else:
-        result = index.to_datetimeindex(unsafe=unsafe)
+        result = index.to_datetimeindex(unsafe=unsafe, time_unit="ns")
 
     assert result.equals(expected)
     np.testing.assert_array_equal(result, expected)
     assert isinstance(result, pd.DatetimeIndex)
 
 
+@requires_cftime
+def test_to_datetimeindex_future_warning() -> None:
+    index = xr.cftime_range("2000", periods=5)
+    expected = pd.date_range("2000", periods=5, unit="ns")
+    with pytest.warns(FutureWarning, match="In a future version"):
+        result = index.to_datetimeindex()
+    assert result.equals(expected)
+    assert result.dtype == expected.dtype
+
+
 @requires_cftime
 @pytest.mark.parametrize("calendar", _ALL_CALENDARS)
-def test_to_datetimeindex_out_of_range(calendar):
+def test_to_datetimeindex_out_of_range(calendar) -> None:
     index = xr.cftime_range("0001", periods=5, calendar=calendar)
-    # todo: suggestion from code review:
-    #  - still warn when converting from a non-standard calendar
-    #  to a proleptic Gregorian calendar
-    #  - also warn when converting from a Gregorian calendar
-    #  to a proleptic Gregorian calendar when dates fall before the reform
-    if calendar in _NON_STANDARD_CALENDARS:
-        with pytest.warns(RuntimeWarning, match="non-standard"):
-            index.to_datetimeindex()
+    with pytest.raises(ValueError, match="0001"):
+        index.to_datetimeindex(time_unit="ns")
+
+
+@requires_cftime
+@pytest.mark.parametrize("unsafe", [False, True])
+def test_to_datetimeindex_gregorian_pre_reform(unsafe) -> None:
+    index = xr.cftime_range("1582", periods=5, calendar="gregorian")
+    if unsafe:
+        result = index.to_datetimeindex(time_unit="us", unsafe=unsafe)
     else:
-        index.to_datetimeindex()
+        with pytest.warns(RuntimeWarning, match="reform"):
+            result = index.to_datetimeindex(time_unit="us", unsafe=unsafe)
+    expected = pd.date_range("1582", periods=5, unit="us")
+    assert result.equals(expected)
+    assert result.dtype == expected.dtype
 
 
 @requires_cftime
 @pytest.mark.parametrize("calendar", ["all_leap", "360_day"])
-def test_to_datetimeindex_feb_29(calendar):
+def test_to_datetimeindex_feb_29(calendar) -> None:
     index = xr.cftime_range("2001-02-28", periods=2, calendar=calendar)
     with pytest.raises(ValueError, match="29"):
-        index.to_datetimeindex()
+        index.to_datetimeindex(time_unit="ns")
 
 
 @pytest.mark.xfail(reason="fails on pandas main branch")
@@ -1269,12 +1285,12 @@ def test_multiindex():
 @requires_cftime
 @pytest.mark.parametrize("freq", ["3663s", "33min", "2h"])
 @pytest.mark.parametrize("method", ["floor", "ceil", "round"])
-def test_rounding_methods_against_datetimeindex(freq, method):
+def test_rounding_methods_against_datetimeindex(freq, method) -> None:
     # for now unit="us" seems good enough
-    expected = pd.date_range("2000-01-02T01:03:51", periods=10, freq="1777s", unit="us")
+    expected = pd.date_range("2000-01-02T01:03:51", periods=10, freq="1777s", unit="ns")
     expected = getattr(expected, method)(freq)
     result = xr.cftime_range("2000-01-02T01:03:51", periods=10, freq="1777s")
-    result = getattr(result, method)(freq).to_datetimeindex()
+    result = getattr(result, method)(freq).to_datetimeindex(time_unit="ns")
     assert result.equals(expected)
 
 
diff --git a/xarray/tests/test_cftimeindex_resample.py b/xarray/tests/test_cftimeindex_resample.py
index 08197010..5510eaac 100644
--- a/xarray/tests/test_cftimeindex_resample.py
+++ b/xarray/tests/test_cftimeindex_resample.py
@@ -97,7 +97,9 @@ def compare_against_pandas(
         ).mean()
     # TODO (benbovy - flexible indexes): update when CFTimeIndex is a xarray Index subclass
     result_cftimeindex["time"] = (
-        result_cftimeindex.xindexes["time"].to_pandas_index().to_datetimeindex()
+        result_cftimeindex.xindexes["time"]
+        .to_pandas_index()
+        .to_datetimeindex(time_unit="ns")
     )
     xr.testing.assert_identical(result_cftimeindex, result_datetimeindex)
 
@@ -181,7 +183,7 @@ def test_calendars(calendar: str) -> None:
     # TODO (benbovy - flexible indexes): update when CFTimeIndex is a xarray Index subclass
     new_pd_index = da_cftime.xindexes["time"].to_pandas_index()
     assert isinstance(new_pd_index, CFTimeIndex)  # shouldn't that be a pd.Index?
-    da_cftime["time"] = new_pd_index.to_datetimeindex()
+    da_cftime["time"] = new_pd_index.to_datetimeindex(time_unit="ns")
     xr.testing.assert_identical(da_cftime, da_datetime)
 
 
diff --git a/xarray/tests/test_duck_array_ops.py b/xarray/tests/test_duck_array_ops.py
index e1306964..9b8b50d9 100644
--- a/xarray/tests/test_duck_array_ops.py
+++ b/xarray/tests/test_duck_array_ops.py
@@ -9,6 +9,7 @@ import pytest
 from numpy import array, nan
 
 from xarray import DataArray, Dataset, cftime_range, concat
+from xarray.coding.times import _NS_PER_TIME_DELTA
 from xarray.core import dtypes, duck_array_ops
 from xarray.core.duck_array_ops import (
     array_notnull_equiv,
@@ -28,6 +29,7 @@ from xarray.core.duck_array_ops import (
     where,
 )
 from xarray.core.extension_array import PandasExtensionArray
+from xarray.core.types import NPDatetimeUnitOptions, PDDatetimeUnitOptions
 from xarray.namedarray.pycompat import array_type
 from xarray.testing import assert_allclose, assert_equal, assert_identical
 from xarray.tests import (
@@ -411,10 +413,11 @@ def assert_dask_array(da, dask):
 @arm_xfail
 @pytest.mark.filterwarnings("ignore:All-NaN .* encountered:RuntimeWarning")
 @pytest.mark.parametrize("dask", [False, True] if has_dask else [False])
-def test_datetime_mean(dask: bool) -> None:
+def test_datetime_mean(dask: bool, time_unit: PDDatetimeUnitOptions) -> None:
     # Note: only testing numpy, as dask is broken upstream
+    dtype = f"M8[{time_unit}]"
     da = DataArray(
-        np.array(["2010-01-01", "NaT", "2010-01-03", "NaT", "NaT"], dtype="M8[ns]"),
+        np.array(["2010-01-01", "NaT", "2010-01-03", "NaT", "NaT"], dtype=dtype),
         dims=["time"],
     )
     if dask:
@@ -846,11 +849,11 @@ def test_multiple_dims(dtype, dask, skipna, func):
 
 
 @pytest.mark.parametrize("dask", [True, False])
-def test_datetime_to_numeric_datetime64(dask):
+def test_datetime_to_numeric_datetime64(dask, time_unit: PDDatetimeUnitOptions):
     if dask and not has_dask:
         pytest.skip("requires dask")
 
-    times = pd.date_range("2000", periods=5, freq="7D").values
+    times = pd.date_range("2000", periods=5, freq="7D").as_unit(time_unit).values
     if dask:
         import dask.array
 
@@ -874,8 +877,8 @@ def test_datetime_to_numeric_datetime64(dask):
         result = duck_array_ops.datetime_to_numeric(
             times, datetime_unit="h", dtype=dtype
         )
-    expected = 24 * np.arange(0, 35, 7).astype(dtype)
-    np.testing.assert_array_equal(result, expected)
+    expected2 = 24 * np.arange(0, 35, 7).astype(dtype)
+    np.testing.assert_array_equal(result, expected2)
 
 
 @requires_cftime
@@ -923,15 +926,18 @@ def test_datetime_to_numeric_cftime(dask):
 
 
 @requires_cftime
-def test_datetime_to_numeric_potential_overflow():
+def test_datetime_to_numeric_potential_overflow(time_unit: PDDatetimeUnitOptions):
     import cftime
 
-    times = pd.date_range("2000", periods=5, freq="7D").values.astype("datetime64[us]")
+    if time_unit == "ns":
+        pytest.skip("out-of-bounds datetime64 overflow")
+    dtype = f"M8[{time_unit}]"
+    times = pd.date_range("2000", periods=5, freq="7D").values.astype(dtype)
     cftimes = cftime_range(
         "2000", periods=5, freq="7D", calendar="proleptic_gregorian"
     ).values
 
-    offset = np.datetime64("0001-01-01")
+    offset = np.datetime64("0001-01-01", time_unit)
     cfoffset = cftime.DatetimeProlepticGregorian(1, 1, 1)
 
     result = duck_array_ops.datetime_to_numeric(
@@ -957,24 +963,33 @@ def test_py_timedelta_to_float():
     assert py_timedelta_to_float(dt.timedelta(days=1e6), "D") == 1e6
 
 
-@pytest.mark.parametrize(
-    "td, expected",
-    ([np.timedelta64(1, "D"), 86400 * 1e9], [np.timedelta64(1, "ns"), 1.0]),
-)
-def test_np_timedelta64_to_float(td, expected):
-    out = np_timedelta64_to_float(td, datetime_unit="ns")
+@pytest.mark.parametrize("np_dt_unit", ["D", "h", "m", "s", "ms", "us", "ns"])
+def test_np_timedelta64_to_float(
+    np_dt_unit: NPDatetimeUnitOptions, time_unit: PDDatetimeUnitOptions
+):
+    # tests any combination of source np.timedelta64 (NPDatetimeUnitOptions) with
+    # np_timedelta_to_float with dedicated target unit (PDDatetimeUnitOptions)
+    td = np.timedelta64(1, np_dt_unit)
+    expected = _NS_PER_TIME_DELTA[np_dt_unit] / _NS_PER_TIME_DELTA[time_unit]
+
+    out = np_timedelta64_to_float(td, datetime_unit=time_unit)
     np.testing.assert_allclose(out, expected)
     assert isinstance(out, float)
 
-    out = np_timedelta64_to_float(np.atleast_1d(td), datetime_unit="ns")
+    out = np_timedelta64_to_float(np.atleast_1d(td), datetime_unit=time_unit)
     np.testing.assert_allclose(out, expected)
 
 
-@pytest.mark.parametrize(
-    "td, expected", ([pd.Timedelta(1, "D"), 86400 * 1e9], [pd.Timedelta(1, "ns"), 1.0])
-)
-def test_pd_timedelta_to_float(td, expected):
-    out = pd_timedelta_to_float(td, datetime_unit="ns")
+@pytest.mark.parametrize("np_dt_unit", ["D", "h", "m", "s", "ms", "us", "ns"])
+def test_pd_timedelta_to_float(
+    np_dt_unit: NPDatetimeUnitOptions, time_unit: PDDatetimeUnitOptions
+):
+    # tests any combination of source pd.Timedelta (NPDatetimeUnitOptions) with
+    # np_timedelta_to_float with dedicated target unit (PDDatetimeUnitOptions)
+    td = pd.Timedelta(1, np_dt_unit)
+    expected = _NS_PER_TIME_DELTA[np_dt_unit] / _NS_PER_TIME_DELTA[time_unit]
+
+    out = pd_timedelta_to_float(td, datetime_unit=time_unit)
     np.testing.assert_allclose(out, expected)
     assert isinstance(out, float)
 
@@ -982,10 +997,11 @@ def test_pd_timedelta_to_float(td, expected):
 @pytest.mark.parametrize(
     "td", [dt.timedelta(days=1), np.timedelta64(1, "D"), pd.Timedelta(1, "D"), "1 day"]
 )
-def test_timedelta_to_numeric(td):
+def test_timedelta_to_numeric(td, time_unit: PDDatetimeUnitOptions):
     # Scalar input
-    out = timedelta_to_numeric(td, "ns")
-    np.testing.assert_allclose(out, 86400 * 1e9)
+    out = timedelta_to_numeric(td, time_unit)
+    expected = _NS_PER_TIME_DELTA["D"] / _NS_PER_TIME_DELTA[time_unit]
+    np.testing.assert_allclose(out, expected)
     assert isinstance(out, float)
 
 
diff --git a/xarray/tests/test_groupby.py b/xarray/tests/test_groupby.py
index 3c7f46f5..2287130e 100644
--- a/xarray/tests/test_groupby.py
+++ b/xarray/tests/test_groupby.py
@@ -1618,6 +1618,8 @@ 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]]],
@@ -1868,7 +1870,7 @@ class TestDataArrayResample:
         def resample_as_pandas(array, *args, **kwargs):
             array_ = array.copy(deep=True)
             if use_cftime:
-                array_["time"] = times.to_datetimeindex()
+                array_["time"] = times.to_datetimeindex(time_unit="ns")
             result = DataArray.from_series(
                 array_.to_series().resample(*args, **kwargs).mean()
             )
@@ -2321,7 +2323,7 @@ class TestDatasetResample:
         def resample_as_pandas(ds, *args, **kwargs):
             ds_ = ds.copy(deep=True)
             if use_cftime:
-                ds_["time"] = times.to_datetimeindex()
+                ds_["time"] = times.to_datetimeindex(time_unit="ns")
             result = Dataset.from_dataframe(
                 ds_.to_dataframe().resample(*args, **kwargs).mean()
             )
diff --git a/xarray/tests/test_missing.py b/xarray/tests/test_missing.py
index eb21cca0..d3eaa70f 100644
--- a/xarray/tests/test_missing.py
+++ b/xarray/tests/test_missing.py
@@ -606,10 +606,12 @@ def test_get_clean_interp_index_cf_calendar(cf_da, calendar):
 @requires_cftime
 @pytest.mark.parametrize("calendar", ["gregorian", "proleptic_gregorian"])
 @pytest.mark.parametrize("freq", ["1D", "1ME", "1YE"])
-def test_get_clean_interp_index_dt(cf_da, calendar, freq):
+def test_get_clean_interp_index_dt(cf_da, calendar, freq) -> None:
     """In the gregorian case, the index should be proportional to normal datetimes."""
     g = cf_da(calendar, freq=freq)
-    g["stime"] = xr.Variable(data=g.time.to_index().to_datetimeindex(), dims=("time",))
+    g["stime"] = xr.Variable(
+        data=g.time.to_index().to_datetimeindex(time_unit="ns"), dims=("time",)
+    )
 
     gi = get_clean_interp_index(g, "time")
     si = get_clean_interp_index(g, "time", use_coordinate="stime")