Repository: pydata/xarray

Issue:
more upstream-dev cftime failures
https://dev.azure.com/xarray/xarray/_build/results?buildId=2116&view=logs&jobId=2280efed-fda1-53bd-9213-1fa8ec9b4fa8&j=2280efed-fda1-53bd-9213-1fa8ec9b4fa8&t=175181ee-1928-5a6b-f537-168f7a8b7c2d


46 failed tests but they all seem to result from the same TypeError

```
=================================== FAILURES ===================================
______________ test_sel_date_scalar_nearest[365_day-sel_kwargs0] _______________

da = <xarray.DataArray (time: 4)>
             0002-02-01 00:00:00],
            dtype='object')
sel_kwargs = {'method': 'nearest'}

    @requires_cftime
    @pytest.mark.parametrize(
        "sel_kwargs",
        [{"method": "nearest"}, {"method": "nearest", "tolerance": timedelta(days=70)}],
    )
    def test_sel_date_scalar_nearest(da, date_type, index, sel_kwargs):
        expected = xr.DataArray(2).assign_coords(time=index[1])
>       result = da.sel(time=date_type(1, 4, 1), **sel_kwargs)

xarray/tests/test_cftimeindex.py:460: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
xarray/core/dataarray.py:1056: in sel
    ds = self._to_temp_dataset().sel(
xarray/core/dataset.py:2056: in sel
    pos_indexers, new_indexes = remap_label_indexers(
xarray/core/coordinates.py:391: in remap_label_indexers
    pos_indexers, new_indexes = indexing.remap_label_indexers(
xarray/core/indexing.py:270: in remap_label_indexers
    idxr, new_idx = convert_label_indexer(index, label, dim, method, tolerance)
xarray/core/indexing.py:189: in convert_label_indexer
    indexer = index.get_loc(
xarray/coding/cftimeindex.py:334: in get_loc
    return pd.Index.get_loc(self, key, method=method, tolerance=tolerance)
/usr/share/miniconda/envs/xarray-tests/lib/python3.8/site-packages/pandas/core/indexes/base.py:2899: in get_loc
    indexer = self.get_indexer([key], method=method, tolerance=tolerance)
/usr/share/miniconda/envs/xarray-tests/lib/python3.8/site-packages/pandas/core/indexes/base.py:2992: in get_indexer
    indexer = self._get_nearest_indexer(target, limit, tolerance)
/usr/share/miniconda/envs/xarray-tests/lib/python3.8/site-packages/pandas/core/indexes/base.py:3076: in _get_nearest_indexer
    left_distances = np.abs(self[left_indexer] - target)
xarray/coding/cftimeindex.py:444: in __sub__
    return CFTimeIndex(np.array(self) - other)
xarray/coding/cftimeindex.py:248: in __new__
    assert_all_valid_date_type(data)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

data = TimedeltaIndex(['-59 days'], dtype='timedelta64[ns]', freq=None)

    def assert_all_valid_date_type(data):
        import cftime
    
        if len(data) > 0:
            sample = data[0]
            date_type = type(sample)
            if not isinstance(sample, cftime.datetime):
>               raise TypeError(
                    "CFTimeIndex requires cftime.datetime "
                    "objects. Got object of {}.".format(date_type)
                )
E               TypeError: CFTimeIndex requires cftime.datetime objects. Got object of <class 'pandas._libs.tslibs.timedeltas.Timedelta'>.

xarray/coding/cftimeindex.py:206: TypeError
```


Discussion:
spencerkclark: This is a little bit trickier than I originally anticipated.  For indexing with dates very distant from the dates in the index, I'm still running into an issue at this step in `pandas.Index._get_nearest_indexer`:
```
left_distances = np.abs(self[left_indexer] - target)
```
Consider the following the example:
```
In [1]: import numpy as np; import pandas as pd; import xarray as xr

In [2]: import cftime

In [3]: times = xr.cftime_range('0001', periods=4)

In [4]: da = xr.DataArray(range(4), coords=[('time', times)])

In [5]: da.sel(time=cftime.DatetimeGregorian(2000, 1, 1), method='nearest')
```

In this example `self[left_indexer]` is `CFTimeIndex([0001-01-04 00:00:00], dtype='object', name='time')`, while `target` is `Index([2000-01-01 00:00:00], dtype='object')`.  The distance between these dates is greater than 292 years, so it cannot be represented in a `TimedeltaIndex`.  

One could argue that we could fall back on using a generic `Index` of dtype object to store the `datetime.timedelta` object produced in this case:
```
In [6]: left_index = xr.CFTimeIndex([cftime.DatetimeGregorian(1, 1, 4)])

In [7]: target = pd.Index([cftime.DatetimeGregorian(2000, 1, 1)])

In [8]: difference = left_index - target

In [9]: difference
Out[9]: Index([-730118 days, 0:00:00], dtype='object')
```
A problem occurs though when we try to take the absolute value of this index.  Pandas (I think reasonably so) tries to detect the datatype of the result and construct a new index.  In doing so it tries to create a `TimedeltaIndex`, but cannot, because the `timedelta` inside is too large:

<details><summary>Result of np.abs(difference)</summary>
<p>

```
In [10]: np.abs(difference)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
~/Software/pandas/pandas/_libs/tslibs/timedeltas.pyx in pandas._libs.tslibs.timedeltas.array_to_timedelta64()

TypeError: Expected unicode, got datetime.timedelta

During handling of the above exception, another exception occurred:

OverflowError                             Traceback (most recent call last)
<ipython-input-17-95776624315c> in <module>
----> 1 np.abs(difference)

~/Software/pandas/pandas/core/indexes/base.py in __array_wrap__(self, result, context)
    628
    629         attrs = self._get_attributes_dict()
--> 630         return Index(result, **attrs)
    631
    632     @cache_readonly

~/Software/pandas/pandas/core/indexes/base.py in __new__(cls, data, dtype, copy, name, tupleize_cols, **kwargs)
    412
    413             if dtype is None:
--> 414                 new_data, new_dtype = _maybe_cast_data_without_dtype(subarr)
    415                 if new_dtype is not None:
    416                     return cls(

~/Software/pandas/pandas/core/indexes/base.py in _maybe_cast_data_without_dtype(subarr)
   5711
   5712         elif inferred.startswith("timedelta"):
-> 5713             data = TimedeltaArray._from_sequence(subarr, copy=False)
   5714             return data, data.dtype
   5715         elif inferred == "period":

~/Software/pandas/pandas/core/arrays/timedeltas.py in _from_sequence(cls, data, dtype, copy, freq, unit)
    212         freq, freq_infer = dtl.maybe_infer_freq(freq)
    213
--> 214         data, inferred_freq = sequence_to_td64ns(data, copy=copy, unit=unit)
    215         freq, freq_infer = dtl.validate_inferred_freq(freq, inferred_freq, freq_infer)
    216

~/Software/pandas/pandas/core/arrays/timedeltas.py in sequence_to_td64ns(data, copy, unit, errors)
    938     if is_object_dtype(data.dtype) or is_string_dtype(data.dtype):
    939         # no need to make a copy, need to convert if string-dtyped
--> 940         data = objects_to_td64ns(data, unit=unit, errors=errors)
    941         copy = False
    942

~/Software/pandas/pandas/core/arrays/timedeltas.py in objects_to_td64ns(data, unit, errors)
   1047     values = np.array(data, dtype=np.object_, copy=False)
   1048
-> 1049     result = array_to_timedelta64(values, unit=unit, errors=errors)
   1050     return result.view("timedelta64[ns]")
   1051

~/Software/pandas/pandas/_libs/tslibs/timedeltas.pyx in pandas._libs.tslibs.timedeltas.array_to_timedelta64()

~/Software/pandas/pandas/_libs/tslibs/timedeltas.pyx in pandas._libs.tslibs.timedeltas.convert_to_timedelta64()

~/Software/pandas/pandas/_libs/tslibs/timedeltas.pyx in pandas._libs.tslibs.timedeltas.delta_to_nanoseconds()

OverflowError: Python int too large to convert to C long
```

