Repository: pydata/xarray

Issue:
Add proper tests with fill_value for test_interpolate_pd_compat-tests
### What is your issue?

In #8125 breaking tests in Upstream-CI was reported. Root cause is upstream #https://github.com/pandas-dev/pandas/issues/54920 which is tackled in https://github.com/pandas-dev/pandas/pull/54927.

To get the tests working, the `fill_value`-kwarg has been removed in #8139. The tests didn't cover the use of `fill_value` for scipy based interpolators (the used np.nan is the default-value) completely. 

This issue is to keep track of adding proper tests for `fill_value`.

Discussion:
kmuehlbauer: @max-sixty I've adapted the test to take fill_value into account and left some additional information on pandas' peculiarities. The upstream errors are unrelated. Would be great if you can review.

Patch:
diff --git a/doc/whats-new.rst b/doc/whats-new.rst
index 67429ed7..9a21bcb7 100644
--- a/doc/whats-new.rst
+++ b/doc/whats-new.rst
@@ -90,6 +90,8 @@ Bug fixes
 - ``.rolling_exp`` functions no longer mistakenly lose non-dimensioned coords
   (:issue:`6528`, :pull:`8114`)
   By `Maximilian Roos <https://github.com/max-sixty>`_.
+- In the event that user-provided datetime64/timedelta64 units and integer dtype encoding parameters conflict with each other, override the units to preserve an integer dtype for most faithful serialization to disk (:issue:`1064`, :pull:`8201`).
+  By `Kai Mühlbauer <https://github.com/kmuehlbauer>`_.
 
 Documentation
 ~~~~~~~~~~~~~
@@ -106,6 +108,9 @@ Internal Changes
   than `.reduce`, as the start of a broader effort to move non-reducing
   functions away from ```.reduce``, (:pull:`8114`).
   By `Maximilian Roos <https://github.com/max-sixty>`_.
+- Test range of fill_value's in test_interpolate_pd_compat (:issue:`8146`, :pull:`8189`).
+  By `Kai Mühlbauer <https://github.com/kmuehlbauer>`_.
+
 
 .. _whats-new.2023.08.0:
 
