Cftime arrays not supported by polyval

[aulemahal]

What happened?

I was trying to use polyval with a cftime coordinate and it failed with TypeError: unsupported operand type(s) for *: 'float' and 'cftime._cftime.DatetimeNoLeap'. The error seems to originate from #6548, where the process transforming coordinates to numerical values was modified. The new _ensure_numeric method seems to ignore the possibility of cftime arrays.

What did you expect to happen?

A polynomial to be evaluated along my coordinate.

Minimal Complete Verifiable Example

import xarray as xr
import numpy as np

# use_cftime=False will work
t = xr.date_range('2001-01-01', periods=100, use_cftime=True, freq='YS')
da = xr.DataArray(np.arange(100) ** 3, dims=('time',), coords={'time': t})
coeffs = da.polyfit('time', 4)
da2 = xr.polyval(da.time, coeffs).polyfit_coefficients

MVCE confirmation

• Minimal example — the example is as focused as reasonably possible to demonstrate the underlying issue in xarray.
• Complete example — the example is self-contained, including all data and the text of any traceback.
• Verifiable example — the example copy & pastes into an IPython prompt or Binder notebook, returning the result.
• New issue — a search of GitHub Issues suggests this is not a duplicate.

Relevant log output

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Input In [5], in <cell line: 4>()
      2 da = xr.DataArray(np.arange(100) ** 3, dims=('time',), coords={'time': t})
      3 coeffs = da.polyfit('time', 4)
----> 4 da2 = xr.polyval(da.time, coeffs).polyfit_coefficients

File ~/Python/xarray/xarray/core/computation.py:1931, in polyval(coord, coeffs, degree_dim)
   1929 res = zeros_like(coord) + coeffs.isel({degree_dim: max_deg}, drop=True)
   1930 for deg in range(max_deg - 1, -1, -1):
-> 1931     res *= coord
   1932     res += coeffs.isel({degree_dim: deg}, drop=True)
   1934 return res

File ~/Python/xarray/xarray/core/_typed_ops.py:103, in DatasetOpsMixin.__imul__(self, other)
    102 def __imul__(self, other):
--> 103     return self._inplace_binary_op(other, operator.imul)

File ~/Python/xarray/xarray/core/dataset.py:6107, in Dataset._inplace_binary_op(self, other, f)
   6105     other = other.reindex_like(self, copy=False)
   6106 g = ops.inplace_to_noninplace_op(f)
-> 6107 ds = self._calculate_binary_op(g, other, inplace=True)
   6108 self._replace_with_new_dims(
   6109     ds._variables,
   6110     ds._coord_names,
   (...)
   6113     inplace=True,
   6114 )
   6115 return self

File ~/Python/xarray/xarray/core/dataset.py:6154, in Dataset._calculate_binary_op(self, f, other, join, inplace)
   6152 else:
   6153     other_variable = getattr(other, "variable", other)
-> 6154     new_vars = {k: f(self.variables[k], other_variable) for k in self.data_vars}
   6155 ds._variables.update(new_vars)
   6156 ds._dims = calculate_dimensions(ds._variables)

File ~/Python/xarray/xarray/core/dataset.py:6154, in <dictcomp>(.0)
   6152 else:
   6153     other_variable = getattr(other, "variable", other)
-> 6154     new_vars = {k: f(self.variables[k], other_variable) for k in self.data_vars}
   6155 ds._variables.update(new_vars)
   6156 ds._dims = calculate_dimensions(ds._variables)

File ~/Python/xarray/xarray/core/_typed_ops.py:402, in VariableOpsMixin.__mul__(self, other)
    401 def __mul__(self, other):
--> 402     return self._binary_op(other, operator.mul)

File ~/Python/xarray/xarray/core/variable.py:2494, in Variable._binary_op(self, other, f, reflexive)
   2491 attrs = self._attrs if keep_attrs else None
   2492 with np.errstate(all="ignore"):
   2493     new_data = (
-> 2494         f(self_data, other_data) if not reflexive else f(other_data, self_data)
   2495     )
   2496 result = Variable(dims, new_data, attrs=attrs)
   2497 return result

TypeError: unsupported operand type(s) for *: 'float' and 'cftime._cftime.DatetimeGregorian'

Anything else we need to know?

I also noticed that since the Horner PR, polyfit and polyval do not use the same function to convert coordinates into numerical values. Isn't this dangerous?

Environment

INSTALLED VERSIONS

commit: None
python: 3.10.4 | packaged by conda-forge | (main, Mar 24 2022, 17:38:57) [GCC 10.3.0]
python-bits: 64
OS: Linux
OS-release: 3.10.0-1160.49.1.el7.x86_64
machine: x86_64
processor: x86_64
byteorder: little
LC_ALL: None
LANG: en_CA.UTF-8
LOCALE: ('en_CA', 'UTF-8')
libhdf5: 1.12.1
libnetcdf: 4.8.1

xarray: 2022.3.1.dev267+gd711d58
pandas: 1.4.2
numpy: 1.21.6
scipy: 1.8.0
netCDF4: 1.5.8
pydap: None
h5netcdf: None
h5py: 3.6.0
Nio: None
zarr: 2.11.3
cftime: 1.6.0
nc_time_axis: 1.4.1
PseudoNetCDF: None
rasterio: None
cfgrib: None
iris: None
bottleneck: 1.3.4
dask: 2022.04.1
distributed: 2022.4.1
matplotlib: 3.5.1
cartopy: 0.20.2
seaborn: None
numbagg: None
fsspec: 2022.3.0
cupy: None
pint: 0.19.2
sparse: 0.13.0
flox: 0.5.0
numpy_groupies: 0.9.15
setuptools: 62.1.0
pip: 22.0.4
conda: None
pytest: None
IPython: 8.2.0
sphinx: 4.5.0

[headtr1ck]

I also noticed that since the Horner PR, polyfit and polyval do not use the same function to convert coordinates into numerical values. Isn't this dangerous?

The old function only supported 1D (index) Arrays.

Also, there are currently ideas to reimplement polyfit, so they should probably use the new method as well.