</p>
</details>


So I'm a little bit back to the drawing board regarding a solution for this.  Part of me is tempted to write an overriding version of `_get_nearest_indexer` for `CFTimeIndex` that basically restores the old behavior; I'm not sure how dangerous that would be to rely on though.
spencerkclark: I added some updates to this PR this morning that in principle would solve the indexing with `method="nearest"` from within xarray.  Unfortunately though, due to the issue I described in https://github.com/pydata/xarray/pull/3764#issuecomment-586597512, I was not able to come up with a solution that did not require overriding the pandas implementation of `_get_nearest_indexer`.  If this seems unacceptable, maybe we can think harder about how we might address this upstream instead (e.g. I think special-casing the change made in https://github.com/pandas-dev/pandas/pull/31511 to just DatetimeIndexes, and preserving the old behavior for everything else, could be sufficient).

In addition, for the time being I xfailed `test_indexing_in_series_getitem`, as I think there is agreement that that would be best addressed upstream, https://github.com/pydata/xarray/issues/3751#issuecomment-587572443.

Finally in the process of doing this, I cleaned up the implementation of `CFTimeIndex.__sub__`, and added some more test coverage; hopefully now it's a little clearer what cases it's supposed to work for and what cases it is not.
dcherian: ping @shoyer for input.
max-sixty: Re master failing, what do we think about `xfail`-ing those tests on master and then adding a PR which undoes the xfail and attempts fixes the issue? 

I'm cognizant that every PR now fails, which less friendly to contributors. The project home page also states ![image](https://user-images.githubusercontent.com/5635139/75507496-429dad00-59af-11ea-84d5-ab7cc6f70d12.png).

max-sixty: Re the specific issue, that's a tough one. It could be worth floating on pandas and seeing if they have thoughts...
shoyer: +1 for marking as xfail. There’s no point in keeping the master branch
failing if we know what the issue is.

On Thu, Feb 27, 2020 at 7:22 PM Maximilian Roos <notifications@github.com>
wrote:

> Re the specific issue, that's a tough one. It could be worth floating on
> pandas and seeing if they have thoughts...
>
> —
> You are receiving this because you were mentioned.
> Reply to this email directly, view it on GitHub
> <https://github.com/pydata/xarray/pull/3764?email_source=notifications&email_token=AAJJFVXTH5QBZHOH47UZIQDRFB7NXA5CNFSM4KTAWZ7KYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOENG2V6Y#issuecomment-592292603>,
> or unsubscribe
> <https://github.com/notifications/unsubscribe-auth/AAJJFVQKSGKMKA35V4JDV4DRFB7NXANCNFSM4KTAWZ7A>
> .
>

pep8speaks: Hello @spencerkclark! Thanks for updating this PR. We checked the lines you've touched for [PEP 8](https://www.python.org/dev/peps/pep-0008) issues, and found:

* In the file [`xarray/coding/cftimeindex.py`](https://github.com/pydata/xarray/blob/2424de0dfb2c57f009766e1cfa2099e5a21b6f6d/xarray/coding/cftimeindex.py):

> [Line 729:27](https://github.com/pydata/xarray/blob/2424de0dfb2c57f009766e1cfa2099e5a21b6f6d/xarray/coding/cftimeindex.py#L729): [W292](https://duckduckgo.com/?q=pep8%20W292) no newline at end of file





Patch:
diff --git a/ci/requirements/py36-min-nep18.yml b/ci/requirements/py36-min-nep18.yml
index c10fdf67..a5eded49 100644
--- a/ci/requirements/py36-min-nep18.yml
+++ b/ci/requirements/py36-min-nep18.yml
@@ -11,7 +11,7 @@ dependencies:
   - msgpack-python=0.6  # remove once distributed is bumped. distributed GH3491
   - numpy=1.17
   - pandas=0.25
-  - pint=0.9  # Actually not enough as it doesn't implement __array_function__yet!
+  - pint=0.11
   - pip
   - pytest
   - pytest-cov
diff --git a/doc/contributing.rst b/doc/contributing.rst
index eb31db24..f581bcd9 100644
--- a/doc/contributing.rst
+++ b/doc/contributing.rst
@@ -51,8 +51,8 @@ Bug reports must:
    <http://github.github.com/github-flavored-markdown/>`_::
 
       ```python
-      >>> from xarray import Dataset
-      >>> df = Dataset(...)
+      >>> import xarray as xr
+      >>> df = xr.Dataset(...)
       ...
       ```
 
@@ -378,8 +378,8 @@ and then running::
 
    pre-commit install
 
-from the root of the xarray repository. You can skip the pre-commit checks with
-``git commit --no-verify``.
+from the root of the xarray repository. You can skip the pre-commit checks
+with ``git commit --no-verify``.
 
 
 Backwards Compatibility
diff --git a/doc/installing.rst b/doc/installing.rst
index dfc2841a..a25bf65e 100644
--- a/doc/installing.rst
+++ b/doc/installing.rst
@@ -11,6 +11,8 @@ Required dependencies
 - `numpy <http://www.numpy.org/>`__ (1.15 or later)
 - `pandas <http://pandas.pydata.org/>`__ (0.25 or later)
 
+.. _optional-dependencies:
+
 Optional dependencies
 ---------------------
 
@@ -24,7 +26,7 @@ For netCDF and IO
 - `h5netcdf <https://github.com/shoyer/h5netcdf>`__: an alternative library for
   reading and writing netCDF4 files that does not use the netCDF-C libraries
 - `pynio <https://www.pyngl.ucar.edu/Nio.shtml>`__: for reading GRIB and other
-  geoscience specific file formats
+  geoscience specific file formats. Note that pynio is not available for Windows.
 - `zarr <http://zarr.readthedocs.io/>`__: for chunked, compressed, N-dimensional arrays.
 - `cftime <https://unidata.github.io/cftime>`__: recommended if you
   want to encode/decode datetimes for non-standard calendars or dates before
@@ -121,16 +123,15 @@ xarray itself is a pure Python package, but its dependencies are not. The
 easiest way to get everything installed is to use conda_. To install xarray
 with its recommended dependencies using the conda command line tool::
 
-    $ conda install xarray dask netCDF4 bottleneck
+    $ conda install -c conda-forge xarray dask netCDF4 bottleneck
 
 .. _conda: http://conda.io/
 
-We recommend using the community maintained `conda-forge <https://conda-forge.github.io/>`__ channel if you need difficult\-to\-build dependencies such as cartopy, pynio or PseudoNetCDF::
-
-    $ conda install -c conda-forge xarray cartopy pynio pseudonetcdf
+If you require other :ref:`optional-dependencies` add them to the line above.
 
-New releases may also appear in conda-forge before being updated in the default
-channel.
+We recommend using the community maintained `conda-forge <https://conda-forge.github.io/>`__ channel,
+as some of the dependencies are difficult to build. New releases may also appear in conda-forge before
+being updated in the default channel.
 
 If you don't use conda, be sure you have the required dependencies (numpy and
 pandas) installed first. Then, install xarray with pip::
diff --git a/doc/related-projects.rst b/doc/related-projects.rst
index edee80b7..57b8da0c 100644
--- a/doc/related-projects.rst
+++ b/doc/related-projects.rst
@@ -63,6 +63,7 @@ Extend xarray capabilities
 - `hypothesis-gufunc <https://hypothesis-gufunc.readthedocs.io/en/latest/>`_: Extension to hypothesis. Makes it easy to write unit tests with xarray objects as input.
 - `nxarray <https://github.com/nxarray/nxarray>`_: NeXus input/output capability for xarray.
 - `xarray_extras <https://github.com/crusaderky/xarray_extras>`_: Advanced algorithms for xarray objects (e.g. integrations/interpolations).
+- `xpublish <https://xpublish.readthedocs.io/>`_: Publish Xarray Datasets via a Zarr compatible REST API.
 - `xrft <https://github.com/rabernat/xrft>`_: Fourier transforms for xarray data.
 - `xr-scipy <https://xr-scipy.readthedocs.io>`_: A lightweight scipy wrapper for xarray.
 - `X-regression <https://github.com/kuchaale/X-regression>`_: Multiple linear regression from Statsmodels library coupled with Xarray library.
diff --git a/doc/whats-new.rst b/doc/whats-new.rst
index bc0e5092..80309dc4 100644
--- a/doc/whats-new.rst
+++ b/doc/whats-new.rst
@@ -32,7 +32,7 @@ New Features
 - Support new h5netcdf backend keyword `phony_dims` (available from h5netcdf
   v0.8.0 for :py:class:`~xarray.backends.H5NetCDFStore`.
   By `Kai Mühlbauer <https://github.com/kmuehlbauer>`_.
-- implement pint support. (:issue:`3594`, :pull:`3706`)
+- Support unit aware arrays with pint. (:issue:`3594`, :pull:`3706`, :pull:`3611`)
   By `Justus Magin <https://github.com/keewis>`_.
 - :py:meth:`Dataset.groupby` and :py:meth:`DataArray.groupby` now raise a 
   `TypeError` on multiple string arguments. Receiving multiple string arguments
@@ -81,6 +81,9 @@ Documentation
 - Fix documentation of :py:class:`DataArray` removing the deprecated mention
   that when omitted, `dims` are inferred from a `coords`-dict. (:pull:`3821`)
   By `Sander van Rijn <https://github.com/sjvrijn>`_.
+- Update the installation instructions: only explicitly list recommended dependencies
+  (:issue:`3756`).
+  By `Mathias Hauser <https://github.com/mathause>`_.
 
 Internal Changes
 ~~~~~~~~~~~~~~~~
@@ -93,6 +96,12 @@ Internal Changes
 - Changed test_open_mfdataset_list_attr to only run with dask installed
   (:issue:`3777`, :pull:`3780`).
   By `Bruno Pagani <https://github.com/ArchangeGabriel>`_.
+- Preserved the ability to index with ``method="nearest"`` with a
+  :py:class:`CFTimeIndex` with pandas versions greater than 1.0.1
+  (:issue:`3751`). By `Spencer Clark <https://github.com/spencerkclark>`_.
+- Greater flexibility and improved test coverage of subtracting various types
+  of objects from a :py:class:`CFTimeIndex`. By `Spencer Clark
+  <https://github.com/spencerkclark>`_.
 - Updated Azure CI MacOS image, given pending removal.
   By `Maximilian Roos <https://github.com/max-sixty>`_
 - Removed xfails for scipy 1.0.1 for tests that append to netCDF files (:pull:`3805`).
diff --git a/xarray/coding/cftimeindex.py b/xarray/coding/cftimeindex.py
index 99f90430..1ea5d3a7 100644
--- a/xarray/coding/cftimeindex.py
+++ b/xarray/coding/cftimeindex.py
@@ -49,6 +49,7 @@ import pandas as pd
 
 from xarray.core.utils import is_scalar
 
+from ..core.common import _contains_cftime_datetimes
 from .times import _STANDARD_CALENDARS, cftime_to_nptime, infer_calendar_name
 
 
@@ -326,6 +327,32 @@ class CFTimeIndex(pd.Index):
             raise KeyError(key)
         return loc
 
+    def _get_nearest_indexer(self, target, limit, tolerance):
+        """Adapted from pandas.Index._get_nearest_indexer"""
+        left_indexer = self.get_indexer(target, "pad", limit=limit)
+        right_indexer = self.get_indexer(target, "backfill", limit=limit)
+        left_distances = abs(self.values[left_indexer] - target.values)
+        right_distances = abs(self.values[right_indexer] - target.values)
+
+        if self.is_monotonic_increasing:
+            condition = (left_distances < right_distances) | (right_indexer == -1)
+        else:
+            condition = (left_distances <= right_distances) | (right_indexer == -1)
+        indexer = np.where(condition, left_indexer, right_indexer)
+
+        if tolerance is not None:
+            indexer = self._filter_indexer_tolerance(target, indexer, tolerance)
+        return indexer
+
+    def _filter_indexer_tolerance(self, target, indexer, tolerance):
+        """Adapted from pandas.Index._filter_indexer_tolerance"""
+        if isinstance(target, pd.Index):
+            distance = abs(self.values[indexer] - target.values)
+        else:
+            distance = abs(self.values[indexer] - target)
+        indexer = np.where(distance <= tolerance, indexer, -1)
+        return indexer
+
     def get_loc(self, key, method=None, tolerance=None):
         """Adapted from pandas.tseries.index.DatetimeIndex.get_loc"""
         if isinstance(key, str):
@@ -427,9 +454,11 @@ class CFTimeIndex(pd.Index):
         return CFTimeIndex(other + np.array(self))
 
     def __sub__(self, other):
-        import cftime
-
-        if isinstance(other, (CFTimeIndex, cftime.datetime)):
+        if _contains_datetime_timedeltas(other):
+            return CFTimeIndex(np.array(self) - other)
+        elif isinstance(other, pd.TimedeltaIndex):
+            return CFTimeIndex(np.array(self) - other.to_pytimedelta())
+        elif _contains_cftime_datetimes(np.array(other)):
             try:
                 return pd.TimedeltaIndex(np.array(self) - np.array(other))
             except OverflowError:
@@ -437,14 +466,17 @@ class CFTimeIndex(pd.Index):
                     "The time difference exceeds the range of values "
                     "that can be expressed at the nanosecond resolution."
                 )
-
-        elif isinstance(other, pd.TimedeltaIndex):
-            return CFTimeIndex(np.array(self) - other.to_pytimedelta())
         else:
-            return CFTimeIndex(np.array(self) - other)
+            return NotImplemented
 
     def __rsub__(self, other):
-        return pd.TimedeltaIndex(other - np.array(self))
+        try:
+            return pd.TimedeltaIndex(other - np.array(self))
+        except OverflowError:
+            raise ValueError(
+                "The time difference exceeds the range of values "
+                "that can be expressed at the nanosecond resolution."
+            )
 
     def to_datetimeindex(self, unsafe=False):
         """If possible, convert this index to a pandas.DatetimeIndex.
@@ -633,6 +665,12 @@ def _parse_array_of_cftime_strings(strings, date_type):
     ).reshape(strings.shape)
 
 
+def _contains_datetime_timedeltas(array):
+    """Check if an input array contains datetime.timedelta objects."""
+    array = np.atleast_1d(array)
+    return isinstance(array[0], timedelta)
+
+
 def _cftimeindex_from_i8(values, date_type, name):
     """Construct a CFTimeIndex from an array of integers.
 
diff --git a/xarray/core/dataarray.py b/xarray/core/dataarray.py
index 7a95aedc..6782070d 100644
--- a/xarray/core/dataarray.py
+++ b/xarray/core/dataarray.py
@@ -1475,20 +1475,23 @@ class DataArray(AbstractArray, DataWithCoords):
 
         Examples
         --------
+
         >>> arr = xr.DataArray(data=[0, 1], dims="x",
-                               coords={"x": ["a", "b"], "y": ("x", [0, 1])})
+        ...                    coords={"x": ["a", "b"], "y": ("x", [0, 1])})
         >>> arr
         <xarray.DataArray (x: 2)>
         array([0, 1])
         Coordinates:
           * x        (x) <U1 'a' 'b'
             y        (x) int64 0 1
+
         >>> arr.swap_dims({"x": "y"})
         <xarray.DataArray (y: 2)>
         array([0, 1])
         Coordinates:
             x        (y) <U1 'a' 'b'
           * y        (y) int64 0 1
+
         >>> arr.swap_dims({"x": "z"})
         <xarray.DataArray (z: 2)>
         array([0, 1])
@@ -1718,7 +1721,7 @@ class DataArray(AbstractArray, DataWithCoords):
         Examples
         --------
 
