Fix interp bug when indexer shares coordinates with array (#3758)

* added test demonstrating interp bug for nd indexes sharing coordinate with array

* fix test so it works with sel

* support shared dimensions in interp

* isort fixes

* update whats new

* Revert "isort fixes"

This reverts commit 5df6c9c.

* test requires scipy

Author: huard

doc/whats-new.rst

@@ -54,6 +54,9 @@ New Features
 
 Bug fixes
 ~~~~~~~~~
+- Fix :py:meth:`Dataset.interp` when indexing array shares coordinates with the
+  indexed variable (:issue:`3252`).
+  By `David Huard <https://github.com/huard>`_.
 
 - Fix alignment with ``join="override"`` when some dimensions are unindexed. (:issue:`3681`).
   By `Deepak Cherian <https://github.com/dcherian>`_.

xarray/core/dataset.py

@@ -2574,6 +2574,17 @@ def interp(
         coords = either_dict_or_kwargs(coords, coords_kwargs, "interp")
         indexers = dict(self._validate_interp_indexers(coords))
 
+        if coords:
+            # This avoids broadcasting over coordinates that are both in
+            # the original array AND in the indexing array. It essentially
+            # forces interpolation along the shared coordinates.
+            sdims = (
+                set(self.dims)
+                .intersection(*[set(nx.dims) for nx in indexers.values()])
+                .difference(coords.keys())
+            )
+            indexers.update({d: self.variables[d] for d in sdims})
+
         obj = self if assume_sorted else self.sortby([k for k in coords])
 
         def maybe_variable(obj, k):

xarray/tests/test_interp.py

@@ -244,6 +244,36 @@ def test_interpolate_nd(case):
     assert_allclose(actual.transpose("y", "z"), expected)
 
 
+@requires_scipy
+def test_interpolate_nd_nd():
+    """Interpolate nd array with an nd indexer sharing coordinates."""
+    # Create original array
+    a = [0, 2]
+    x = [0, 1, 2]
+    da = xr.DataArray(
+        np.arange(6).reshape(2, 3), dims=("a", "x"), coords={"a": a, "x": x}
+    )
+
+    # Create indexer into `a` with dimensions (y, x)
+    y = [10]
+    c = {"x": x, "y": y}
+    ia = xr.DataArray([[1, 2, 2]], dims=("y", "x"), coords=c)
+    out = da.interp(a=ia)
+    expected = xr.DataArray([[1.5, 4, 5]], dims=("y", "x"), coords=c)
+    xr.testing.assert_allclose(out.drop_vars("a"), expected)
+
+    # If the *shared* indexing coordinates do not match, interp should fail.
+    with pytest.raises(ValueError):
+        c = {"x": [1], "y": y}
+        ia = xr.DataArray([[1]], dims=("y", "x"), coords=c)
+        da.interp(a=ia)
+
+    with pytest.raises(ValueError):
+        c = {"x": [5, 6, 7], "y": y}
+        ia = xr.DataArray([[1]], dims=("y", "x"), coords=c)
+        da.interp(a=ia)
+
+
 @pytest.mark.parametrize("method", ["linear"])
 @pytest.mark.parametrize("case", [0, 1])
 def test_interpolate_scalar(method, case):