diff --git a/pyproject.toml b/pyproject.toml
index cb51c6ea..25263928 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -96,7 +96,6 @@ module = [
   "fsspec.*",
   "h5netcdf.*",
   "h5py.*",
-  "importlib_metadata.*",
   "iris.*",
   "matplotlib.*",
   "mpl_toolkits.*",
diff --git a/xarray/__init__.py b/xarray/__init__.py
index b63b0d81..1fd3b0c4 100644
--- a/xarray/__init__.py
+++ b/xarray/__init__.py
@@ -1,3 +1,5 @@
+from importlib.metadata import version as _version
+
 from xarray import testing, tutorial
 from xarray.backends.api import (
     load_dataarray,
@@ -41,12 +43,6 @@ from xarray.core.parallel import map_blocks
 from xarray.core.variable import IndexVariable, Variable, as_variable
 from xarray.util.print_versions import show_versions
 
-try:
-    from importlib.metadata import version as _version
-except ImportError:
-    # if the fallback library is missing, we are doomed.
-    from importlib_metadata import version as _version
-
 try:
     __version__ = _version("xarray")
 except Exception:
diff --git a/xarray/coding/times.py b/xarray/coding/times.py
index 79efbecf..2822f02d 100644
--- a/xarray/coding/times.py
+++ b/xarray/coding/times.py
@@ -656,8 +656,22 @@ def cast_to_int_if_safe(num) -> np.ndarray:
     return num
 
 
+def _division(deltas, delta, floor):
+    if floor:
+        # calculate int64 floor division
+        # to preserve integer dtype if possible (GH 4045, GH7817).
+        num = deltas // delta.astype(np.int64)
+        num = num.astype(np.int64, copy=False)
+    else:
+        num = deltas / delta
+    return num
+
+
 def encode_cf_datetime(
-    dates, units: str | None = None, calendar: str | None = None
+    dates,
+    units: str | None = None,
+    calendar: str | None = None,
+    dtype: np.dtype | None = None,
 ) -> tuple[np.ndarray, str, str]:
     """Given an array of datetime objects, returns the tuple `(num, units,
     calendar)` suitable for a CF compliant time variable.
@@ -689,6 +703,12 @@ def encode_cf_datetime(
         time_units, ref_date = _unpack_time_units_and_ref_date(units)
         time_delta = _time_units_to_timedelta64(time_units)
 
+        # Wrap the dates in a DatetimeIndex to do the subtraction to ensure
+        # an OverflowError is raised if the ref_date is too far away from
+        # dates to be encoded (GH 2272).
+        dates_as_index = pd.DatetimeIndex(dates.ravel())
+        time_deltas = dates_as_index - ref_date
+
         # retrieve needed units to faithfully encode to int64
         needed_units, data_ref_date = _unpack_time_units_and_ref_date(data_units)
         if data_units != units:
@@ -697,26 +717,32 @@ def encode_cf_datetime(
             if ref_delta > np.timedelta64(0, "ns"):
                 needed_units = _infer_time_units_from_diff(ref_delta)
 
-        # Wrap the dates in a DatetimeIndex to do the subtraction to ensure
-        # an OverflowError is raised if the ref_date is too far away from
-        # dates to be encoded (GH 2272).
-        dates_as_index = pd.DatetimeIndex(dates.ravel())
-        time_deltas = dates_as_index - ref_date
-
         # needed time delta to encode faithfully to int64
         needed_time_delta = _time_units_to_timedelta64(needed_units)
-        if time_delta <= needed_time_delta:
-            # calculate int64 floor division
-            # to preserve integer dtype if possible (GH 4045, GH7817).
-            num = time_deltas // time_delta.astype(np.int64)
-            num = num.astype(np.int64, copy=False)
-        else:
-            emit_user_level_warning(
-                f"Times can't be serialized faithfully with requested units {units!r}. "
-                f"Resolution of {needed_units!r} needed. "
-                f"Serializing timeseries to floating point."
-            )
-            num = time_deltas / time_delta
+
+        floor_division = True
+        if time_delta > needed_time_delta:
+            floor_division = False
+            if dtype is None:
+                emit_user_level_warning(
+                    f"Times can't be serialized faithfully to int64 with requested units {units!r}. "
+                    f"Resolution of {needed_units!r} needed. Serializing times to floating point instead. "
+                    f"Set encoding['dtype'] to integer dtype to serialize to int64. "
+                    f"Set encoding['dtype'] to floating point dtype to silence this warning."
+                )
+            elif np.issubdtype(dtype, np.integer):
+                new_units = f"{needed_units} since {format_timestamp(ref_date)}"
+                emit_user_level_warning(
+                    f"Times can't be serialized faithfully to int64 with requested units {units!r}. "
+                    f"Serializing with units {new_units!r} instead. "
+                    f"Set encoding['dtype'] to floating point dtype to serialize with units {units!r}. "
+                    f"Set encoding['units'] to {new_units!r} to silence this warning ."
+                )
+                units = new_units
+                time_delta = needed_time_delta
+                floor_division = True
+
+        num = _division(time_deltas, time_delta, floor_division)
         num = num.values.reshape(dates.shape)
 
     except (OutOfBoundsDatetime, OverflowError, ValueError):
@@ -728,7 +754,9 @@ def encode_cf_datetime(
     return (num, units, calendar)
 
 
-def encode_cf_timedelta(timedeltas, units: str | None = None) -> tuple[np.ndarray, str]:
+def encode_cf_timedelta(
+    timedeltas, units: str | None = None, dtype: np.dtype | None = None
+) -> tuple[np.ndarray, str]:
     data_units = infer_timedelta_units(timedeltas)
 
     if units is None:
@@ -744,18 +772,29 @@ def encode_cf_timedelta(timedeltas, units: str | None = None) -> tuple[np.ndarra
 
     # needed time delta to encode faithfully to int64
     needed_time_delta = _time_units_to_timedelta64(needed_units)
-    if time_delta <= needed_time_delta:
-        # calculate int64 floor division
-        # to preserve integer dtype if possible
-        num = time_deltas // time_delta.astype(np.int64)
-        num = num.astype(np.int64, copy=False)
-    else:
-        emit_user_level_warning(
-            f"Timedeltas can't be serialized faithfully with requested units {units!r}. "
-            f"Resolution of {needed_units!r} needed. "
-            f"Serializing timedeltas to floating point."
-        )
-        num = time_deltas / time_delta
+
+    floor_division = True
+    if time_delta > needed_time_delta:
+        floor_division = False
+        if dtype is None:
+            emit_user_level_warning(
+                f"Timedeltas can't be serialized faithfully to int64 with requested units {units!r}. "
+                f"Resolution of {needed_units!r} needed. Serializing timeseries to floating point instead. "
+                f"Set encoding['dtype'] to integer dtype to serialize to int64. "
+                f"Set encoding['dtype'] to floating point dtype to silence this warning."
+            )
+        elif np.issubdtype(dtype, np.integer):
+            emit_user_level_warning(
+                f"Timedeltas can't be serialized faithfully with requested units {units!r}. "
+                f"Serializing with units {needed_units!r} instead. "
+                f"Set encoding['dtype'] to floating point dtype to serialize with units {units!r}. "
+                f"Set encoding['units'] to {needed_units!r} to silence this warning ."
+            )
+            units = needed_units
+            time_delta = needed_time_delta
+            floor_division = True
+
+    num = _division(time_deltas, time_delta, floor_division)
     num = num.values.reshape(timedeltas.shape)
     return (num, units)
 
@@ -772,7 +811,8 @@ class CFDatetimeCoder(VariableCoder):
 
             units = encoding.pop("units", None)
             calendar = encoding.pop("calendar", None)
-            (data, units, calendar) = encode_cf_datetime(data, units, calendar)
+            dtype = encoding.get("dtype", None)
+            (data, units, calendar) = encode_cf_datetime(data, units, calendar, dtype)
 
             safe_setitem(attrs, "units", units, name=name)
             safe_setitem(attrs, "calendar", calendar, name=name)
@@ -807,7 +847,9 @@ class CFTimedeltaCoder(VariableCoder):
         if np.issubdtype(variable.data.dtype, np.timedelta64):
             dims, data, attrs, encoding = unpack_for_encoding(variable)
 
-            data, units = encode_cf_timedelta(data, encoding.pop("units", None))
+            data, units = encode_cf_timedelta(
+                data, encoding.pop("units", None), encoding.get("dtype", None)
+            )
             safe_setitem(attrs, "units", units, name=name)
 
             return Variable(dims, data, attrs, encoding, fastpath=True)
diff --git a/xarray/core/computation.py b/xarray/core/computation.py
index 971f036b..bae779af 100644
--- a/xarray/core/computation.py
+++ b/xarray/core/computation.py
@@ -893,7 +893,7 @@ def apply_ufunc(
     dataset_fill_value: object = _NO_FILL_VALUE,
     keep_attrs: bool | str | None = None,
     kwargs: Mapping | None = None,
-    dask: str = "forbidden",
+    dask: Literal["forbidden", "allowed", "parallelized"] = "forbidden",
     output_dtypes: Sequence | None = None,
     output_sizes: Mapping[Any, int] | None = None,
     meta: Any = None,
diff --git a/xarray/core/coordinates.py b/xarray/core/coordinates.py
index 97ba383e..0c85b2a2 100644
--- a/xarray/core/coordinates.py
+++ b/xarray/core/coordinates.py
@@ -23,7 +23,7 @@ from xarray.core.indexes import (
     create_default_index_implicit,
 )
 from xarray.core.merge import merge_coordinates_without_align, merge_coords
-from xarray.core.types import Self, T_DataArray, T_Xarray
+from xarray.core.types import DataVars, Self, T_DataArray, T_Xarray
 from xarray.core.utils import (
     Frozen,
     ReprObject,
@@ -937,7 +937,7 @@ def assert_coordinate_consistent(obj: T_Xarray, coords: Mapping[Any, Variable])
 
 
 def create_coords_with_default_indexes(
-    coords: Mapping[Any, Any], data_vars: Mapping[Any, Any] | None = None
+    coords: Mapping[Any, Any], data_vars: DataVars | None = None
 ) -> Coordinates:
     """Returns a Coordinates object from a mapping of coordinates (arbitrary objects).
 
diff --git a/xarray/core/dataarray.py b/xarray/core/dataarray.py
index 73464c07..724a5fc2 100644
--- a/xarray/core/dataarray.py
+++ b/xarray/core/dataarray.py
@@ -128,7 +128,7 @@ def _check_coords_dims(shape, coords, dims):
                 f"dimensions {dims}"
             )
 
-        for d, s in zip(v.dims, v.shape):
+        for d, s in v.sizes.items():
             if s != sizes[d]:
                 raise ValueError(
                     f"conflicting sizes for dimension {d!r}: "
@@ -136,13 +136,6 @@ def _check_coords_dims(shape, coords, dims):
                     f"coordinate {k!r}"
                 )
 
-        if k in sizes and v.shape != (sizes[k],):
-            raise ValueError(
-                f"coordinate {k!r} is a DataArray dimension, but "
-                f"it has shape {v.shape!r} rather than expected shape {sizes[k]!r} "
-                "matching the dimension size"
-            )
-
 
 def _infer_coords_and_dims(
     shape, coords, dims
diff --git a/xarray/core/dataset.py b/xarray/core/dataset.py
index 9d771f03..44016e87 100644
--- a/xarray/core/dataset.py
+++ b/xarray/core/dataset.py
@@ -124,7 +124,7 @@ if TYPE_CHECKING:
     from xarray.backends.api import T_NetcdfEngine, T_NetcdfTypes
     from xarray.core.dataarray import DataArray
     from xarray.core.groupby import DatasetGroupBy
-    from xarray.core.merge import CoercibleMapping, CoercibleValue
+    from xarray.core.merge import CoercibleMapping, CoercibleValue, _MergeResult
     from xarray.core.parallelcompat import ChunkManagerEntrypoint
     from xarray.core.resample import DatasetResample
     from xarray.core.rolling import DatasetCoarsen, DatasetRolling
@@ -133,6 +133,7 @@ if TYPE_CHECKING:
         CoarsenBoundaryOptions,
         CombineAttrsOptions,
         CompatOptions,
+        DataVars,
         DatetimeLike,
         DatetimeUnitOptions,
         Dims,
@@ -404,7 +405,7 @@ def _initialize_curvefit_params(params, p0, bounds, func_args):
     return param_defaults, bounds_defaults
 
 
-def merge_data_and_coords(data_vars, coords):
+def merge_data_and_coords(data_vars: DataVars, coords) -> _MergeResult:
     """Used in Dataset.__init__."""
     if isinstance(coords, Coordinates):
         coords = coords.copy()
@@ -666,7 +667,7 @@ class Dataset(
         self,
         # could make a VariableArgs to use more generally, and refine these
         # categories
-        data_vars: Mapping[Any, Any] | None = None,
+        data_vars: DataVars | None = None,
         coords: Mapping[Any, Any] | None = None,
         attrs: Mapping[Any, Any] | None = None,
     ) -> None:
@@ -1220,9 +1221,7 @@ class Dataset(
         else:
             return replaced
 
-    def copy(
-        self, deep: bool = False, data: Mapping[Any, ArrayLike] | None = None
-    ) -> Self:
+    def copy(self, deep: bool = False, data: DataVars | None = None) -> Self:
         """Returns a copy of this dataset.
 
         If `deep=True`, a deep copy is made of each of the component variables.
@@ -1324,7 +1323,7 @@ class Dataset(
     def _copy(
         self,
         deep: bool = False,
-        data: Mapping[Any, ArrayLike] | None = None,
+        data: DataVars | None = None,
         memo: dict[int, Any] | None = None,
     ) -> Self:
         if data is None:
diff --git a/xarray/core/types.py b/xarray/core/types.py
index e9e700b0..6b6f9300 100644
--- a/xarray/core/types.py
+++ b/xarray/core/types.py
@@ -187,6 +187,9 @@ OrderedDims = Union[str, Sequence[Union[Hashable, "ellipsis"]], "ellipsis", None
 T_Chunks = Union[int, dict[Any, Any], Literal["auto"], None]
 T_NormalizedChunks = tuple[tuple[int, ...], ...]
 
+DataVars = Mapping[Any, Any]
+
+
 ErrorOptions = Literal["raise", "ignore"]
 ErrorOptionsWithWarn = Literal["raise", "warn", "ignore"]
 


Test patch:
diff --git a/xarray/tests/test_coding_times.py b/xarray/tests/test_coding_times.py
index 079e432b..5f76a4a2 100644
--- a/xarray/tests/test_coding_times.py
+++ b/xarray/tests/test_coding_times.py
@@ -30,7 +30,7 @@ from xarray.coding.times import (
 from xarray.coding.variables import SerializationWarning
 from xarray.conventions import _update_bounds_attributes, cf_encoder
 from xarray.core.common import contains_cftime_datetimes
-from xarray.testing import assert_allclose, assert_equal, assert_identical
+from xarray.testing import assert_equal, assert_identical
 from xarray.tests import (
     FirstElementAccessibleArray,
     arm_xfail,
@@ -1036,7 +1036,7 @@ def test_encode_cf_datetime_defaults_to_correct_dtype(
         pytest.skip("Nanosecond frequency is not valid for cftime dates.")
     times = date_range("2000", periods=3, freq=freq)
     units = f"{encoding_units} since 2000-01-01"
-    encoded, _, _ = coding.times.encode_cf_datetime(times, units)
+    encoded, _units, _ = coding.times.encode_cf_datetime(times, units)
 
     numpy_timeunit = coding.times._netcdf_to_numpy_timeunit(encoding_units)
     encoding_units_as_timedelta = np.timedelta64(1, numpy_timeunit)
@@ -1212,6 +1212,7 @@ def test_contains_cftime_lazy() -> None:
         ("1677-09-21T00:12:43.145224193", "ns", np.int64, None, False),
         ("1677-09-21T00:12:43.145225", "us", np.int64, None, False),
         ("1970-01-01T00:00:01.000001", "us", np.int64, None, False),
+        ("1677-09-21T00:21:52.901038080", "ns", np.float32, 20.0, True),
     ],
 )
 def test_roundtrip_datetime64_nanosecond_precision(
@@ -1261,14 +1262,52 @@ def test_roundtrip_datetime64_nanosecond_precision_warning() -> None:
     ]
     units = "days since 1970-01-10T01:01:00"
     needed_units = "hours"
-    encoding = dict(_FillValue=20, units=units)
+    new_units = f"{needed_units} since 1970-01-10T01:01:00"
+
+    encoding = dict(dtype=None, _FillValue=20, units=units)
     var = Variable(["time"], times, encoding=encoding)
-    wmsg = (
-        f"Times can't be serialized faithfully with requested units {units!r}. "
-        f"Resolution of {needed_units!r} needed. "
-    )
-    with pytest.warns(UserWarning, match=wmsg):
+    with pytest.warns(UserWarning, match=f"Resolution of {needed_units!r} needed."):
+        encoded_var = conventions.encode_cf_variable(var)
+    assert encoded_var.dtype == np.float64
+    assert encoded_var.attrs["units"] == units
+    assert encoded_var.attrs["_FillValue"] == 20.0
+
+    decoded_var = conventions.decode_cf_variable("foo", encoded_var)
+    assert_identical(var, decoded_var)
+
+    encoding = dict(dtype="int64", _FillValue=20, units=units)
+    var = Variable(["time"], times, encoding=encoding)
+    with pytest.warns(
+        UserWarning, match=f"Serializing with units {new_units!r} instead."
+    ):
+        encoded_var = conventions.encode_cf_variable(var)
+    assert encoded_var.dtype == np.int64
+    assert encoded_var.attrs["units"] == new_units
+    assert encoded_var.attrs["_FillValue"] == 20
+
+    decoded_var = conventions.decode_cf_variable("foo", encoded_var)
+    assert_identical(var, decoded_var)
+
+    encoding = dict(dtype="float64", _FillValue=20, units=units)
+    var = Variable(["time"], times, encoding=encoding)
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
+        encoded_var = conventions.encode_cf_variable(var)
+    assert encoded_var.dtype == np.float64
+    assert encoded_var.attrs["units"] == units
+    assert encoded_var.attrs["_FillValue"] == 20.0
+
+    decoded_var = conventions.decode_cf_variable("foo", encoded_var)
+    assert_identical(var, decoded_var)
+
+    encoding = dict(dtype="int64", _FillValue=20, units=new_units)
+    var = Variable(["time"], times, encoding=encoding)
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
         encoded_var = conventions.encode_cf_variable(var)
+    assert encoded_var.dtype == np.int64
+    assert encoded_var.attrs["units"] == new_units
+    assert encoded_var.attrs["_FillValue"] == 20
 
     decoded_var = conventions.decode_cf_variable("foo", encoded_var)
     assert_identical(var, decoded_var)
@@ -1309,14 +1348,18 @@ def test_roundtrip_timedelta64_nanosecond_precision_warning() -> None:
     needed_units = "hours"
     wmsg = (
         f"Timedeltas can't be serialized faithfully with requested units {units!r}. "
-        f"Resolution of {needed_units!r} needed. "
+        f"Serializing with units {needed_units!r} instead."
     )
-    encoding = dict(_FillValue=20, units=units)
+    encoding = dict(dtype=np.int64, _FillValue=20, units=units)
     var = Variable(["time"], timedelta_values, encoding=encoding)
     with pytest.warns(UserWarning, match=wmsg):
         encoded_var = conventions.encode_cf_variable(var)
+    assert encoded_var.dtype == np.int64
+    assert encoded_var.attrs["units"] == needed_units
+    assert encoded_var.attrs["_FillValue"] == 20
     decoded_var = conventions.decode_cf_variable("foo", encoded_var)
-    assert_allclose(var, decoded_var)
+    assert_identical(var, decoded_var)
+    assert decoded_var.encoding["dtype"] == np.int64
 
 
 def test_roundtrip_float_times() -> None:
diff --git a/xarray/tests/test_computation.py b/xarray/tests/test_computation.py
index b75e80db..87f8328e 100644
--- a/xarray/tests/test_computation.py
+++ b/xarray/tests/test_computation.py
@@ -1190,7 +1190,7 @@ def test_apply_dask() -> None:
 
     # unknown setting for dask array handling
     with pytest.raises(ValueError):
-        apply_ufunc(identity, array, dask="unknown")
+        apply_ufunc(identity, array, dask="unknown")  # type: ignore
 
     def dask_safe_identity(x):
         return apply_ufunc(identity, x, dask="allowed")
diff --git a/xarray/tests/test_dataarray.py b/xarray/tests/test_dataarray.py
index 76dc4345..11ebc4da 100644
--- a/xarray/tests/test_dataarray.py
+++ b/xarray/tests/test_dataarray.py
@@ -38,6 +38,7 @@ from xarray.core.coordinates import Coordinates
 from xarray.core.indexes import Index, PandasIndex, filter_indexes_from_coords
 from xarray.core.types import QueryEngineOptions, QueryParserOptions
 from xarray.core.utils import is_scalar
+from xarray.testing import _assert_internal_invariants
 from xarray.tests import (
     InaccessibleArray,
     ReturnItem,
@@ -415,9 +416,6 @@ class TestDataArray:
         with pytest.raises(ValueError, match=r"conflicting MultiIndex"):
             DataArray(np.random.rand(4, 4), [("x", self.mindex), ("level_1", range(4))])
 
-        with pytest.raises(ValueError, match=r"matching the dimension size"):
-            DataArray(data, coords={"x": 0}, dims=["x", "y"])
-
     def test_constructor_from_self_described(self) -> None:
         data = [[-0.1, 21], [0, 2]]
         expected = DataArray(
@@ -7112,3 +7110,35 @@ class TestStackEllipsis:
         da = DataArray([[1, 2], [1, 2]], dims=("x", "y"))
         with pytest.raises(ValueError):
             da.stack(flat=...)  # type: ignore
+
+
+def test_nD_coord_dataarray() -> None:
+    # should succeed
+    da = DataArray(
+        np.ones((2, 4)),
+        dims=("x", "y"),
+        coords={
+            "x": (("x", "y"), np.arange(8).reshape((2, 4))),
+            "y": ("y", np.arange(4)),
+        },
+    )
+    _assert_internal_invariants(da, check_default_indexes=True)
+
+    da2 = DataArray(np.ones(4), dims=("y"), coords={"y": ("y", np.arange(4))})
+    da3 = DataArray(np.ones(4), dims=("z"))
+
+    _, actual = xr.align(da, da2)
+    assert_identical(da2, actual)
+
+    expected = da.drop_vars("x")
+    _, actual = xr.broadcast(da, da2)
+    assert_identical(expected, actual)
+
+    actual, _ = xr.broadcast(da, da3)
+    expected = da.expand_dims(z=4, axis=-1)
+    assert_identical(actual, expected)
+
+    da4 = DataArray(np.ones((2, 4)), coords={"x": 0}, dims=["x", "y"])
+    _assert_internal_invariants(da4, check_default_indexes=True)
+    assert "x" not in da4.xindexes
+    assert "x" in da4.coords
diff --git a/xarray/tests/test_missing.py b/xarray/tests/test_missing.py
index fe2cdc58..c57d84c9 100644
--- a/xarray/tests/test_missing.py
+++ b/xarray/tests/test_missing.py
@@ -92,26 +92,36 @@ def make_interpolate_example_data(shape, frac_nan, seed=12345, non_uniform=False
     return da, df
 
 
+@pytest.mark.parametrize("fill_value", [None, np.nan, 47.11])
+@pytest.mark.parametrize(
+    "method", ["linear", "nearest", "zero", "slinear", "quadratic", "cubic"]
+)
 @requires_scipy
-def test_interpolate_pd_compat():
+def test_interpolate_pd_compat(method, fill_value) -> None:
     shapes = [(8, 8), (1, 20), (20, 1), (100, 100)]
     frac_nans = [0, 0.5, 1]
-    methods = ["linear", "nearest", "zero", "slinear", "quadratic", "cubic"]
 
-    for shape, frac_nan, method in itertools.product(shapes, frac_nans, methods):
+    for shape, frac_nan in itertools.product(shapes, frac_nans):
         da, df = make_interpolate_example_data(shape, frac_nan)
 
         for dim in ["time", "x"]:
-            actual = da.interpolate_na(method=method, dim=dim, fill_value=np.nan)
+            actual = da.interpolate_na(method=method, dim=dim, fill_value=fill_value)
+            # need limit_direction="both" here, to let pandas fill
+            # in both directions instead of default forward direction only
             expected = df.interpolate(
                 method=method,
                 axis=da.get_axis_num(dim),
+                limit_direction="both",
+                fill_value=fill_value,
             )
-            # Note, Pandas does some odd things with the left/right fill_value
-            # for the linear methods. This next line inforces the xarray
-            # fill_value convention on the pandas output. Therefore, this test
-            # only checks that interpolated values are the same (not nans)
-            expected.values[pd.isnull(actual.values)] = np.nan
+
+            if method == "linear":
+                # Note, Pandas does not take left/right fill_value into account
+                # for the numpy linear methods.
+                # see https://github.com/pandas-dev/pandas/issues/55144
+                # This aligns the pandas output with the xarray output
+                expected.values[pd.isnull(actual.values)] = np.nan
+                expected.values[actual.values == fill_value] = fill_value
 
             np.testing.assert_allclose(actual.values, expected.values)