-        >>> arr = DataArray(np.arange(6).reshape(2, 3),
+        >>> arr = xr.DataArray(np.arange(6).reshape(2, 3),
         ...                 coords=[('x', ['a', 'b']), ('y', [0, 1, 2])])
         >>> arr
         <xarray.DataArray (x: 2, y: 3)>
@@ -1768,7 +1771,7 @@ class DataArray(AbstractArray, DataWithCoords):
         Examples
         --------
 
-        >>> arr = DataArray(np.arange(6).reshape(2, 3),
+        >>> arr = xr.DataArray(np.arange(6).reshape(2, 3),
         ...                 coords=[('x', ['a', 'b']), ('y', [0, 1, 2])])
         >>> arr
         <xarray.DataArray (x: 2, y: 3)>
@@ -1817,7 +1820,7 @@ class DataArray(AbstractArray, DataWithCoords):
         Examples
         --------
         >>> import xarray as xr
-        >>> arr = DataArray(np.arange(6).reshape(2, 3),
+        >>> arr = xr.DataArray(np.arange(6).reshape(2, 3),
         ...                 coords=[('x', ['a', 'b']), ('y', [0, 1, 2])])
         >>> data = xr.Dataset({'a': arr, 'b': arr.isel(y=0)})
         >>> data
@@ -2623,7 +2626,7 @@ class DataArray(AbstractArray, DataWithCoords):
         """
         Access plotting functions for DataArray's
 
-        >>> d = DataArray([[1, 2], [3, 4]])
+        >>> d = xr.DataArray([[1, 2], [3, 4]])
 
         For convenience just call this directly
 
@@ -2849,18 +2852,20 @@ class DataArray(AbstractArray, DataWithCoords):
         --------
 
         >>> da_vals = np.arange(6 * 5 * 4).reshape((6, 5, 4))
-        >>> da = DataArray(da_vals, dims=['x', 'y', 'z'])
+        >>> da = xr.DataArray(da_vals, dims=['x', 'y', 'z'])
         >>> dm_vals = np.arange(4)
-        >>> dm = DataArray(dm_vals, dims=['z'])
+        >>> dm = xr.DataArray(dm_vals, dims=['z'])
 
         >>> dm.dims
         ('z')
+
         >>> da.dims
         ('x', 'y', 'z')
 
         >>> dot_result = da.dot(dm)
         >>> dot_result.dims
         ('x', 'y')
+
         """
         if isinstance(other, Dataset):
             raise NotImplementedError(
diff --git a/xarray/core/dataset.py b/xarray/core/dataset.py
index f286236d..a4d20a79 100644
--- a/xarray/core/dataset.py
+++ b/xarray/core/dataset.py
@@ -1011,7 +1011,7 @@ class Dataset(Mapping, ImplementsDatasetReduce, DataWithCoords):
 
         >>> da = xr.DataArray(np.random.randn(2, 3))
         >>> ds = xr.Dataset({'foo': da, 'bar': ('x', [-1, 2])},
-                            coords={'x': ['one', 'two']})
+        ...                 coords={'x': ['one', 'two']})
         >>> ds.copy()
         <xarray.Dataset>
         Dimensions:  (dim_0: 2, dim_1: 3, x: 2)
@@ -1021,6 +1021,7 @@ class Dataset(Mapping, ImplementsDatasetReduce, DataWithCoords):
         Data variables:
             foo      (dim_0, dim_1) float64 -0.8079 0.3897 -1.862 -0.6091 -1.051 -0.3003
             bar      (x) int64 -1 2
+
         >>> ds_0 = ds.copy(deep=False)
         >>> ds_0['foo'][0, 0] = 7
         >>> ds_0
@@ -1032,6 +1033,7 @@ class Dataset(Mapping, ImplementsDatasetReduce, DataWithCoords):
         Data variables:
             foo      (dim_0, dim_1) float64 7.0 0.3897 -1.862 -0.6091 -1.051 -0.3003
             bar      (x) int64 -1 2
+
         >>> ds
         <xarray.Dataset>
         Dimensions:  (dim_0: 2, dim_1: 3, x: 2)
@@ -1055,6 +1057,7 @@ class Dataset(Mapping, ImplementsDatasetReduce, DataWithCoords):
         Data variables:
             foo      (dim_0, dim_1) int64 0 1 2 3 4 5
             bar      (x) <U1 'a' 'b'
+
         >>> ds
         <xarray.Dataset>
         Dimensions:  (dim_0: 2, dim_1: 3, x: 2)
@@ -2883,7 +2886,7 @@ class Dataset(Mapping, ImplementsDatasetReduce, DataWithCoords):
         Examples
         --------
         >>> ds = xr.Dataset(data_vars={"a": ("x", [5, 7]), "b": ("x", [0.1, 2.4])},
-                            coords={"x": ["a", "b"], "y": ("x", [0, 1])})
+        ...                 coords={"x": ["a", "b"], "y": ("x", [0, 1])})
         >>> ds
         <xarray.Dataset>
         Dimensions:  (x: 2)
@@ -2893,6 +2896,7 @@ class Dataset(Mapping, ImplementsDatasetReduce, DataWithCoords):
         Data variables:
             a        (x) int64 5 7
             b        (x) float64 0.1 2.4
+
         >>> ds.swap_dims({"x": "y"})
         <xarray.Dataset>
         Dimensions:  (y: 2)
@@ -2902,6 +2906,7 @@ class Dataset(Mapping, ImplementsDatasetReduce, DataWithCoords):
         Data variables:
             a        (y) int64 5 7
             b        (y) float64 0.1 2.4
+
         >>> ds.swap_dims({"x": "z"})
         <xarray.Dataset>
         Dimensions:  (z: 2)
@@ -3341,7 +3346,7 @@ class Dataset(Mapping, ImplementsDatasetReduce, DataWithCoords):
 
         Examples
         --------
-        >>> data = Dataset(
+        >>> data = xr.Dataset(
         ...     data_vars={'a': (('x', 'y'), [[0, 1, 2], [3, 4, 5]]),
         ...                'b': ('x', [6, 7])},
         ...     coords={'y': ['u', 'v', 'w']}
diff --git a/xarray/core/rolling.py b/xarray/core/rolling.py
index 61178cfb..5f633abb 100644
--- a/xarray/core/rolling.py
+++ b/xarray/core/rolling.py
@@ -231,21 +231,22 @@ class DataArrayRolling(Rolling):
 
         Examples
         --------
-        >>> da = DataArray(np.arange(8).reshape(2, 4), dims=('a', 'b'))
-        >>>
+        >>> da = xr.DataArray(np.arange(8).reshape(2, 4), dims=('a', 'b'))
+
         >>> rolling = da.rolling(b=3)
         >>> rolling.construct('window_dim')
         <xarray.DataArray (a: 2, b: 4, window_dim: 3)>
         array([[[np.nan, np.nan, 0], [np.nan, 0, 1], [0, 1, 2], [1, 2, 3]],
                [[np.nan, np.nan, 4], [np.nan, 4, 5], [4, 5, 6], [5, 6, 7]]])
         Dimensions without coordinates: a, b, window_dim
-        >>>
+
         >>> rolling = da.rolling(b=3, center=True)
         >>> rolling.construct('window_dim')
         <xarray.DataArray (a: 2, b: 4, window_dim: 3)>
         array([[[np.nan, 0, 1], [0, 1, 2], [1, 2, 3], [2, 3, np.nan]],
                [[np.nan, 4, 5], [4, 5, 6], [5, 6, 7], [6, 7, np.nan]]])
         Dimensions without coordinates: a, b, window_dim
+
         """
 
         from .dataarray import DataArray
@@ -278,26 +279,26 @@ class DataArrayRolling(Rolling):
 
         Examples
         --------
-        >>> da = DataArray(np.arange(8).reshape(2, 4), dims=('a', 'b'))
-        >>>
+        >>> da = xr.DataArray(np.arange(8).reshape(2, 4), dims=('a', 'b'))
         >>> rolling = da.rolling(b=3)
         >>> rolling.construct('window_dim')
         <xarray.DataArray (a: 2, b: 4, window_dim: 3)>
         array([[[np.nan, np.nan, 0], [np.nan, 0, 1], [0, 1, 2], [1, 2, 3]],
                [[np.nan, np.nan, 4], [np.nan, 4, 5], [4, 5, 6], [5, 6, 7]]])
         Dimensions without coordinates: a, b, window_dim
-        >>>
+
         >>> rolling.reduce(np.sum)
         <xarray.DataArray (a: 2, b: 4)>
         array([[nan, nan,  3.,  6.],
                [nan, nan, 15., 18.]])
         Dimensions without coordinates: a, b
-        >>>
+
         >>> rolling = da.rolling(b=3, min_periods=1)
         >>> rolling.reduce(np.nansum)
         <xarray.DataArray (a: 2, b: 4)>
         array([[ 0.,  1.,  3.,  6.],
                [ 4.,  9., 15., 18.]])
+
         """
         rolling_dim = utils.get_temp_dimname(self.obj.dims, "_rolling_dim")
         windows = self.construct(rolling_dim)


Test patch:
diff --git a/conftest.py b/conftest.py
index 25dc2849..712af1d3 100644
--- a/conftest.py
+++ b/conftest.py
@@ -21,3 +21,14 @@ def pytest_runtest_setup(item):
         pytest.skip(
             "set --run-network-tests to run test requiring an " "internet connection"
         )
+
+
+@pytest.fixture(autouse=True)
+def add_standard_imports(doctest_namespace):
+    import numpy as np
+    import pandas as pd
+    import xarray as xr
+
+    doctest_namespace["np"] = np
+    doctest_namespace["pd"] = pd
+    doctest_namespace["xr"] = xr
diff --git a/xarray/tests/test_cftimeindex.py b/xarray/tests/test_cftimeindex.py
index 8d83b833..43d6d7b0 100644
--- a/xarray/tests/test_cftimeindex.py
+++ b/xarray/tests/test_cftimeindex.py
@@ -451,10 +451,21 @@ def test_sel_date_scalar(da, date_type, index):
 
 
 @pytest.mark.xfail(reason="https://github.com/pydata/xarray/issues/3751")
+@requires_cftime
+def test_sel_date_distant_date(da, date_type, index):
+    expected = xr.DataArray(4).assign_coords(time=index[3])
+    result = da.sel(time=date_type(2000, 1, 1), method="nearest")
+    assert_identical(result, expected)
+
+
 @requires_cftime
 @pytest.mark.parametrize(
     "sel_kwargs",
-    [{"method": "nearest"}, {"method": "nearest", "tolerance": timedelta(days=70)}],
+    [
+        {"method": "nearest"},
+        {"method": "nearest", "tolerance": timedelta(days=70)},
+        {"method": "nearest", "tolerance": timedelta(days=1800000)},
+    ],
 )
 def test_sel_date_scalar_nearest(da, date_type, index, sel_kwargs):
     expected = xr.DataArray(2).assign_coords(time=index[1])
@@ -738,7 +749,7 @@ def test_timedeltaindex_add_cftimeindex(calendar):
 
 
 @requires_cftime
-def test_cftimeindex_sub(index):
+def test_cftimeindex_sub_timedelta(index):
     date_type = index.date_type
     expected_dates = [
         date_type(1, 1, 2),
@@ -753,6 +764,27 @@ def test_cftimeindex_sub(index):
     assert isinstance(result, CFTimeIndex)
 
 
+@requires_cftime
+@pytest.mark.parametrize(
+    "other",
+    [np.array(4 * [timedelta(days=1)]), np.array(timedelta(days=1))],
+    ids=["1d-array", "scalar-array"],
+)
+def test_cftimeindex_sub_timedelta_array(index, other):
+    date_type = index.date_type
+    expected_dates = [
+        date_type(1, 1, 2),
+        date_type(1, 2, 2),
+        date_type(2, 1, 2),
+        date_type(2, 2, 2),
+    ]
+    expected = CFTimeIndex(expected_dates)
+    result = index + timedelta(days=2)
+    result = result - other
+    assert result.equals(expected)
+    assert isinstance(result, CFTimeIndex)
+
+
 @requires_cftime
 @pytest.mark.parametrize("calendar", _CFTIME_CALENDARS)
 def test_cftimeindex_sub_cftimeindex(calendar):
@@ -784,6 +816,14 @@ def test_cftime_datetime_sub_cftimeindex(calendar):
     assert isinstance(result, pd.TimedeltaIndex)
 
 
+@requires_cftime
+@pytest.mark.parametrize("calendar", _CFTIME_CALENDARS)
+def test_distant_cftime_datetime_sub_cftimeindex(calendar):
+    a = xr.cftime_range("2000", periods=5, calendar=calendar)
+    with pytest.raises(ValueError, match="difference exceeds"):
+        a.date_type(1, 1, 1) - a
+
+
 @requires_cftime
 @pytest.mark.parametrize("calendar", _CFTIME_CALENDARS)
 def test_cftimeindex_sub_timedeltaindex(calendar):
@@ -795,6 +835,25 @@ def test_cftimeindex_sub_timedeltaindex(calendar):
     assert isinstance(result, CFTimeIndex)
 
 
+@requires_cftime
+@pytest.mark.parametrize("calendar", _CFTIME_CALENDARS)
+def test_cftimeindex_sub_index_of_cftime_datetimes(calendar):
+    a = xr.cftime_range("2000", periods=5, calendar=calendar)
+    b = pd.Index(a.values)
+    expected = a - a
+    result = a - b
+    assert result.equals(expected)
+    assert isinstance(result, pd.TimedeltaIndex)
+
+
+@requires_cftime
+@pytest.mark.parametrize("calendar", _CFTIME_CALENDARS)
+def test_cftimeindex_sub_not_implemented(calendar):
+    a = xr.cftime_range("2000", periods=5, calendar=calendar)
+    with pytest.raises(TypeError, match="unsupported operand"):
+        a - 1
+
+
 @requires_cftime
 def test_cftimeindex_rsub(index):
     with pytest.raises(TypeError):
diff --git a/xarray/tests/test_units.py b/xarray/tests/test_units.py
index 9f63ebb1..bef3af62 100644
--- a/xarray/tests/test_units.py
+++ b/xarray/tests/test_units.py
@@ -1,3 +1,4 @@
+import functools
 import operator
 from distutils.version import LooseVersion
 
@@ -8,6 +9,7 @@ import pytest
 import xarray as xr
 from xarray.core import formatting
 from xarray.core.npcompat import IS_NEP18_ACTIVE
+from xarray.testing import assert_allclose, assert_identical
 
 from .test_variable import VariableSubclassobjects
 
@@ -70,53 +72,17 @@ def array_strip_units(array):
         return array
 
 
-def array_attach_units(data, unit, convert_from=None):
-    try:
-        unit, convert_from = unit
-    except TypeError:
-        pass
-
+def array_attach_units(data, unit):
     if isinstance(data, Quantity):
-        if not convert_from:
-            raise ValueError(
-                "cannot attach unit {unit} to quantity ({data.units})".format(
-                    unit=unit, data=data
-                )
-            )
-        elif isinstance(convert_from, unit_registry.Unit):
-            data = data.magnitude
-        elif convert_from is True:  # intentionally accept exactly true
-            if data.check(unit):
-                convert_from = data.units
-                data = data.magnitude
-            else:
-                raise ValueError(
-                    "cannot convert quantity ({data.units}) to {unit}".format(
-                        unit=unit, data=data
-                    )
-                )
-        else:
-            raise ValueError(
-                "cannot convert from invalid unit {convert_from}".format(
-                    convert_from=convert_from
-                )
-            )
+        raise ValueError(f"cannot attach unit {unit} to quantity {data}")
 
-    # to make sure we also encounter the case of "equal if converted"
-    if convert_from is not None:
-        quantity = (data * convert_from).to(
-            unit
-            if isinstance(unit, unit_registry.Unit)
-            else unit_registry.dimensionless
-        )
-    else:
-        try:
-            quantity = data * unit
-        except np.core._exceptions.UFuncTypeError:
-            if unit != 1:
-                raise
+    try:
+        quantity = data * unit
+    except np.core._exceptions.UFuncTypeError:
+        if isinstance(unit, unit_registry.Unit):
+            raise
 
-            quantity = data
+        quantity = data
 
     return quantity
 
@@ -241,6 +207,11 @@ def attach_units(obj, units):
 
 
 def convert_units(obj, to):
+    # preprocess
+    to = {
+        key: None if not isinstance(value, unit_registry.Unit) else value
+        for key, value in to.items()
+    }
     if isinstance(obj, xr.Dataset):
         data_vars = {
             name: convert_units(array.variable, {None: to.get(name)})
@@ -282,6 +253,7 @@ def convert_units(obj, to):
 
 
 def assert_units_equal(a, b):
+    __tracebackhide__ = True
     assert extract_units(a) == extract_units(b)
 
 
@@ -414,9 +386,8 @@ class function:
         return f"function_{self.name}"
 
 
-@pytest.mark.xfail(reason="test bug: apply_ufunc should not be called that way")
 def test_apply_ufunc_dataarray(dtype):
-    func = function(
+    func = functools.partial(
         xr.apply_ufunc, np.mean, input_core_dims=[["x"]], kwargs={"axis": -1}
     )
 
@@ -427,12 +398,12 @@ def test_apply_ufunc_dataarray(dtype):
     expected = attach_units(func(strip_units(data_array)), extract_units(data_array))
     actual = func(data_array)
 
-    assert_equal_with_units(expected, actual)
+    assert_units_equal(expected, actual)
+    assert_identical(expected, actual)
 
 
-@pytest.mark.xfail(reason="test bug: apply_ufunc should not be called that way")
 def test_apply_ufunc_dataset(dtype):
-    func = function(
+    func = functools.partial(
         xr.apply_ufunc, np.mean, input_core_dims=[["x"]], kwargs={"axis": -1}
     )
 
@@ -450,10 +421,10 @@ def test_apply_ufunc_dataset(dtype):
     expected = attach_units(func(strip_units(ds)), extract_units(ds))
     actual = func(ds)
 
-    assert_equal_with_units(expected, actual)
+    assert_units_equal(expected, actual)
+    assert_identical(expected, actual)
 
 
-@pytest.mark.xfail(reason="blocked by `reindex` / `where`")
 @pytest.mark.parametrize(
     "unit,error",
     (
@@ -475,36 +446,40 @@ def test_apply_ufunc_dataset(dtype):
         "coords",
     ),
 )
-@pytest.mark.parametrize("fill_value", (np.float64(10), np.float64(np.nan)))
+@pytest.mark.parametrize("fill_value", (10, np.nan))
 def test_align_dataarray(fill_value, variant, unit, error, dtype):
     original_unit = unit_registry.m
 
     variants = {
-        "data": (unit, 1, 1),
-        "dims": (original_unit, unit, 1),
-        "coords": (original_unit, 1, unit),
+        "data": (unit, original_unit, original_unit),
+        "dims": (original_unit, unit, original_unit),
+        "coords": (original_unit, original_unit, unit),
     }
     data_unit, dim_unit, coord_unit = variants.get(variant)
 
     array1 = np.linspace(0, 10, 2 * 5).reshape(2, 5).astype(dtype) * original_unit
     array2 = np.linspace(0, 8, 2 * 5).reshape(2, 5).astype(dtype) * data_unit
     x = np.arange(2) * original_unit
-    x_a1 = np.array([10, 5]) * original_unit
-    x_a2 = np.array([10, 5]) * coord_unit
 
     y1 = np.arange(5) * original_unit
     y2 = np.arange(2, 7) * dim_unit
+    y_a1 = np.array([3, 5, 7, 8, 9]) * original_unit
+    y_a2 = np.array([7, 8, 9, 11, 13]) * coord_unit
 
-    data_array1 = xr.DataArray(
-        data=array1, coords={"x": x, "x_a": ("x", x_a1), "y": y1}, dims=("x", "y")
-    )
-    data_array2 = xr.DataArray(
-        data=array2, coords={"x": x, "x_a": ("x", x_a2), "y": y2}, dims=("x", "y")
-    )
+    coords1 = {"x": x, "y": y1}
+    coords2 = {"x": x, "y": y2}
+    if variant == "coords":
+        coords1["y_a"] = ("y", y_a1)
+        coords2["y_a"] = ("y", y_a2)
+
+    data_array1 = xr.DataArray(data=array1, coords=coords1, dims=("x", "y"))
+    data_array2 = xr.DataArray(data=array2, coords=coords2, dims=("x", "y"))
 
     fill_value = fill_value * data_unit
     func = function(xr.align, join="outer", fill_value=fill_value)
-    if error is not None:
+    if error is not None and not (
+        np.isnan(fill_value) and not isinstance(fill_value, Quantity)
+    ):
         with pytest.raises(error):
             func(data_array1, data_array2)
 
@@ -524,15 +499,19 @@ def test_align_dataarray(fill_value, variant, unit, error, dtype):
         **stripped_kwargs,
     )
     expected_a = attach_units(expected_a, units_a)
-    expected_b = convert_units(attach_units(expected_b, units_a), units_b)
+    if isinstance(array2, Quantity):
+        expected_b = convert_units(attach_units(expected_b, units_a), units_b)
+    else:
+        expected_b = attach_units(expected_b, units_b)
 
     actual_a, actual_b = func(data_array1, data_array2)
 
-    assert_equal_with_units(expected_a, actual_a)
-    assert_equal_with_units(expected_b, actual_b)
+    assert_units_equal(expected_a, actual_a)
+    assert_allclose(expected_a, actual_a)
+    assert_units_equal(expected_b, actual_b)
+    assert_allclose(expected_b, actual_b)
 
 
-@pytest.mark.xfail(reason="blocked by `reindex` / `where`")
 @pytest.mark.parametrize(
     "unit,error",
     (
@@ -558,31 +537,37 @@ def test_align_dataarray(fill_value, variant, unit, error, dtype):
 def test_align_dataset(fill_value, unit, variant, error, dtype):
     original_unit = unit_registry.m
 
-    variants = {"data": (unit, 1, 1), "dims": (1, unit, 1), "coords": (1, 1, unit)}
+    variants = {
+        "data": (unit, original_unit, original_unit),
+        "dims": (original_unit, unit, original_unit),
+        "coords": (original_unit, original_unit, unit),
+    }
     data_unit, dim_unit, coord_unit = variants.get(variant)
 
     array1 = np.linspace(0, 10, 2 * 5).reshape(2, 5).astype(dtype) * original_unit
     array2 = np.linspace(0, 10, 2 * 5).reshape(2, 5).astype(dtype) * data_unit
 
     x = np.arange(2) * original_unit
-    x_a1 = np.array([10, 5]) * original_unit
-    x_a2 = np.array([10, 5]) * coord_unit
 
     y1 = np.arange(5) * original_unit
     y2 = np.arange(2, 7) * dim_unit
+    y_a1 = np.array([3, 5, 7, 8, 9]) * original_unit
+    y_a2 = np.array([7, 8, 9, 11, 13]) * coord_unit
 
-    ds1 = xr.Dataset(
-        data_vars={"a": (("x", "y"), array1)},
-        coords={"x": x, "x_a": ("x", x_a1), "y": y1},
-    )
-    ds2 = xr.Dataset(
-        data_vars={"a": (("x", "y"), array2)},
-        coords={"x": x, "x_a": ("x", x_a2), "y": y2},
-    )
+    coords1 = {"x": x, "y": y1}
+    coords2 = {"x": x, "y": y2}
+    if variant == "coords":
+        coords1["y_a"] = ("y", y_a1)
+        coords2["y_a"] = ("y", y_a2)
+
+    ds1 = xr.Dataset(data_vars={"a": (("x", "y"), array1)}, coords=coords1)
+    ds2 = xr.Dataset(data_vars={"a": (("x", "y"), array2)}, coords=coords2)
 
     fill_value = fill_value * data_unit
     func = function(xr.align, join="outer", fill_value=fill_value)
-    if error is not None:
+    if error is not None and not (
+        np.isnan(fill_value) and not isinstance(fill_value, Quantity)
+    ):
         with pytest.raises(error):
             func(ds1, ds2)
 
@@ -600,12 +585,17 @@ def test_align_dataset(fill_value, unit, variant, error, dtype):
         strip_units(ds1), strip_units(convert_units(ds2, units_a)), **stripped_kwargs
     )
     expected_a = attach_units(expected_a, units_a)
-    expected_b = convert_units(attach_units(expected_b, units_a), units_b)
+    if isinstance(array2, Quantity):
+        expected_b = convert_units(attach_units(expected_b, units_a), units_b)
+    else:
+        expected_b = attach_units(expected_b, units_b)
 
     actual_a, actual_b = func(ds1, ds2)
 
-    assert_equal_with_units(expected_a, actual_a)
-    assert_equal_with_units(expected_b, actual_b)
+    assert_units_equal(expected_a, actual_a)
+    assert_allclose(expected_a, actual_a)
+    assert_units_equal(expected_b, actual_b)
+    assert_allclose(expected_b, actual_b)
 
 
 def test_broadcast_dataarray(dtype):
@@ -615,28 +605,53 @@ def test_broadcast_dataarray(dtype):
     a = xr.DataArray(data=array1, dims="x")
     b = xr.DataArray(data=array2, dims="y")
 
-    expected_a, expected_b = tuple(
-        attach_units(elem, extract_units(a))
-        for elem in xr.broadcast(strip_units(a), strip_units(b))
-    )
+    units_a = extract_units(a)
+    units_b = extract_units(b)
+    expected_a, expected_b = xr.broadcast(strip_units(a), strip_units(b))
+    expected_a = attach_units(expected_a, units_a)
+    expected_b = convert_units(attach_units(expected_b, units_a), units_b)
+
     actual_a, actual_b = xr.broadcast(a, b)
 
-    assert_equal_with_units(expected_a, actual_a)
-    assert_equal_with_units(expected_b, actual_b)
+    assert_units_equal(expected_a, actual_a)
+    assert_identical(expected_a, actual_a)
+    assert_units_equal(expected_b, actual_b)
+    assert_identical(expected_b, actual_b)
 
 
 def test_broadcast_dataset(dtype):
     array1 = np.linspace(0, 10, 2) * unit_registry.Pa
     array2 = np.linspace(0, 10, 3) * unit_registry.Pa
 
-    ds = xr.Dataset(data_vars={"a": ("x", array1), "b": ("y", array2)})
+    x1 = np.arange(2)
+    y1 = np.arange(3)
+
+    x2 = np.arange(2, 4)
+    y2 = np.arange(3, 6)
 
-    (expected,) = tuple(
-        attach_units(elem, extract_units(ds)) for elem in xr.broadcast(strip_units(ds))
+    ds = xr.Dataset(
+        data_vars={"a": ("x", array1), "b": ("y", array2)}, coords={"x": x1, "y": y1}
+    )
+    other = xr.Dataset(
+        data_vars={
+            "a": ("x", array1.to(unit_registry.hPa)),
+            "b": ("y", array2.to(unit_registry.hPa)),
+        },
+        coords={"x": x2, "y": y2},
     )
-    (actual,) = xr.broadcast(ds)
 
-    assert_equal_with_units(expected, actual)
+    units_a = extract_units(ds)
+    units_b = extract_units(other)
+    expected_a, expected_b = xr.broadcast(strip_units(ds), strip_units(other))
+    expected_a = attach_units(expected_a, units_a)
+    expected_b = attach_units(expected_b, units_b)
+
+    actual_a, actual_b = xr.broadcast(ds, other)
+
+    assert_units_equal(expected_a, actual_a)
+    assert_identical(expected_a, actual_a)
+    assert_units_equal(expected_b, actual_b)
+    assert_identical(expected_b, actual_b)
 
 
 @pytest.mark.parametrize(
@@ -706,7 +721,8 @@ def test_combine_by_coords(variant, unit, error, dtype):
     )
     actual = xr.combine_by_coords([ds, other])
 
-    assert_equal_with_units(expected, actual)
+    assert_units_equal(expected, actual)
+    assert_identical(expected, actual)
 
 
 @pytest.mark.parametrize(
@@ -717,12 +733,7 @@ def test_combine_by_coords(variant, unit, error, dtype):
             unit_registry.dimensionless, DimensionalityError, id="dimensionless"
         ),
         pytest.param(unit_registry.s, DimensionalityError, id="incompatible_unit"),
-        pytest.param(
-            unit_registry.mm,
-            None,
-            id="compatible_unit",
-            marks=pytest.mark.xfail(reason="wrong order of arguments to `where`"),
-        ),
+        pytest.param(unit_registry.mm, None, id="compatible_unit"),
         pytest.param(unit_registry.m, None, id="identical_unit"),
     ),
     ids=repr,
@@ -810,7 +821,8 @@ def test_combine_nested(variant, unit, error, dtype):
     )
     actual = func([[ds1, ds2], [ds3, ds4]])
 
-    assert_equal_with_units(expected, actual)
+    assert_units_equal(expected, actual)
+    assert_identical(expected, actual)
 
 
 @pytest.mark.parametrize(
@@ -862,7 +874,8 @@ def test_concat_dataarray(variant, unit, error, dtype):
     )
     actual = xr.concat([arr1, arr2], dim="x")
 
-    assert_equal_with_units(expected, actual)
+    assert_units_equal(expected, actual)
+    assert_identical(expected, actual)
 
 
 @pytest.mark.parametrize(
@@ -912,10 +925,10 @@ def test_concat_dataset(variant, unit, error, dtype):
     )
     actual = xr.concat([ds1, ds2], dim="x")
 
-    assert_equal_with_units(expected, actual)
+    assert_units_equal(expected, actual)
+    assert_identical(expected, actual)
 
 
-@pytest.mark.xfail(reason="blocked by `reindex` / `where`")
 @pytest.mark.parametrize(
     "unit,error",
     (
@@ -948,64 +961,81 @@ def test_merge_dataarray(variant, unit, error, dtype):
     data_unit, dim_unit, coord_unit = variants.get(variant)
 
     array1 = np.linspace(0, 1, 2 * 3).reshape(2, 3).astype(dtype) * original_unit
+    x1 = np.arange(2) * original_unit
+    y1 = np.arange(3) * original_unit
+    u1 = np.linspace(10, 20, 2) * original_unit
+    v1 = np.linspace(10, 20, 3) * original_unit
+
     array2 = np.linspace(1, 2, 2 * 4).reshape(2, 4).astype(dtype) * data_unit
-    array3 = np.linspace(0, 2, 3 * 4).reshape(3, 4).astype(dtype) * data_unit
+    x2 = np.arange(2, 4) * dim_unit
+    z2 = np.arange(4) * original_unit
+    u2 = np.linspace(20, 30, 2) * coord_unit
+    w2 = np.linspace(10, 20, 4) * original_unit
 
-    x = np.arange(2) * original_unit
-    y = np.arange(3) * original_unit
-    z = np.arange(4) * original_unit
-    u = np.linspace(10, 20, 2) * original_unit
-    v = np.linspace(10, 20, 3) * original_unit
-    w = np.linspace(10, 20, 4) * original_unit
+    array3 = np.linspace(0, 2, 3 * 4).reshape(3, 4).astype(dtype) * data_unit
+    y3 = np.arange(3, 6) * dim_unit
+    z3 = np.arange(4, 8) * dim_unit
+    v3 = np.linspace(10, 20, 3) * coord_unit
+    w3 = np.linspace(10, 20, 4) * coord_unit
 
     arr1 = xr.DataArray(
         name="a",
         data=array1,
-        coords={"x": x, "y": y, "u": ("x", u), "v": ("y", v)},
+        coords={"x": x1, "y": y1, "u": ("x", u1), "v": ("y", v1)},
         dims=("x", "y"),
     )
     arr2 = xr.DataArray(
-        name="b",
+        name="a",
         data=array2,
-        coords={
-            "x": np.arange(2, 4) * dim_unit,
-            "z": z,
-            "u": ("x", np.linspace(20, 30, 2) * coord_unit),
-            "w": ("z", w),
-        },
+        coords={"x": x2, "z": z2, "u": ("x", u2), "w": ("z", w2)},
         dims=("x", "z"),
     )
     arr3 = xr.DataArray(
-        name="c",
+        name="a",
         data=array3,
-        coords={
-            "y": np.arange(3, 6) * dim_unit,
-            "z": np.arange(4, 8) * dim_unit,
-            "v": ("y", np.linspace(10, 20, 3) * coord_unit),
-            "w": ("z", np.linspace(10, 20, 4) * coord_unit),
-        },
+        coords={"y": y3, "z": z3, "v": ("y", v3), "w": ("z", w3)},
         dims=("y", "z"),
     )
 
-    func = function(xr.merge)
     if error is not None:
         with pytest.raises(error):
-            func([arr1, arr2, arr3])
+            xr.merge([arr1, arr2, arr3])
 
         return
 
-    units = {name: original_unit for name in list("abcuvwxyz")}
+    units = {name: original_unit for name in list("axyzuvw")}
+
     convert_and_strip = lambda arr: strip_units(convert_units(arr, units))
-    expected = attach_units(
-        func([strip_units(arr1), convert_and_strip(arr2), convert_and_strip(arr3)]),
-        units,
+    expected_units = {
+        "a": original_unit,
+        "u": original_unit,
+        "v": original_unit,
+        "w": original_unit,
+        "x": original_unit,
+        "y": original_unit,
+        "z": original_unit,
+    }
+
+    expected = convert_units(
+        attach_units(
+            xr.merge(
+                [
+                    convert_and_strip(arr1),
+                    convert_and_strip(arr2),
+                    convert_and_strip(arr3),
+                ]
+            ),
+            units,
+        ),
+        expected_units,
     )
-    actual = func([arr1, arr2, arr3])
 
-    assert_equal_with_units(expected, actual)
+    actual = xr.merge([arr1, arr2, arr3])
+
+    assert_units_equal(expected, actual)
+    assert_allclose(expected, actual)
 
 
-@pytest.mark.xfail(reason="blocked by `reindex` / `where`")
 @pytest.mark.parametrize(
     "unit,error",
     (
@@ -1046,7 +1076,7 @@ def test_merge_dataset(variant, unit, error, dtype):
 
     ds1 = xr.Dataset(
         data_vars={"a": (("y", "x"), array1), "b": (("y", "x"), array2)},
-        coords={"x": x, "y": y, "z": ("x", z)},
+        coords={"x": x, "y": y, "u": ("x", z)},
     )
     ds2 = xr.Dataset(
         data_vars={
@@ -1056,18 +1086,18 @@ def test_merge_dataset(variant, unit, error, dtype):
         coords={
             "x": np.arange(3) * dim_unit,
             "y": np.arange(2, 4) * dim_unit,
-            "z": ("x", np.arange(-3, 0) * coord_unit),
+            "u": ("x", np.arange(-3, 0) * coord_unit),
         },
     )
     ds3 = xr.Dataset(
         data_vars={
-            "a": (("y", "x"), np.zeros_like(array1) * np.nan * data_unit),
-            "b": (("y", "x"), np.zeros_like(array2) * np.nan * data_unit),
+            "a": (("y", "x"), np.full_like(array1, np.nan) * data_unit),
+            "b": (("y", "x"), np.full_like(array2, np.nan) * data_unit),
         },
         coords={
             "x": np.arange(3, 6) * dim_unit,
             "y": np.arange(4, 6) * dim_unit,
-            "z": ("x", np.arange(3, 6) * coord_unit),
+            "u": ("x", np.arange(3, 6) * coord_unit),
         },
     )
 
@@ -1080,12 +1110,20 @@ def test_merge_dataset(variant, unit, error, dtype):
 
     units = extract_units(ds1)
     convert_and_strip = lambda ds: strip_units(convert_units(ds, units))
-    expected = attach_units(
-        func([strip_units(ds1), convert_and_strip(ds2), convert_and_strip(ds3)]), units
+    expected_units = {name: original_unit for name in list("abxyzu")}
+    expected = convert_units(
+        attach_units(
+            func(
+                [convert_and_strip(ds1), convert_and_strip(ds2), convert_and_strip(ds3)]
+            ),
+            units,
+        ),
+        expected_units,
     )
     actual = func([ds1, ds2, ds3])
 
-    assert_equal_with_units(expected, actual)
+    assert_units_equal(expected, actual)
+    assert_allclose(expected, actual)
 
 
 @pytest.mark.parametrize("func", (xr.zeros_like, xr.ones_like))
@@ -1094,10 +1132,12 @@ def test_replication_dataarray(func, dtype):
     data_array = xr.DataArray(data=array, dims="x")
 
     numpy_func = getattr(np, func.__name__)
-    expected = xr.DataArray(data=numpy_func(array), dims="x")
+    units = extract_units(numpy_func(data_array))
+    expected = attach_units(func(data_array), units)
     actual = func(data_array)
 
-    assert_equal_with_units(expected, actual)
+    assert_units_equal(expected, actual)
+    assert_identical(expected, actual)
 
 
 @pytest.mark.parametrize("func", (xr.zeros_like, xr.ones_like))
@@ -1114,12 +1154,13 @@ def test_replication_dataset(func, dtype):
     )
 
     numpy_func = getattr(np, func.__name__)
-    expected = ds.copy(
-        data={name: numpy_func(array.data) for name, array in ds.data_vars.items()}
-    )
+    units = extract_units(ds.map(numpy_func))
+    expected = attach_units(func(strip_units(ds)), units)
+
     actual = func(ds)
 
-    assert_equal_with_units(expected, actual)
+    assert_units_equal(expected, actual)
+    assert_identical(expected, actual)
 
 
 @pytest.mark.xfail(
@@ -1158,7 +1199,8 @@ def test_replication_full_like_dataarray(unit, error, dtype):
     )
     actual = xr.full_like(data_array, fill_value=fill_value)
 
-    assert_equal_with_units(expected, actual)
+    assert_units_equal(expected, actual)
+    assert_identical(expected, actual)
 
 
 @pytest.mark.xfail(
@@ -1208,7 +1250,8 @@ def test_replication_full_like_dataset(unit, error, dtype):
     )
     actual = xr.full_like(ds, fill_value=fill_value)
 
-    assert_equal_with_units(expected, actual)
+    assert_units_equal(expected, actual)
+    assert_identical(expected, actual)
 
 
 @pytest.mark.parametrize(
@@ -1250,7 +1293,8 @@ def test_where_dataarray(fill_value, unit, error, dtype):
     )
     actual = xr.where(cond, x, fill_value)
 
-    assert_equal_with_units(expected, actual)
+    assert_units_equal(expected, actual)
+    assert_identical(expected, actual)
 
 
 @pytest.mark.parametrize(
@@ -1294,7 +1338,8 @@ def test_where_dataset(fill_value, unit, error, dtype):
     )
     actual = xr.where(cond, ds, fill_value)
 
-    assert_equal_with_units(expected, actual)
+    assert_units_equal(expected, actual)
+    assert_identical(expected, actual)
 
 
 def test_dot_dataarray(dtype):
@@ -1315,7 +1360,8 @@ def test_dot_dataarray(dtype):
     )
     actual = xr.dot(data_array, other)
 
-    assert_equal_with_units(expected, actual)
+    assert_units_equal(expected, actual)
+    assert_identical(expected, actual)
 
 
 def delete_attrs(*to_delete):