issue_comments

@clyne Let me rephrase my question: how do you feel about xarray providing a polyfit/polyval implementation essentially duplicating GeoCat's implementation ?

@maboualidev Is your objective to integrate the GeoCat implementation into xarray or keep it standalone ?

On my end, I'll submit a PR to add support for non-standard calendars to xarray.core.missing.get_clean_interp, which you'd then be able to use to get x values from coordinates.

@maboualidev Nice ! I see you're storing the residuals in the DataArray attributes. From my perspective, it would be useful to have those directly as DataArrays. Thoughts ?

So it looks like there are multiple inspirations to draw from. Here is what I could gather.

• xscale.signal.fitting.polyfit(obj, deg=1, dim=None, coord=None) supports chunking along the fitting dimension using dask.array.linalg.lstsq. No explicit missing data handling.
• xyzpy.signal.xr_polyfit(obj, dim, ix=None, deg=0.5, poly='hermite') applies np.polynomial.polynomial.polyfit using xr.apply_ufunc along dim with the help of numba. Also supports other types of polynomial (legendre, chebyshev, ...). Missing values are masked out 1D wise.
• geocat.comp.ndpolyfit(x: Iterable, y: Iterable, deg: int, axis: int = 0, **kwargs) -> (xr.DataArray, da.Array) reorders the array to apply np.polyfit along dim, returns the full outputs (residuals, rank, etc) as DataArray attributes. Missing values are masked out in bulk if possible, 1D-wise otherwise.

There does not seem to be matching polyval implementations for any of those nor support for indexing along a time dimension with a non-standard calendar.

My current implementation is pretty naive. It's just calling numpy.polyfit using dask.array.apply_along_axis. Happy to put that in a PR as a starting point, but there are a couple of questions I had: * How to return the full output (residuals, rank, singular_values, rcond) ? A tuple of dataarrays or a dataset ? * Do we want to use the dask least square functionality to allow for chunking within the x dimension ? Then it's not just a simple wrapper around polyfit. * Should we use np.polyfit or np.polynomial.polynomial.polyfit ?

Thanks, it seems to work !

Started to work on this and facing some issues with the x-coordinate when its a datetime. For standard calendars, I can use pd.to_numeric(da.time), but for non-standard calendars, it's not clear how to go ahead. If I use xr.coding.times.encode_cf_datetime(coord), the coefficients we'll find will only make sense in the polyval function if we use the same time encoding.