Repository: pydata/xarray

Issue:
datetime interpolation doesn't work
#### Code Sample, a copy-pastable example if possible

This code doesn't work anymore on master.

```python
a = xr.DataArray(np.arange(21).reshape(3, 7), 
                 dims=['x', 'time'], 
                 coords={'x': [1, 2, 3], 'time': pd.date_range('01-01-2001', periods=7, freq='D')})
xi = xr.DataArray(np.linspace(1, 3, 50),  
                  dims=['time'], 
                  coords={'time': pd.date_range('01-01-2001', periods=50, freq='H')})
a.interp(x=xi, time=xi.time)
```

#### Problem description

The above code now raises the error

```
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-26-dda3a6d5725b> in <module>
      6                   dims=['time'],
      7                   coords={'time': pd.date_range('01-01-2001', periods=50, freq='H')})
----> 8 a.interp(x=xi, time=xi.time)

~/work/python/xarray/xarray/core/dataarray.py in interp(self, coords, method, assume_sorted, kwargs, **coords_kwargs)
   1032         ds = self._to_temp_dataset().interp(
   1033             coords, method=method, kwargs=kwargs, assume_sorted=assume_sorted,
-> 1034             **coords_kwargs)
   1035         return self._from_temp_dataset(ds)
   1036 

~/work/python/xarray/xarray/core/dataset.py in interp(self, coords, method, assume_sorted, kwargs, **coords_kwargs)
   2008                         in indexers.items() if k in var.dims}
   2009                     variables[name] = missing.interp(
-> 2010                         var, var_indexers, method, **kwargs)
   2011                 elif all(d not in indexers for d in var.dims):
   2012                     # keep unrelated object array

~/work/python/xarray/xarray/core/missing.py in interp(var, indexes_coords, method, **kwargs)
    468     new_dims = broadcast_dims + list(destination[0].dims)
    469     interped = interp_func(var.transpose(*original_dims).data,
--> 470                            x, destination, method, kwargs)
    471 
    472     result = Variable(new_dims, interped, attrs=var.attrs)

~/work/python/xarray/xarray/core/missing.py in interp_func(var, x, new_x, method, kwargs)
    535                              new_axis=new_axis, drop_axis=drop_axis)
    536 
--> 537     return _interpnd(var, x, new_x, func, kwargs)
    538 
    539 

~/work/python/xarray/xarray/core/missing.py in _interpnd(var, x, new_x, func, kwargs)
    558     var = var.transpose(range(-len(x), var.ndim - len(x)))
    559     # stack new_x to 1 vector, with reshape
--> 560     xi = np.stack([x1.values.ravel() for x1 in new_x], axis=-1)
    561     rslt = func(x, var, xi, **kwargs)
    562     # move back the interpolation axes to the last position

~/work/python/xarray/xarray/core/missing.py in <listcomp>(.0)
    558     var = var.transpose(range(-len(x), var.ndim - len(x)))
    559     # stack new_x to 1 vector, with reshape
--> 560     xi = np.stack([x1.values.ravel() for x1 in new_x], axis=-1)
    561     rslt = func(x, var, xi, **kwargs)
    562     # move back the interpolation axes to the last position

AttributeError: 'numpy.ndarray' object has no attribute 'values'
```

I think the issue is this line which returns a numpy array instead of a Variable. This was added in the `coarsen` PR (cc @fujiisoup)
https://github.com/pydata/xarray/blob/d4c46829b283ab7e7b7db8b86dae77861ce68f3c/xarray/core/utils.py#L636

Discussion:
shoyer: @spencerkclark what convention do you use for missing values in cftime arrays?
spencerkclark: > what convention do you use for missing values in cftime arrays?

Currently there is nothing built in to xarray to support missing values for cftime arrays.  I think we would have to come up with a convention.
shoyer: OK, in that case let's only worry about NaT for datetime64/timedelta64
arrays.

On Fri, Jan 11, 2019 at 2:37 PM Spencer Clark <notifications@github.com>
wrote:

> what convention do you use for missing values in cftime arrays?
>
> Currently there is nothing built in to xarray to support missing values
> for cftime arrays. I think we would have to come up with a convention.
>
> —
> You are receiving this because you commented.
> Reply to this email directly, view it on GitHub
> <https://github.com/pydata/xarray/pull/2668#issuecomment-453679270>, or mute
> the thread
> <https://github.com/notifications/unsubscribe-auth/ABKS1iSjXGAEXm6e5gqGCY390bNdMIiMks5vCRISgaJpZM4Z8P5F>
> .
>

shoyer: (at least for now, until we can figure out how to handle missing values)

On Fri, Jan 11, 2019 at 2:44 PM Stephan Hoyer <shoyer@gmail.com> wrote:

> OK, in that case let's only worry about NaT for datetime64/timedelta64
> arrays.
>
> On Fri, Jan 11, 2019 at 2:37 PM Spencer Clark <notifications@github.com>
> wrote:
>
>> what convention do you use for missing values in cftime arrays?
>>
>> Currently there is nothing built in to xarray to support missing values
>> for cftime arrays. I think we would have to come up with a convention.
>>
>> —
>> You are receiving this because you commented.
>> Reply to this email directly, view it on GitHub
>> <https://github.com/pydata/xarray/pull/2668#issuecomment-453679270>, or mute
>> the thread
>> <https://github.com/notifications/unsubscribe-auth/ABKS1iSjXGAEXm6e5gqGCY390bNdMIiMks5vCRISgaJpZM4Z8P5F>
>> .
>>
>

