document limitations of cftime arithmetics

doc/user-guide/weather-climate.rst

@@ -279,3 +279,90 @@ For data indexed by a :py:class:`~xarray.CFTimeIndex` xarray currently supports:
 .. _precision range: https://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html#timestamp-limitations
 .. _ISO 8601 standard: https://en.wikipedia.org/wiki/ISO_8601
 .. _partial datetime string indexing: https://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html#partial-string-indexing
+
+.. _cftime_arithmetic_limitations:
+
+Arithmetic limitations with ``cftime`` objects
+----------------------------------------------
+
+A current limitation when working with non-standard calendars and :py:class:`cftime.datetime`
+objects is that they support arithmetic with :py:class:`datetime.timedelta`, but **not** with :py:class:`numpy.timedelta64`.
+
+This means that certain xarray operations (such as :py:meth:`~xarray.DataArray.diff`)
+may produce ``timedelta64`` results that cannot be directly combined with ``cftime`` coordinates.
+
+For example, let's define a time axis using ``cftime`` objects:
+
+.. jupyter-execute::
+
+    import xarray as xr
+    import numpy as np
+    import pandas as pd
+    import cftime
+
+    time = xr.DataArray(
+        xr.date_range("2000", periods=3, freq="MS", use_cftime=True),
+        dims="time",
+    )
+
+If you want to compute, e.g., midpoints in the time intervals, this will not work:
+
+.. code-block:: python
+
+    # Attempt to compute midpoints
+    time[:-1] + 0.5 * time.diff("time")
+
+and result in an error like this:
+
+.. code-block:: none
+
+   UFuncTypeError: ufunc 'add' cannot use operands with types dtype('O') and dtype('<m8[ns]')
+
+This is because :py:meth:`~xarray.DataArray.diff` returns ``timedelta64``, which is not
+compatible with ``cftime.datetime``. These limitations stem from the ``cftime`` library itself;
+arithmetic between ``cftime.datetime`` and ``numpy.timedelta64`` is not implemented.
+
+**Workarounds**
+
+One possible workaround is to use :py:meth:`numpy.diff` on the underlying values,
+which returns ``datetime.timedelta`` objects that are compatible:
+
+.. jupyter-execute::
+
+    time[:-1] + 0.5 * np.diff(time.values)
+
+or you can convert ``timedelta64`` values to Python ``timedelta`` objects explicitly,
+for example via :py:meth:`pandas.to_timedelta`:
+
+.. jupyter-execute::
+
+    td = pd.to_timedelta(time.diff("time").values).to_pytimedelta()
+    time[:-1] + 0.5 * td
+
+You could also consider a more verbose, but also more robust, workaround that adds the two
+DataArrays via :py:meth:`~xarray.apply_ufunc`:
+
+.. jupyter-execute::
+
+    def add_cftime_and_timedelta64(cftime_da, timedelta64_da, join="inner"):
+        def func(cftime_array, timedelta64_array):
+            shape = cftime_array.shape
+            cftime_array = cftime_array.ravel()
+            timedelta64_array = timedelta64_array.ravel()
+            timedelta_array = pd.to_timedelta(timedelta64_array).to_pytimedelta()
+            return (cftime_array + timedelta_array).reshape(shape)
+
+        return xr.apply_ufunc(
+            func,
+            cftime_da,
+            timedelta64_da,
+            dask="parallelized",
+            output_dtypes=[cftime_da.dtype],
+            join=join
+        )
+
+
+    add_cftime_and_timedelta64(time, 0.5 * time.diff("time"))
+
+This function has the advantage that it preserves the coordinate alignment, as well as multi-dimensional
+and dask compatibility features of xarray.

doc/whats-new.rst

@@ -53,6 +53,8 @@ Bug Fixes
 Documentation
 ~~~~~~~~~~~~~
 
+- Added section on the `limitations of cftime arithmetic <https://docs.xarray.dev/en/stable/user-guide/weather-climate.html#arithmetic-limitations-with-non-standard-calendars>`_ (:pull:`10653`).
+  By `Lars Buntemeyer <https://github.com/larsbuntemeyer>`_
 
 Internal Changes
 ~~~~~~~~~~~~~~~~