What's the best way to perform a convolution (or rolling mean) on an irregularly-spaced grid?

[tomchor]

I'm trying to perform a convolution on 2D data that's on an irregularly-spaced grid. Reading past issues, I came across #4196, which sounds like I could do what I want, but it looks like it doesn't take into account the irregular spacing. For example this MWE:

In [29]: import numpy as np
    ...: import xarray as xr
    ...: import random
    ...: np.random.seed(10)
    ...: 
    ...: da = xr.DataArray(np.random.randn(5, 10), dims=['x', 'y'], coords=dict(x=[0, 1, 3, 6, 10],
    ...:                                                                        y=range(10)))
    ...: kernel = xr.DataArray(np.random.randn(3, 3), dims=['kx', 'ky'], coords=dict(kx=[-1, 0, 1],
    ...:                                                                             ky=[-1, 0, 1]))
    ...: 
    ...: da.rolling(x=3, y=3).construct(x='kx', y='ky').dot(kernel)
    ...: 
Out[29]: 
<xarray.DataArray (x: 5, y: 10)>
array([[        nan,         nan,         nan,         nan,         nan,
                nan,         nan,         nan,         nan,         nan],
       [        nan,         nan,         nan,         nan,         nan,
                nan,         nan,         nan,         nan,         nan],
       [        nan,         nan, -4.23539009, -4.41080965,  2.63043495,
         8.13440411,  4.15214765,  3.36171765,  1.18381928,  4.24895863],
       [        nan,         nan, -3.1161138 , -2.53048603,  2.35874729,
         2.8133469 ,  6.22844301,  2.01160326,  0.96701393,  2.32143665],
       [        nan,         nan,  0.14328192, -3.18686044, -2.57300755,
        -1.72151673, -0.51058699, -1.21852662,  0.21368092,  1.50704157]])
Coordinates:
  * x        (x) int64 0 1 3 6 10
  * y        (y) int64 0 1 2 3 4 5 6 7 8 9

gives me the result regardless of whether or not the directions are irregularly-spaced.

The kernel itself is not important, so I'd be happy with a rolling mean that takes the grid spacing into account (which should be equivalent to a top-hat kernel).

Is there a solution already implemented for this?

[TomNicholas]

If you're asking if there is a DataArray.rolling.weighted.mean() then there currently is not one implemented.

Can you instead do this by passing a custom reduction function to Dataset.rolling.reduce? Your custom function would need to look at the values of the (irregularly spaced) coordinate variables as well as at the data, and compute your kernel/mean correctly using that information.

[tomchor]

Thanks for the answer! I actually wasn't aware of Dataset.rolling.reduce, but that seems like it could solve the problem. I can't fully understand how to use it though. For example, let's say I try to construct a Gaussian kernel using the rolling operation below:

In [25]: import numpy as np
    ...: import xarray as xr
    ...: import random
    ...: np.random.seed(10)
    ...: 
    ...: da = xr.DataArray(np.random.randn(5, 10), dims=['x', 'y'], coords=dict(x=[0, 1, 3, 6, 10],
    ...:                                                                        y=range(10)))
    ...: 
    ...: rolling = da.rolling(x=3, y=3, center=True).construct(x="kx", y="ky")
    ...: 

I'm not sure how to access the values of x and y at each point x and y, since:

In [26]: rolling.isel(x=1, y=1)
Out[26]: 
<xarray.DataArray (kx: 3, ky: 3)>
array([[ 1.3315865 ,  0.71527897, -1.54540029],
       [ 0.43302619,  1.20303737, -0.96506567],
       [-1.97772828, -1.7433723 ,  0.26607016]])
Coordinates:
    x        int64 1
    y        int64 1
Dimensions without coordinates: kx, ky

Sorry if this isn't very clear. I'm having a hard time finding the right way to formulate my question!

In any case, I'll keep trying and maybe come up with a better way to voice my doubts :)

[dcherian]

I think your first approach of using .dot is right. The kernel needs to take care of the irregular spacing so it needs to be 3D: kx, ky, x. You could construct the appropriate weights by working with da.x.rolling(x=3).construct.

You could also iterate over the rolling object for a much slower, but simpler to implement method. Probably good for debugging.