Repository: pydata/xarray

Issue:
Dataset.interp drops boolean variables
<!-- Please include a self-contained copy-pastable example that generates the issue if possible.

Please be concise with code posted. See guidelines below on how to provide a good bug report:

- Craft Minimal Bug Reports: http://matthewrocklin.com/blog/work/2018/02/28/minimal-bug-reports
- Minimal Complete Verifiable Examples: https://stackoverflow.com/help/mcve

Bug reports that follow these guidelines are easier to diagnose, and so are often handled much more quickly.
-->

**What happened**:
`Dataset.interp` silently drops boolean variables.

**What you expected to happen**:
If I'm interpolating a group of variables I expect to get all of them back in the correct shape with relevant values in them.

If the variables are boolean or object arrays I don't expect it to do linear interpolation because it doesn't make sense but stepwise interpolation like nearest or zero order interpolation should be fine to expect.

**Minimal Complete Verifiable Example**:

```python
import numpy as np
a = np.arange(0, 5)
b = np.core.defchararray.add("long_variable_name", a.astype(str))
coords = dict(time=da.array([0, 1]))
data_vars = dict()
for v in b:
    data_vars[v] = xr.DataArray(
        name=v,
        data=np.array([0, 1]).astype(bool),
        dims=["time"],
        coords=coords,
    )
ds1 = xr.Dataset(data_vars)

# Print raw data:
print(ds1)
Out[3]: 
<xarray.Dataset>
Dimensions:              (time: 2)
Coordinates:
  * time                 (time) int32 0 1
Data variables:
    long_variable_name0  (time) bool False True
    long_variable_name1  (time) bool False True
    long_variable_name2  (time) bool False True
    long_variable_name3  (time) bool False True
    long_variable_name4  (time) bool False True

# Interpolate:
ds1 = ds1.interp(
    time=da.array([0, 0.5, 1, 2]),
    assume_sorted=True,
    method="nearest",
    kwargs=dict(fill_value="extrapolate"),
)

# Print interpolated data:
<xarray.Dataset>
Dimensions:  (time: 4)
Coordinates:
  * time     (time) float64 0.0 0.5 1.0 2.0
Data variables:
    *empty*
```

**Anything else we need to know?**:
`ds.interp_like `use `ds.reindex` in these cases which seems like a good choice in `ds.interp` as well. But I think that both `ds.interp` and `ds.interp_like` should fill by default with nearest value instead of np.nan because we're still requesting interpolation.

**Environment**:

<details><summary>Output of <tt>xr.show_versions()</tt></summary>

xr.show_versions()

INSTALLED VERSIONS
------------------
commit: None
python: 3.8.5 (default, Sep  3 2020, 21:29:08) [MSC v.1916 64 bit (AMD64)]
python-bits: 64
OS: Windows
libhdf5: 1.10.4
libnetcdf: None

xarray: 0.16.2
pandas: 1.1.5
numpy: 1.17.5
scipy: 1.4.1
netCDF4: None
pydap: None
h5netcdf: None
h5py: 2.10.0
Nio: None
zarr: None
cftime: None
nc_time_axis: None
PseudoNetCDF: None
rasterio: None
cfgrib: None
iris: None
bottleneck: 1.3.2
dask: 2020.12.0
distributed: 2020.12.0
matplotlib: 3.3.2
cartopy: None
seaborn: 0.11.1
numbagg: None
pint: None
setuptools: 51.0.0.post20201207
pip: 20.3.3
conda: 4.9.2
pytest: 6.2.1
IPython: 7.19.0
sphinx: 3.4.0


</details>


