search

xarray-contrib / pint-xarray

Interface for using pint with xarray, providing convenience accessors
https://pint-xarray.readthedocs.io/en/latest/
Apache License 2.0
101 stars 12 forks source link

Extracting annual amplitude and phase of xarray dataset for each pixel using xarray.DataArray.curvefit #183

Closed mohseniaref closed 2 years ago

mohseniaref commented 2 years ago

I was wondering how can I extract annual amplitude and phase of xarray time series using xarray.DataArray.curvefit. We fit a 1d function in time to return annual and seasonal amplitude and phase with dims (x, y)

We can use formula similar to https://stats.stackexchange.com/questions/77543/how-do-i-get-the-amplitude-and-phase-for-sine-wave-from-lm-summary but I have difficulty to extract a0,a1,a2 pixelwise for the dataset.I have tried also to convert time to julian day to use instead of ds.time in ds.curvefit but it didn't work. I really appreciate if you can help me. Similar problem with numpy is solved in https://stackoverflow.com/questions/15094619/fitting-a-3d-array-of-data-to-a-1d-function-with-numpy-or-scipy

ds = xr.tutorial.open_dataset('air_temperature')
ds['air2'] = ds.air.copy()
timejulian=ds.time.dt.strftime('%y%j')
def timeseries_function_season(x, a0, a1, a2):
    return a0+(a1*np.cos((2*np.pi/365)*x)+a2*np.sin((2*np.pi/365)*x))
dn = ds.curvefit('time', func=timeseries_function_season)
keewis commented 2 years ago

err... how does that relate to pint-xarray?

As far as I can tell, you already asked exactly this question in pydata/xarray#6968, and I think we should keep that discussion there.

A different issue would be whether curve_fit propagates units. I would guess it doesn't since it uses scipy.optimize underneath, so I'd probably put it in the same category in Working with numpy-like arrays / Missing features as interp.