Plans for a N-dimensional convolution operator for xarray DataArray ?

[julien....@gmail.com]

Hi all,

I am starting to think about the best approach for implementing spatial filtering methods to xarray DataArray with dask.

I was wondering whether one of you had plans for implementing a N-dimensional convolution operator for xarray.
I was thinking of something close to scipy.signal.convolve or one of its equivalent.

I guess a straightforward approach is to use a numpy implementation and map_block/map_overlap methods, right ? 

thanks for your comments/feedbacks.

cheers.

--
J.

[Ryan Abernathey]

Julien,

The big question is whether you want to convolve __across_ dask blocks or __within__ dask blocks.

If you are convolving within the blocks (e.g. you want to do smoothing in space and your dataset is chunked in time), it is easy: you use dask.map_blocks on the underlying xarray.Dataset.data objects. I have an example notebook for this:

https://gist.github.com/rabernat/3e2c3645d4666674cf2b

I use the astropy convolution operators because they deal with missing data better. (Note: this application motivated a pull request to astropy to release the GIL for better parallel performance. This was merged in version 1.0.4.)

The one downside of this approach is that you have to work with the low-level dask objects. There has been talk about providing a higher-level dask wrapper for this functionality, but it has not been implemented: https://github.com/pydata/xarray/issues/585

If you want to convolve across dask blocks (e.g. data is chunked in all dimensions), that would have to be implemented in dask. It doesn't sound easy to me, but maybe Matt or Stephen could say otherwise.

-Ryan

--
You received this message because you are subscribed to the Google Groups "xarray" group.
To unsubscribe from this group and stop receiving emails from it, send an email to xarray+un...@googlegroups.com.
To post to this group, send email to xar...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/xarray/8a9031e4-a9d1-4cee-ad61-c722606c3f74%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[julien....@gmail.com]

thanks Ryan for your reply and the references.

I have actually made some tests for convolution across dask blocks

(Your notebook  http://nbviewer.jupyter.org/gist/rabernat/3e2c3645d4666674cf2b was actually very useful)

The preliminary results are available here :
https://github.com/lesommer/notebooks/blob/master/OptimizingPyClimLanczosFilter.ipynb

it seems that using map_block with a depth specified from the width of the averaging window works ok.

But this is done at the dask level and I was wondering how to do the same with xarray dask.

[Matthew Rocklin]

Right, relevant dask.array documentation for mapping with a slight block overlap is here:

• API: http://dask.pydata.org/en/latest/array-api.html#dask.array.Array.map_overlap
  
• Internal: http://dask.pydata.org/en/latest/array-ghost.html

To view this discussion on the web visit https://groups.google.com/d/msgid/xarray/b2ede4a4-e679-4377-ac9b-6e05dfd83fd1%40googlegroups.com.

[Stephan Hoyer]

Wrapping dask array functions for xarray is usually pretty straightforward, e.g.,

data = orig_data_array.chunk().data

new_data = da.map_blocks(data, ...)

new_data_array = xr.DataArray(new_data, orig_data_array.coords)

If necessary, you can get axis numbers by calling the get_axis_num method.

To view this discussion on the web visit https://groups.google.com/d/msgid/xarray/CAJ8oX-Gp9zrg7Ck3SDy9WWHe3D1jYn_xkbr_cer8ozHLULB6qg%40mail.gmail.com.

[julien....@gmail.com]

Thank you all for your answers. I will proceed with this approach.
--
J.