Repository: pydata/xarray

Issue:
Temporal polyfit does not use the same coordinate as polyval
### What happened?

When fitting a variable along a temporal dimension, results from computing the fitted "trend" with  `xr.polyval` were very different from the original variable. It seems that the function is not evaluated on the same coordinate the fit was made on.

### What did you expect to happen?

The fit should be made on the same coordinate as the evaluation.

### Minimal Complete Verifiable Example

```Python
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt

da = xr.DataArray(np.arange(100), dims=('time',), coords={'time': xr.date_range('2001-01-01', periods=100, freq='YS')})

fit = da.polyfit('time', deg=3)

val = xr.polyval(da.time, fit.polyfit_coefficients)

print(da[0], val[0])
# 0,  31.00174342
# The data is linear, the fit should be exact.

plt.plot(da.time, da, label='Original')
plt.plot(da.time, val, label='Fit')
```


### MVCE confirmation

- [X] Minimal example — the example is as focused as reasonably possible to demonstrate the underlying issue in xarray.
- [X] Complete example — the example is self-contained, including all data and the text of any traceback.
- [X] Verifiable example — the example copy & pastes into an IPython prompt or [Binder notebook](https://mybinder.org/v2/gh/pydata/xarray/main?urlpath=lab/tree/doc/examples/blank_template.ipynb), returning the result.
- [X] New issue — a search of GitHub Issues suggests this is not a duplicate.
- [X] Recent environment — the issue occurs with the latest version of xarray and its dependencies.

### Relevant log output

_No response_

### Anything else we need to know?

`xr.polyval` transforms the temporal coordinate here : https://github.com/pydata/xarray/blob/dbb98b4a40d1679800f7f85d0dad59ef60b5b790/xarray/core/computation.py#L2070

And `da.polyfit` does that here : https://github.com/pydata/xarray/blob/dbb98b4a40d1679800f7f85d0dad59ef60b5b790/xarray/core/dataset.py#L9069-L9077

The results are similar, however, in `polyfit` the new `x` starts at 0, while this offsetting is not done in `polyval`.

This bug was introduced in #9369 I believe. I added a review there woops :sweat:.  Would it be possible to use the same function in both `polyfit` and `polyval` ?

I can make a PR. My intuition would be to use `xr.core.computation._ensure_numeric` in `da.polyfit`.

### Environment

<details>

INSTALLED VERSIONS
------------------
commit: None
python: 3.12.7 | packaged by conda-forge | (main, Oct  4 2024, 16:05:46) [GCC 13.3.0]
python-bits: 64
OS: Linux
OS-release: 6.11.5-200.fc40.x86_64
machine: x86_64
processor: 
byteorder: little
LC_ALL: None
LANG: fr_CA.UTF-8
LOCALE: ('fr_CA', 'UTF-8')
libhdf5: 1.14.3
libnetcdf: 4.9.2

xarray: 2024.10.0
pandas: 2.2.3
numpy: 2.0.2
scipy: 1.14.1
netCDF4: 1.7.1
pydap: None
h5netcdf: 1.4.0
h5py: 3.12.1
zarr: None
cftime: 1.6.4
nc_time_axis: 1.4.1
iris: None
bottleneck: 1.4.2
dask: 2024.10.0
distributed: 2024.10.0
matplotlib: 3.9.2
cartopy: None
seaborn: None
numbagg: None
fsspec: 2024.10.0
cupy: None
pint: 0.24.3
sparse: 0.16.0a10.dev3+ga73b20d
flox: 0.9.12
numpy_groupies: 0.11.2
setuptools: 75.1.0
pip: 24.2
conda: None
pytest: 8.3.3
mypy: 1.13.0
IPython: 8.29.0
sphinx: 8.1.3


</details>


Discussion:
aulemahal: I guess I should also notice @Karl-Krauth. This changes code affected by your PR on how coordinates are transformed to float, but if I understand all of this correctly, it kinda "reverts" to the old result for datetime objects : 0 on 1970-01-01.
aulemahal: @headtr1ck  I would argue that this is more a regression fix than a breaking change ?  The behaviour for datetime coordinates changed in 2024.10.0 and broke the symmetry between `polyfit` and `polyval`. This makes it work again.

Patch:
diff --git a/doc/whats-new.rst b/doc/whats-new.rst
index 52f75c1a..659063f2 100644
--- a/doc/whats-new.rst
+++ b/doc/whats-new.rst
@@ -37,6 +37,8 @@ Bug fixes
 
 - Fix inadvertent deep-copying of child data in DataTree.
   By `Stephan Hoyer <https://github.com/shoyer>`_.
+- Fix regression in the interoperability of :py:meth:`DataArray.polyfit` and :py:meth:`xr.polyval` for date-time coordinates. (:pull:`9691`).
+  By `Pascal Bourgault <https://github.com/aulemahal>`_.
 
 Documentation
 ~~~~~~~~~~~~~
diff --git a/xarray/core/dataset.py b/xarray/core/dataset.py
index b483cfea..bc9360a8 100644
--- a/xarray/core/dataset.py
+++ b/xarray/core/dataset.py
@@ -60,7 +60,7 @@ from xarray.core.common import (
     _contains_datetime_like_objects,
     get_chunksizes,
 )
-from xarray.core.computation import unify_chunks
+from xarray.core.computation import _ensure_numeric, unify_chunks
 from xarray.core.coordinates import (
     Coordinates,
     DatasetCoordinates,
@@ -87,7 +87,6 @@ from xarray.core.merge import (
     merge_coordinates_without_align,
     merge_core,
 )
-from xarray.core.missing import _floatize_x
 from xarray.core.options import OPTIONS, _get_keep_attrs
 from xarray.core.types import (
     Bins,
@@ -9066,22 +9065,14 @@ class Dataset(
         variables = {}
         skipna_da = skipna
 
-        x: Any = self.coords[dim].variable
-        x = _floatize_x((x,), (x,))[0][0]
-
-        try:
-            x = x.values.astype(np.float64)
-        except TypeError as e:
-            raise TypeError(
-                f"Dim {dim!r} must be castable to float64, got {type(x).__name__}."
-            ) from e
+        x = np.asarray(_ensure_numeric(self.coords[dim]).astype(np.float64))
 
         xname = f"{self[dim].name}_"
         order = int(deg) + 1
         lhs = np.vander(x, order)
 
         if rcond is None:
-            rcond = x.shape[0] * np.finfo(x.dtype).eps
+            rcond = x.shape[0] * np.finfo(x.dtype).eps  # type: ignore[assignment]
 
         # Weights:
         if w is not None:


Test patch:
diff --git a/xarray/tests/test_dataset.py b/xarray/tests/test_dataset.py
index c89cfa85..95285b1c 100644
--- a/xarray/tests/test_dataset.py
+++ b/xarray/tests/test_dataset.py
@@ -6720,6 +6720,34 @@ class TestDataset:
             ds.var1.polyfit("dim2", 10, full=True)
             assert len(ws) == 1
 
+    def test_polyfit_polyval(self) -> None:
+        da = xr.DataArray(
+            np.arange(1, 10).astype(np.float64), dims=["x"], coords=dict(x=np.arange(9))
+        )
+
+        out = da.polyfit("x", 3, full=False)
+        da_fitval = xr.polyval(da.x, out.polyfit_coefficients)
+        # polyval introduces very small errors (1e-16 here)
+        xr.testing.assert_allclose(da_fitval, da)
+
+        da = da.assign_coords(x=xr.date_range("2001-01-01", periods=9, freq="YS"))
+        out = da.polyfit("x", 3, full=False)
+        da_fitval = xr.polyval(da.x, out.polyfit_coefficients)
+        xr.testing.assert_allclose(da_fitval, da, rtol=1e-3)
+
+    @requires_cftime
+    def test_polyfit_polyval_cftime(self) -> None:
+        da = xr.DataArray(
+            np.arange(1, 10).astype(np.float64),
+            dims=["x"],
+            coords=dict(
+                x=xr.date_range("2001-01-01", periods=9, freq="YS", calendar="noleap")
+            ),
+        )
+        out = da.polyfit("x", 3, full=False)
+        da_fitval = xr.polyval(da.x, out.polyfit_coefficients)
+        np.testing.assert_allclose(da_fitval, da)
+
     @staticmethod
     def _test_data_var_interior(
         original_data_var, padded_data_var, padded_dim_name, expected_pad_values
@@ -7230,7 +7258,7 @@ def test_differentiate_datetime(dask) -> None:
     assert np.allclose(actual, 1.0)
 
 
-@pytest.mark.skipif(not has_cftime, reason="Test requires cftime.")
+@requires_cftime
 @pytest.mark.parametrize("dask", [True, False])
 def test_differentiate_cftime(dask) -> None:
     rs = np.random.RandomState(42)