pep8speaks: Hello @fujiisoup! Thanks for updating the PR.











 - In the file [`xarray/tests/test_duck_array_ops.py`](https://github.com/pydata/xarray/blob/fca735cb9a9a1e0d11e7a202f39fc363451953ce/xarray/tests/test_duck_array_ops.py), following are the PEP8 issues :

> [Line 607:80](https://github.com/pydata/xarray/blob/fca735cb9a9a1e0d11e7a202f39fc363451953ce/xarray/tests/test_duck_array_ops.py#L607): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (85 > 79 characters)
> [Line 612:80](https://github.com/pydata/xarray/blob/fca735cb9a9a1e0d11e7a202f39fc363451953ce/xarray/tests/test_duck_array_ops.py#L612): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (83 > 79 characters)






fujiisoup: I couldn't make `datetime_to_numeric` compatible both `np.ndarray` and `Variable`. I think it would be cleaner if `Variable` has a (private) method to do this, like `Variable._to_numeric(offset=None, datetime_unit=None)`.
fujiisoup: Thanks, @spencerkclark 

It looks I need more work for this (probably a week). 
It is not working for cftime stuffs, but I do not yet understand how it is handled in xarray.
spencerkclark: > It looks I need more work for this (probably a week).
> It is not working for cftime stuffs, but I do not yet understand how it is handled in xarray.

Thanks @fujiisoup, no worries.  There was only one place where I saw a potential issue for cftime dates, which I noted above.  I made a PR to your branch with the fix (https://github.com/fujiisoup/xarray/pull/9); let me know if there were other places with issues and I can sort those out too.
fujiisoup: Thanks, @spencerkclark, for the PR.
It was super helpful.

I will finish it up tonight (hopefully).
fujiisoup: @spencerkclark 

Sorry for leaving it for a long time.
I was a little stressed these days...

> Should we strive to have compatibility in `datetime_to_numeric` there? 

As I'm not a heavy user of `datetime` arrays, I personally do not mind if we do not have dask compatibility for `datetime_to_numeric`.
Probably it would be OK to leave it for the later development if someone wants this feature.
spencerkclark: > As I'm not a heavy user of `datetime` arrays, I personally do not mind if we do not have dask compatibility for `datetime_to_numeric`.

No worries @fujiisoup, this is totally fine with me.  I think the most common use case for datetime arrays is in indexes, which currently need to be stored in memory anyway, so I don't think the dask use case would come up very often.

Patch:
diff --git a/xarray/core/dataset.py b/xarray/core/dataset.py
index 636cff2c..7bb08584 100644
--- a/xarray/core/dataset.py
+++ b/xarray/core/dataset.py
@@ -28,6 +28,7 @@ from .coordinates import (
     DatasetCoordinates, LevelCoordinatesSource, assert_coordinate_consistent,
     remap_label_indexers,
 )
+from .duck_array_ops import datetime_to_numeric
 from .indexes import Indexes, default_indexes, isel_variable_and_index
 from .merge import (
     dataset_merge_method, dataset_update_method, merge_data_and_coords,
@@ -35,9 +36,9 @@ from .merge import (
 from .options import OPTIONS, _get_keep_attrs
 from .pycompat import dask_array_type
 from .utils import (
-    Frozen, SortedKeysDict, _check_inplace, datetime_to_numeric,
-    decode_numpy_dict_values, either_dict_or_kwargs, hashable,
-    maybe_wrap_array)
+    Frozen, SortedKeysDict, _check_inplace,
+    decode_numpy_dict_values, either_dict_or_kwargs, ensure_us_time_resolution,
+    hashable, maybe_wrap_array)
 from .variable import IndexVariable, Variable, as_variable, broadcast_variables
 if TYPE_CHECKING:
     from .dataarray import DataArray
@@ -3997,15 +3998,14 @@ class Dataset(Mapping, ImplementsDatasetReduce, DataWithCoords):
                 datetime_unit, _ = np.datetime_data(coord_var.dtype)
             elif datetime_unit is None:
                 datetime_unit = 's'  # Default to seconds for cftime objects
-            coord_var = datetime_to_numeric(
-                coord_var, datetime_unit=datetime_unit)
+            coord_var = coord_var._to_numeric(datetime_unit=datetime_unit)
 
         variables = OrderedDict()
         for k, v in self.variables.items():
             if (k in self.data_vars and dim in v.dims and
                     k not in self.coords):
                 if _contains_datetime_like_objects(v):
-                    v = datetime_to_numeric(v, datetime_unit=datetime_unit)
+                    v = v._to_numeric(datetime_unit=datetime_unit)
                 grad = duck_array_ops.gradient(
                     v.data, coord_var, edge_order=edge_order,
                     axis=v.get_axis_num(dim))
diff --git a/xarray/core/duck_array_ops.py b/xarray/core/duck_array_ops.py
index 36c40902..4d6d716a 100644
--- a/xarray/core/duck_array_ops.py
+++ b/xarray/core/duck_array_ops.py
@@ -11,7 +11,7 @@ import warnings
 import numpy as np
 import pandas as pd
 
-from . import dask_array_ops, dtypes, npcompat, nputils, utils
+from . import dask_array_ops, dtypes, npcompat, nputils
 from .nputils import nanfirst, nanlast
 from .pycompat import dask_array_type
 
@@ -289,14 +289,51 @@ cumsum_1d.numeric_only = True
 _mean = _create_nan_agg_method('mean')
 
 
+def datetime_to_numeric(array, offset=None, datetime_unit=None, dtype=float):
+    """Convert an array containing datetime-like data to an array of floats.
+
+    Parameters
+    ----------
+    da : array
+        Input data
+    offset: Scalar with the same type of array or None
+        If None, subtract minimum values to reduce round off error
+    datetime_unit: None or any of {'Y', 'M', 'W', 'D', 'h', 'm', 's', 'ms',
+        'us', 'ns', 'ps', 'fs', 'as'}
+    dtype: target dtype
+
+    Returns
+    -------
+    array
+    """
+    if offset is None:
+        offset = array.min()
+    array = array - offset
+
+    if not hasattr(array, 'dtype'):  # scalar is converted to 0d-array
+        array = np.array(array)
+
+    if array.dtype.kind in 'O':
+        # possibly convert object array containing datetime.timedelta
+        array = np.asarray(pd.Series(array.ravel())).reshape(array.shape)
+
+    if datetime_unit:
+        array = array / np.timedelta64(1, datetime_unit)
+
+    # convert np.NaT to np.nan
+    if array.dtype.kind in 'mM':
+        return np.where(isnull(array), np.nan, array.astype(dtype))
+    return array.astype(dtype)
+
+
 def mean(array, axis=None, skipna=None, **kwargs):
     """ inhouse mean that can handle datatime dtype """
     array = asarray(array)
-    if array.dtype.kind == 'M':
+    if array.dtype.kind in 'Mm':
         offset = min(array)
         # xarray always uses datetime[ns] for datetime
         dtype = 'timedelta64[ns]'
-        return _mean(utils.datetime_to_numeric(array, offset), axis=axis,
+        return _mean(datetime_to_numeric(array, offset), axis=axis,
                      skipna=skipna, **kwargs).astype(dtype) + offset
     else:
         return _mean(array, axis=axis, skipna=skipna, **kwargs)
diff --git a/xarray/core/missing.py b/xarray/core/missing.py
index 4c9435e0..50c42020 100644
--- a/xarray/core/missing.py
+++ b/xarray/core/missing.py
@@ -9,8 +9,8 @@ import pandas as pd
 from . import utils
 from .common import _contains_datetime_like_objects
 from .computation import apply_ufunc
-from .duck_array_ops import dask_array_type
-from .utils import OrderedSet, datetime_to_numeric, is_scalar
+from .duck_array_ops import dask_array_type, datetime_to_numeric
+from .utils import OrderedSet, is_scalar
 from .variable import Variable, broadcast_variables
 
 
@@ -411,10 +411,9 @@ def _floatize_x(x, new_x):
             # We assume that the most of the bits are used to represent the
             # offset (min(x)) and the variation (x - min(x)) can be
             # represented by float.
-            xmin = x[i].min()
-            x[i] = datetime_to_numeric(x[i], offset=xmin, dtype=np.float64)
-            new_x[i] = datetime_to_numeric(
-                new_x[i], offset=xmin, dtype=np.float64)
+            xmin = x[i].values.min()
+            x[i] = x[i]._to_numeric(offset=xmin, dtype=np.float64)
+            new_x[i] = new_x[i]._to_numeric(offset=xmin, dtype=np.float64)
     return x, new_x
 
 
diff --git a/xarray/core/utils.py b/xarray/core/utils.py
index b8e81869..053a45f0 100644
--- a/xarray/core/utils.py
+++ b/xarray/core/utils.py
@@ -603,40 +603,6 @@ class HiddenKeyDict(MutableMapping):
         return len(self._data) - num_hidden
 
 
-def datetime_to_numeric(array, offset=None, datetime_unit=None, dtype=float):
-    """Convert an array containing datetime-like data to an array of floats.
-
-    Parameters
-    ----------
-    da : array
-        Input data
-    offset: Scalar with the same type of array or None
-        If None, subtract minimum values to reduce round off error
-    datetime_unit: None or any of {'Y', 'M', 'W', 'D', 'h', 'm', 's', 'ms',
-        'us', 'ns', 'ps', 'fs', 'as'}
-    dtype: target dtype
-
-    Returns
-    -------
-    array
-    """
-    from . import duck_array_ops
-
-    if offset is None:
-        offset = array.min()
-    array = array - offset
-
-    if datetime_unit:
-        array = array / np.timedelta64(1, datetime_unit)
-    # convert np.NaT to np.nan
-    if array.dtype.kind in 'mM':
-        if hasattr(array, 'isnull'):
-            return np.where(array.isnull(), np.nan, array.astype(dtype))
-        return np.where(duck_array_ops.isnull(array), np.nan,
-                        array.astype(dtype))
-    return array
-
-
 def get_temp_dimname(dims, new_dim):
     """ Get an new dimension name based on new_dim, that is not used in dims.
     If the same name exists, we add an underscore(s) in the head.
diff --git a/xarray/core/variable.py b/xarray/core/variable.py
index a35f8cf0..85eab294 100644
--- a/xarray/core/variable.py
+++ b/xarray/core/variable.py
@@ -1782,6 +1782,14 @@ class Variable(common.AbstractArray, arithmetic.SupportsArithmetic,
             return self
         return func
 
+    def _to_numeric(self, offset=None, datetime_unit=None, dtype=float):
+        """ A (private) method to convert datetime array to numeric dtype
+        See duck_array_ops.datetime_to_numeric
+        """
+        numeric_array = duck_array_ops.datetime_to_numeric(
+            self.data, offset, datetime_unit, dtype)
+        return type(self)(self.dims, numeric_array, self._attrs)
+
 
 ops.inject_all_ops_and_reduce_methods(Variable)
 


Test patch:
diff --git a/xarray/tests/test_dataset.py b/xarray/tests/test_dataset.py
index 463c6756..e4ffdad4 100644
--- a/xarray/tests/test_dataset.py
+++ b/xarray/tests/test_dataset.py
@@ -16,7 +16,7 @@ from xarray import (
     ALL_DIMS, DataArray, Dataset, IndexVariable, MergeError, Variable, align,
     backends, broadcast, open_dataset, set_options)
 from xarray.core import dtypes, indexing, npcompat, utils
-from xarray.core.common import full_like
+from xarray.core.common import duck_array_ops, full_like
 from xarray.core.pycompat import integer_types
 
 from . import (
@@ -4676,7 +4676,7 @@ def test_differentiate_datetime(dask):
     actual = da.differentiate('x', edge_order=1, datetime_unit='D')
     expected_x = xr.DataArray(
         npcompat.gradient(
-            da, utils.datetime_to_numeric(da['x'], datetime_unit='D'),
+            da, da['x'].variable._to_numeric(datetime_unit='D'),
             axis=0, edge_order=1), dims=da.dims, coords=da.coords)
     assert_equal(expected_x, actual)
 
@@ -4710,7 +4710,7 @@ def test_differentiate_cftime(dask):
 
     actual = da.differentiate('time', edge_order=1, datetime_unit='D')
     expected_data = npcompat.gradient(
-        da, utils.datetime_to_numeric(da['time'], datetime_unit='D'),
+        da, da['time'].variable._to_numeric(datetime_unit='D'),
         axis=0, edge_order=1)
     expected = xr.DataArray(expected_data, coords=da.coords, dims=da.dims)
     assert_equal(expected, actual)
@@ -4789,7 +4789,8 @@ def test_trapz_datetime(dask, which_datetime):
 
     actual = da.integrate('time', datetime_unit='D')
     expected_data = np.trapz(
-        da, utils.datetime_to_numeric(da['time'], datetime_unit='D'), axis=0)
+        da, duck_array_ops.datetime_to_numeric(da['time'], datetime_unit='D'),
+        axis=0)
     expected = xr.DataArray(
         expected_data, dims=['y'],
         coords={k: v for k, v in da.coords.items() if 'time' not in v.dims})
diff --git a/xarray/tests/test_duck_array_ops.py b/xarray/tests/test_duck_array_ops.py
index ba7f6ba5..ab3cafed 100644
--- a/xarray/tests/test_duck_array_ops.py
+++ b/xarray/tests/test_duck_array_ops.py
@@ -7,16 +7,17 @@ import pandas as pd
 import pytest
 from numpy import array, nan
 
-from xarray import DataArray, Dataset, concat
+from xarray import DataArray, Dataset, concat, cftime_range
 from xarray.core import dtypes, duck_array_ops
 from xarray.core.duck_array_ops import (
     array_notnull_equiv, concatenate, count, first, gradient, last, mean,
     rolling_window, stack, where)
 from xarray.core.pycompat import dask_array_type
-from xarray.testing import assert_allclose, assert_equal
+from xarray.testing import assert_allclose, assert_equal, assert_identical
 
 from . import (
-    assert_array_equal, has_dask, has_np113, raises_regex, requires_dask)
+    assert_array_equal, has_dask, has_np113, raises_regex, requires_cftime,
+    requires_dask)
 
 
 class TestOps(object):
@@ -569,3 +570,42 @@ def test_docs():
             indicated dimension(s) removed.
         """)
     assert actual == expected
+
+
+def test_datetime_to_numeric_datetime64():
+    times = pd.date_range('2000', periods=5, freq='7D').values
+    result = duck_array_ops.datetime_to_numeric(times, datetime_unit='h')
+    expected = 24 * np.arange(0, 35, 7)
+    np.testing.assert_array_equal(result, expected)
+
+    offset = times[1]
+    result = duck_array_ops.datetime_to_numeric(
+        times, offset=offset, datetime_unit='h')
+    expected = 24 * np.arange(-7, 28, 7)
+    np.testing.assert_array_equal(result, expected)
+
+    dtype = np.float32
+    result = duck_array_ops.datetime_to_numeric(
+        times, datetime_unit='h', dtype=dtype)
+    expected = 24 * np.arange(0, 35, 7).astype(dtype)
+    np.testing.assert_array_equal(result, expected)
+
+
+@requires_cftime
+def test_datetime_to_numeric_cftime():
+    times = cftime_range('2000', periods=5, freq='7D').values
+    result = duck_array_ops.datetime_to_numeric(times, datetime_unit='h')
+    expected = 24 * np.arange(0, 35, 7)
+    np.testing.assert_array_equal(result, expected)
+
+    offset = times[1]
+    result = duck_array_ops.datetime_to_numeric(
+        times, offset=offset, datetime_unit='h')
+    expected = 24 * np.arange(-7, 28, 7)
+    np.testing.assert_array_equal(result, expected)
+
+    dtype = np.float32
+    result = duck_array_ops.datetime_to_numeric(
+        times, datetime_unit='h', dtype=dtype)
+    expected = 24 * np.arange(0, 35, 7).astype(dtype)
+    np.testing.assert_array_equal(result, expected)
diff --git a/xarray/tests/test_interp.py b/xarray/tests/test_interp.py
index d01929f1..0d92f937 100644
--- a/xarray/tests/test_interp.py
+++ b/xarray/tests/test_interp.py
@@ -571,3 +571,16 @@ def test_cftime_to_non_cftime_error():
 
     with pytest.raises(TypeError):
         da.interp(time=0.5)
+
+
+@requires_scipy
+def test_datetime_interp_noerror():
+    # GH:2667
+    a = xr.DataArray(
+        np.arange(21).reshape(3, 7), dims=['x', 'time'],
+        coords={'x': [1, 2, 3],
+                'time': pd.date_range('01-01-2001', periods=7, freq='D')})
+    xi = xr.DataArray(
+        np.linspace(1, 3, 50), dims=['time'],
+        coords={'time': pd.date_range('01-01-2001', periods=50, freq='H')})
+    a.interp(x=xi, time=xi.time)  # should not raise an error
diff --git a/xarray/tests/test_utils.py b/xarray/tests/test_utils.py
index 09152bac..e98ab5cd 100644
--- a/xarray/tests/test_utils.py
+++ b/xarray/tests/test_utils.py
@@ -279,42 +279,3 @@ def test_either_dict_or_kwargs():
 
     with pytest.raises(ValueError, match=r'foo'):
         result = either_dict_or_kwargs(dict(a=1), dict(a=1), 'foo')
-
-
-def test_datetime_to_numeric_datetime64():
-    times = pd.date_range('2000', periods=5, freq='7D')
-    da = xr.DataArray(times, coords=[times], dims=['time'])
-    result = utils.datetime_to_numeric(da, datetime_unit='h')
-    expected = 24 * xr.DataArray(np.arange(0, 35, 7), coords=da.coords)
-    assert_identical(result, expected)
-
-    offset = da.isel(time=1)
-    result = utils.datetime_to_numeric(da, offset=offset, datetime_unit='h')
-    expected = 24 * xr.DataArray(np.arange(-7, 28, 7), coords=da.coords)
-    assert_identical(result, expected)
-
-    dtype = np.float32
-    result = utils.datetime_to_numeric(da, datetime_unit='h', dtype=dtype)
-    expected = 24 * xr.DataArray(
-        np.arange(0, 35, 7), coords=da.coords).astype(dtype)
-    assert_identical(result, expected)
-
-
-@requires_cftime
-def test_datetime_to_numeric_cftime():
-    times = xr.cftime_range('2000', periods=5, freq='7D')
-    da = xr.DataArray(times, coords=[times], dims=['time'])
-    result = utils.datetime_to_numeric(da, datetime_unit='h')
-    expected = 24 * xr.DataArray(np.arange(0, 35, 7), coords=da.coords)
-    assert_identical(result, expected)
-
-    offset = da.isel(time=1)
-    result = utils.datetime_to_numeric(da, offset=offset, datetime_unit='h')
-    expected = 24 * xr.DataArray(np.arange(-7, 28, 7), coords=da.coords)
-    assert_identical(result, expected)
-
-    dtype = np.float32
-    result = utils.datetime_to_numeric(da, datetime_unit='h', dtype=dtype)
-    expected = 24 * xr.DataArray(
-        np.arange(0, 35, 7), coords=da.coords).astype(dtype)
-    assert_identical(result, expected)