Averaging timestamps with non-nanosecond precision

[sfinkens]

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

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

• Minimal example — the example is as focused as reasonably possible to demonstrate the underlying issue in xarray.
• Complete example — the example is self-contained, including all data and the text of any traceback.
• Verifiable example — the example copy & pastes into an IPython prompt or Binder notebook, returning the result.
• New issue — a search of GitHub Issues suggests this is not a duplicate.
• Recent environment — the issue occurs with the latest version of xarray and its dependencies.

Relevant log output

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

Anything else we need to know?

No response

Environment

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

[kmuehlbauer]

Thanks for the report @sfinkens.

This really indicates it's interpreting the int64 values as nanosecond values. The 1800 nanoseconds should be seconds in this case, which would be the expected half hour

[kmuehlbauer]

How did I miss this? I've thought I've grep'ed through the codebase with nanosecond 🤷.

xarray/xarray/core/duck_array_ops.py

Lines 715 to 732 in 609412d

	if dtypes.is_datetime_like(array.dtype):
	offset = _datetime_nanmin(array)
	# xarray always uses np.datetime64[ns] for np.datetime64 data
	dtype = "timedelta64[ns]"
	return (
	_mean(
	datetime_to_numeric(array, offset), axis=axis, skipna=skipna, **kwargs
	).astype(dtype)
	+ offset
	)
	elif _contains_cftime_datetimes(array):
	offset = min(array)
	timedeltas = datetime_to_numeric(array, offset, datetime_unit="us")
	mean_timedeltas = _mean(timedeltas, axis=axis, skipna=skipna, **kwargs)
	return _to_pytimedelta(mean_timedeltas, unit="us") + offset
	else:
	return _mean(array, axis=axis, skipna=skipna, **kwargs)

[kmuehlbauer]

@sfinkens Would you mind checking your workflows with #9977?

There might still be some issues lurking in some parts of the computation when using non-nanosecond datetimes/timedeltas. See also issue #9937.