primap-community / primap2

The next generation of the PRIMAP climate policy analysis suite
https://primap2.readthedocs.io
Apache License 2.0
9 stars 2 forks source link

xarray polyfit time coordinate treatment #293

Open JGuetschow opened 1 week ago

JGuetschow commented 1 week ago

When trying to solve #292 I found an unexpected behaviour in xarray's polyfit routine. This issue is to discuss how we deal with that

In issue #292 I found the reason for the problem with using xr.polyfit and xr.polyval. It's not (only) the time resolution as I first suspected but a coordinate recomputation in xr.polyfit. The function _floatize_x (https://github.com/pydata/xarray/blob/91962d6aec380cb83fe80b2afdfa556efdd817a3/xarray/core/missing.py#L585) shifts the time coordinate such that the minimal value is 0 to avoid accuracy problems as np.datetime64 uses 64 bit integers which can not be represented by 64 bit floats with the necessary accuracy. All further computations are done with the shifted coordinates and thus the coefficients are also computed in relation to the shifted coordinates. 6As 1970 is zero this case worked while other dates did not work. The following code takes the shifting of _floatize_x into account:

def test_temp_polyval():
    test_ts = xr.DataArray(
        np.linspace(6, 12, 11),
        coords={"time": pd.date_range("1956-01-01", "1966-01-01", freq="YS")},
        dims="time",
        name="test_ts",
    )

    time_to_eval = np.datetime64('1957-01-01')

    fit = test_ts.polyfit(dim='time', deg=1, skipna=True)
    value = xr.polyval(
        test_ts.coords['time'].loc[{'time': [time_to_eval]}]-test_ts.coords['time'].data[0],
        fit.polyfit_coefficients
    )

    value.name='test_ts' # for assertion

    assert_aligned_equal(
        test_ts.loc[{'time': [time_to_eval]}],
        value,
        rtol=1e-03,
    )
    return None

I still think this is very weird behaviour. Maybe I missed a hint when reading the docs.

Originally posted by @JGuetschow in https://github.com/primap-community/primap2/issues/292#issuecomment-2468361809

I personally think this is bug in xarray, but before I open an issue, maybe you can add your view @mikapfl?

JGuetschow commented 6 days ago

This behaviour is only in xarray 2024.10.0 and not in 2024.9.0. I have excluded 2024.10.0 for now.

JGuetschow commented 6 days ago

Already fixed in xarray just not released yet. As we have to change the allowed xarray versions I'll keep this open until the bug has been released.