Calculating slope(using polyfit) between sliding window of two xarray

[mohseniaref]

I am trying two calculate slope between sliding window(3x4) of two xarray dataset.you can imagine we want to calculate slope of dx and dy and put the value in the center of sliding window for x=3,y=4. I really appreciate if someone could help me to calculate slope or any function between sliding window of xarray in x,y dimension. Finally ,I want to interpolate the slope results into similar grid dataset of original dataset like dv

https://stackoverflow.com/questions/73596059/calculating-slope-between-sliding-window-of-two-xarray

dv = xr.tutorial.open_dataset('air_temperature_gradient')
dx=dv.Tair.mean(dim='time')
de = dv.dTdx
def linear_trend(x, y):
    pf = np.polyfit(x, y, 1)
    return xr.DataArray(pf[0])
slopes = xr.apply_ufunc(linear_trend,
                        de.rolling(lat=3, lon=4, center=True).construct('window_dim'), dx.rolling(lat=3, lon=4, center=True).construct('window_dim'),
                        vectorize=True,
                        )

[mathause]

Can you use de.rolling(...).construct(...) and then use the resulting array to calculate the trend. Not very memory-efficient but I don't see a way to use the rolling object directly.

[mathause]

The error message is not very helpful - you need to provide the names of the new dimensions:

de.rolling(lat=3, lon=4, center=True).construct(lat="lat_window", lon="lon_window")

[mohseniaref]

Good Point! Still couldn't solve it .The dataset is part of tutroial in xarray so you can also try it.

for instance

de.rolling(lat=3, lon=4, center=True).construct(
    lat="lat_window", lon="lon_window")

you will have

xarray.DataArray'dTdx'time: 2920lat: 25lon: 53lat_window: 3lon_window: 4

I have tried following

import xarray as xr 
import numpy as np 
dv = xr.tutorial.open_dataset('air_temperature_gradient')
dx = dv.Tair.mean(dim='time')
de = dv.dTdx
slopes = de.rolling(lat=3, lon=4, center=True).construct(lat="lat_window", lon="lon_window").polyfit(dx.rolling(lat=3, lon=4, center=True).construct(lat="lat_window", lon="lon_window"),deg=1)

I got following error message after running the command

TypeError                                 Traceback (most recent call last)
/tmp/ipykernel_14835/846679853.py in <module>
 de = dv.dTdx

 slopes = de.rolling(lat=3, lon=4, center=True).construct(lat="lat_window", lon="lon_window").polyfit(dx.rolling(lat=3, lon=4, center=True).construct(lat="lat_window", lon="lon_window"),deg=1)

~/anaconda3/envs/code/lib/python3.8/site-packages/xarray/core/dataarray.py in polyfit(self, dim, deg, skipna, rcond, w, full, cov)
   4401         xarray.polyval
   4402         """
-> 4403         return self._to_temp_dataset().polyfit(
   4404             dim, deg, skipna=skipna, rcond=rcond, w=w, full=full, cov=cov
   4405         )

~/anaconda3/envs/code/lib/python3.8/site-packages/xarray/core/dataset.py in polyfit(self, dim, deg, skipna, rcond, w, full, cov)
   7477         skipna_da = skipna
   7478 
-> 7479         x = get_clean_interp_index(self, dim, strict=False)
   7480         xname = f"{self[dim].name}_"
   7481         order = int(deg) + 1

~/anaconda3/envs/code/lib/python3.8/site-packages/xarray/core/missing.py in get_clean_interp_index(arr, dim, use_coordinate, strict)
    259 
    260     if use_coordinate is True:
--> 261         index = arr.get_index(dim)
...
--> 467         return key in self.mapping
    468 
    469     def __repr__(self) -> str:

TypeError: unhashable type: 'DataArray'

[mathause]

polyfit only works along one dimension and the x-values need to be coordinates i.e. you need to pass the name of the dimension you want to apply it along.

Overall, I think I don't understand what you are trying to achieve - can you try to first solve the problem in 1D? Or with pure numpy arrays?

[mohseniaref]

Ah..I see, Yes, for sure the numpy is faster .Even one idea is using from skimage.util import view_as_windows .I was wondering if it can be done first with xarray

My problem is 2D or 3D .I have tried first with following function and then I have tried with polyfit!

from skimage.util import view_as_windows

I modified the code but still didn't work!

def linear_trend(x, y):
    idx = np.isfinite(x) & np.isfinite(y)
    pf = np.polyfit(x[idx], y[idx], 1)
    return xr.DataArray(pf[0])


slopes = xr.apply_ufunc(linear_trend,de.rolling(lat=3, lon=4, center=True).construct(lat="lat_window", lon="lon_window"), dx.rolling(
                            lat=3, lon=4, center=True).construct(lat="lat_window", lon="lon_window"),
                        vectorize=True,
                        )
slopes

let's imagine you have two image from surface temprature and elevation of an area. you can easily calculate regression between two image and extract slope and intercept.I want to do this instead of whole image for smaller window siliding over image and calculate just for center of sliding window slope value. I want to do something like this with time series of surface deformation and elevation. Does it help ?

[mathause]

You may want to go through the apply_ufunc tutorial and check the trick on outputting what apply_ufunc feeds your function, i.e.:

def linear_trend(x, y):
    print(f"{x.shape=}")
    print(f"{y.shape=}")
    ...

You may have to play with the input_core_dims.