Discussion:
pep8speaks: Hello @Illviljan! Thanks for opening 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/core/dataset.py`](https://github.com/pydata/xarray/blob/ca2ef3a461400e92774c66ef65477b83019dd88e/xarray/core/dataset.py):

> [Line 2895:36](https://github.com/pydata/xarray/blob/ca2ef3a461400e92774c66ef65477b83019dd88e/xarray/core/dataset.py#L2895): [E999](https://duckduckgo.com/?q=pep8%20E999) SyntaxError: invalid syntax
> [Line 2896:21](https://github.com/pydata/xarray/blob/ca2ef3a461400e92774c66ef65477b83019dd88e/xarray/core/dataset.py#L2896): [E116](https://duckduckgo.com/?q=pep8%20E116) unexpected indentation (comment)
> [Line 2897:21](https://github.com/pydata/xarray/blob/ca2ef3a461400e92774c66ef65477b83019dd88e/xarray/core/dataset.py#L2897): [E116](https://duckduckgo.com/?q=pep8%20E116) unexpected indentation (comment)
> [Line 2898:21](https://github.com/pydata/xarray/blob/ca2ef3a461400e92774c66ef65477b83019dd88e/xarray/core/dataset.py#L2898): [E113](https://duckduckgo.com/?q=pep8%20E113) unexpected indentation


dcherian: This seems like a good change but shouldn't we do this only when `method="nearest"`?
dcherian: > How do you linearly interpolate strings, objects or booleans?

This is why we drop them now. my point is that doing "nearest" when the user has specified `method="linear"` is v. confusing. Perhaps we should have `method="linear+nearest"` (or something better) for what you propose: linear for numeric, nearest for everything else. 

Q: Does the current version of xarray just ignore strings, objects & booleans when `method="nearest"`?
Illviljan: Well, as a user I was very confused why my boolean variables silently disappeared. 
And I don't get why `ds.interp` has any mandate to drop variables, if it can't handle interpolating a variable with the method an error should've been raised.
But I'm a lazy user too so I of course would prefer to not have to deal with errors or figure out which variables were dropped and do multiple interpolation steps and finally merge the different datasets back together again. This is why I prefer ds.interp to do the only current method that makes sense for non-numerics rather than dropping them.

`ds.interp` currently only interpolates `if var.dtype.kind in "uifc"`, so it ignores bools, strings, and objects for any method. Which made sense I think because the scipy interpolants always attempts to force to float. This is why I use `reindex_variables` as it both retains the dtype and handles strings/objects just fine.

I've added a parameter now `method_for_non_numerics` for these cases now, so now you can choose which reindex method you want.
max-sixty: Let's discuss this at the next dev meeting.
Illviljan: @dcherian, renamed the parameter and fixed some errors due to the #5102 merge.

Docs failure seems unrelated.
dcherian: Thanks @Illviljan ! Great PR. Thank you for your patience

Patch:
diff --git a/xarray/core/dataset.py b/xarray/core/dataset.py
index 706ccbde..19af7f6c 100644
--- a/xarray/core/dataset.py
+++ b/xarray/core/dataset.py
@@ -2853,6 +2853,7 @@ class Dataset(DataWithCoords, DatasetArithmetic, Mapping):
         method: str = "linear",
         assume_sorted: bool = False,
         kwargs: Mapping[str, Any] = None,
+        method_non_numeric: str = "nearest",
         **coords_kwargs: Any,
     ) -> "Dataset":
         """Multidimensional interpolation of Dataset.
@@ -2877,6 +2878,9 @@ class Dataset(DataWithCoords, DatasetArithmetic, Mapping):
             Additional keyword arguments passed to scipy's interpolator. Valid
             options and their behavior depend on if 1-dimensional or
             multi-dimensional interpolation is used.
+        method_non_numeric : {"nearest", "pad", "ffill", "backfill", "bfill"}, optional
+            Method for non-numeric types. Passed on to :py:meth:`Dataset.reindex`.
+            ``"nearest"`` is used by default.
         **coords_kwargs : {dim: coordinate, ...}, optional
             The keyword arguments form of ``coords``.
             One of coords or coords_kwargs must be provided.
@@ -3034,6 +3038,7 @@ class Dataset(DataWithCoords, DatasetArithmetic, Mapping):
             }
 
         variables: Dict[Hashable, Variable] = {}
+        to_reindex: Dict[Hashable, Variable] = {}
         for name, var in obj._variables.items():
             if name in indexers:
                 continue
@@ -3043,20 +3048,45 @@ class Dataset(DataWithCoords, DatasetArithmetic, Mapping):
             else:
                 use_indexers = validated_indexers
 
-            if var.dtype.kind in "uifc":
+            dtype_kind = var.dtype.kind
+            if dtype_kind in "uifc":
+                # For normal number types do the interpolation:
                 var_indexers = {k: v for k, v in use_indexers.items() if k in var.dims}
                 variables[name] = missing.interp(var, var_indexers, method, **kwargs)
+            elif dtype_kind in "ObU" and (use_indexers.keys() & var.dims):
+                # For types that we do not understand do stepwise
+                # interpolation to avoid modifying the elements.
+                # Use reindex_variables instead because it supports
+                # booleans and objects and retains the dtype but inside
+                # this loop there might be some duplicate code that slows it
+                # down, therefore collect these signals and run it later:
+                to_reindex[name] = var
             elif all(d not in indexers for d in var.dims):
-                # keep unrelated object array
+                # For anything else we can only keep variables if they
+                # are not dependent on any coords that are being
+                # interpolated along:
                 variables[name] = var
 
+        if to_reindex:
+            # Reindex variables:
+            variables_reindex = alignment.reindex_variables(
+                variables=to_reindex,
+                sizes=obj.sizes,
+                indexes=obj.xindexes,
+                indexers={k: v[-1] for k, v in validated_indexers.items()},
+                method=method_non_numeric,
+            )[0]
+            variables.update(variables_reindex)
+
+        # Get the coords that also exist in the variables:
         coord_names = obj._coord_names & variables.keys()
+        # Get the indexes that are not being interpolated along:
         indexes = {k: v for k, v in obj.xindexes.items() if k not in indexers}
         selected = self._replace_with_new_dims(
             variables.copy(), coord_names, indexes=indexes
         )
 
-        # attach indexer as coordinate
+        # Attach indexer as coordinate
         variables.update(indexers)
         for k, v in indexers.items():
             assert isinstance(v, Variable)
@@ -3077,6 +3107,7 @@ class Dataset(DataWithCoords, DatasetArithmetic, Mapping):
         method: str = "linear",
         assume_sorted: bool = False,
         kwargs: Mapping[str, Any] = None,
+        method_non_numeric: str = "nearest",
     ) -> "Dataset":
         """Interpolate this object onto the coordinates of another object,
         filling the out of range values with NaN.
@@ -3098,6 +3129,9 @@ class Dataset(DataWithCoords, DatasetArithmetic, Mapping):
             values.
         kwargs : dict, optional
             Additional keyword passed to scipy's interpolator.
+        method_non_numeric : {"nearest", "pad", "ffill", "backfill", "bfill"}, optional
+            Method for non-numeric types. Passed on to :py:meth:`Dataset.reindex`.
+            ``"nearest"`` is used by default.
 
         Returns
         -------
@@ -3133,7 +3167,13 @@ class Dataset(DataWithCoords, DatasetArithmetic, Mapping):
             # We do not support interpolation along object coordinate.
             # reindex instead.
             ds = self.reindex(object_coords)
-        return ds.interp(numeric_coords, method, assume_sorted, kwargs)
+        return ds.interp(
+            coords=numeric_coords,
+            method=method,
+            assume_sorted=assume_sorted,
+            kwargs=kwargs,
+            method_non_numeric=method_non_numeric,
+        )
 
     # Helper methods for rename()
     def _rename_vars(self, name_dict, dims_dict):


Test patch:
diff --git a/xarray/tests/test_interp.py b/xarray/tests/test_interp.py
index 97a6d236..ab023dc1 100644
--- a/xarray/tests/test_interp.py
+++ b/xarray/tests/test_interp.py
@@ -416,15 +416,19 @@ def test_errors(use_dask):
 
 @requires_scipy
 def test_dtype():
-    ds = xr.Dataset(
-        {"var1": ("x", [0, 1, 2]), "var2": ("x", ["a", "b", "c"])},
-        coords={"x": [0.1, 0.2, 0.3], "z": ("x", ["a", "b", "c"])},
-    )
-    actual = ds.interp(x=[0.15, 0.25])
-    assert "var1" in actual
-    assert "var2" not in actual
-    # object array should be dropped
-    assert "z" not in actual.coords
+    data_vars = dict(
+        a=("time", np.array([1, 1.25, 2])),
+        b=("time", np.array([True, True, False], dtype=bool)),
+        c=("time", np.array(["start", "start", "end"], dtype=str)),
+    )
+    time = np.array([0, 0.25, 1], dtype=float)
+    expected = xr.Dataset(data_vars, coords=dict(time=time))
+    actual = xr.Dataset(
+        {k: (dim, arr[[0, -1]]) for k, (dim, arr) in data_vars.items()},
+        coords=dict(time=time[[0, -1]]),
+    )
+    actual = actual.interp(time=time, method="linear")
+    assert_identical(expected, actual)
 
 
 @requires_scipy