Xarray: Convolution operation

[clausmichele]

I would like to perform 2d convolution operation with dask-data stored in xarray, defining the kernel.

Related issue are:
https://github.com/pydata/xarray/issues/1142
https://github.com/pydata/xarray/issues/2010
https://github.com/pydata/xarray/issues/1323

Referring to Dask, they have a specific image-oriented API which allows for convolutions:
http://image.dask.org/en/latest/dask_image.ndfilters.html

And this guy already implemented convolution operation based on xarray and dask:
https://github.com/serazing
https://github.com/serazing/xscale
https://xscale.readthedocs.io/en/latest/api.html

Are there plans to implement this feature?
In the meanwhile I'll try dask-image and xscale.

[fujiisoup]

One thing I would like to implement in somday is multi-dimensional rolling operation.
The 1-dimensional convolution can be done with rolling -> construct -> dot, as can be seen in the doc page
(see the last paragraph of http://xarray.pydata.org/en/stable/computation.html#rolling-window-operations)

This is can be extended to multiple dimensions, but it may not be straightforward.

[clausmichele]

Could you please give an example on how to perform a 2d convolution with a 3x3 kernel?

[fujiisoup]

Maybe we can have a simpler API for convolution operation, though.

In [1]: import numpy as np 
   ...: import xarray as xr 
   ...:  
   ...: da = xr.DataArray(np.random.randn(15, 30), dims=['x', 'y']) 
   ...: kernel = xr.DataArray(np.random.randn(3, 3), dims=['kx', 'ky']) 
   ...:  
   ...: da.rolling(x=3, y=3).construct(x='kx', y='ky').dot(kernel)              
Out[1]: 
<xarray.DataArray (x: 15, y: 30)>
array([[            nan,             nan,             nan,
                    nan,             nan,             nan,
                    nan,             nan,             nan,
                    nan,             nan,             nan,
                    nan,             nan,             nan,
                    nan,             nan,             nan,
                    nan,             nan,             nan,
                    nan,             nan,             nan,
                    nan,             nan,             nan,
                    nan,             nan,             nan],
       [            nan,             nan,             nan,
                    nan,             nan,             nan,
                    nan,             nan,             nan,
                    nan,             nan,             nan,
                    nan,             nan,             nan,
                    nan,             nan,             nan,
                    nan,             nan,             nan,
                    nan,             nan,             nan,
                    nan,             nan,             nan,
                    nan,             nan,             nan],
...
       [            nan,             nan, -2.30319699e-01,
         3.98542408e-01,  7.65734275e+00, -3.78602564e-01,
        -3.79670552e+00, -4.63870114e+00,  3.34264622e-02,
        -3.12097772e+00, -5.76697267e+00,  1.19804861e+00,
        -8.94696248e-01,  2.29308845e+00, -6.39524525e-01,
         4.63574750e+00,  9.72065650e-01, -2.79080617e-01,
        -4.08284408e-01,  4.09498738e+00,  2.21513156e+00,
         2.46188185e-01, -1.30140822e+00, -4.70525588e+00,
        -4.60012056e+00,  2.33333189e-01, -2.86204413e-01,
        -5.63190762e-01,  9.31915537e-01,  7.84188609e-01],
       [            nan,             nan,  1.04286238e+00,
        -1.51693719e+00,  2.49199283e+00,  1.74931359e-01,
        -4.26361392e+00, -1.85066273e-01, -2.45780660e+00,
        -3.20920459e+00, -4.13765502e+00, -3.64119127e+00,
         1.13819179e-01, -2.10588083e-01, -2.58307399e-02,
        -6.73602885e-01,  1.51186293e+00,  2.22395020e+00,
         3.59169613e+00,  4.44203028e+00,  3.15528384e-01,
        -2.30913656e+00,  3.07864240e+00, -9.21743416e-01,
        -2.87995499e+00, -1.92025700e+00, -3.95047208e-01,
         4.60378793e+00,  1.11828099e+00,  4.29419626e-01]])
Dimensions without coordinates: x, y

[fujiisoup]

Maybe we can keep this issue open.

da.convolve(kernel, x='kx', y='ky', mode='same')

would be a possible API?

The contribution will be very much appreciated ;)

[fujiisoup]

Or maybe we can convolve over the shared dimensions.

da = xr.DataArray(np.random.randn(15, 30), dims=['x', 'y']) 
kernel = xr.DataArray(np.random.randn(3, 3), dims=['x', 'y']) 
da.convolve(kernel, mode='same')

Other dimensions maybe broadcasted.