Minor fix to combine_by_coords to allow for the combination of CFTimeIndexes separated by large time intervals

doc/whats-new.rst

@@ -34,6 +34,9 @@ New Features
 
 Bug fixes
 ~~~~~~~~~
+- Fix :py:meth:`xarray.combine_by_coords` when combining cftime coordinates
+  which span long time intervals (:issue:`3535`).  By `Spencer Clark
+  <https://github.com/spencerkclark>`_.
 - Fix plotting with transposed 2D non-dimensional coordinates. (:issue:`3138`, :pull:`3441`)
   By `Deepak Cherian <https://github.com/dcherian>`_.
 - Fix issue with Dask-backed datasets raising a ``KeyError`` on some computations involving ``map_blocks`` (:pull:`3598`)

xarray/core/combine.py

@@ -88,7 +88,7 @@ def _infer_concat_order_from_coords(datasets):
                 # with the same value have the same coord values throughout.
                 if any(index.size == 0 for index in indexes):
                     raise ValueError("Cannot handle size zero dimensions")
-                first_items = pd.Index([index.take([0]) for index in indexes])
+                first_items = pd.Index([index[0] for index in indexes])
 
                 # Sort datasets along dim
                 # We want rank but with identical elements given identical

xarray/tests/test_combine.py

@@ -22,7 +22,7 @@
     _new_tile_id,
 )
 
-from . import assert_equal, assert_identical, raises_regex
+from . import assert_equal, assert_identical, raises_regex, requires_cftime
 from .test_dataset import create_test_data
 
 
@@ -877,3 +877,25 @@ def test_auto_combine_without_coords(self):
         objs = [Dataset({"foo": ("x", [0])}), Dataset({"foo": ("x", [1])})]
         with pytest.warns(FutureWarning, match="supplied do not have global"):
             auto_combine(objs)
+
+
+@requires_cftime
+def test_combine_by_coords_distant_cftime_dates():
+    # Regression test for https://github.com/pydata/xarray/issues/3535
+    import cftime
+
+    time_1 = [cftime.DatetimeGregorian(4500, 12, 31)]
+    time_2 = [cftime.DatetimeGregorian(4600, 12, 31)]
+    time_3 = [cftime.DatetimeGregorian(5100, 12, 31)]
+
+    da_1 = DataArray([0], dims=["time"], coords=[time_1], name="a").to_dataset()
+    da_2 = DataArray([1], dims=["time"], coords=[time_2], name="a").to_dataset()
+    da_3 = DataArray([2], dims=["time"], coords=[time_3], name="a").to_dataset()
+
+    result = combine_by_coords([da_1, da_2, da_3])
+
+    expected_time = np.concatenate([time_1, time_2, time_3])
+    expected = DataArray(
+        [0, 1, 2], dims=["time"], coords=[expected_time], name="a"
+    ).to_dataset()
+    assert_identical(result, expected)