Feature: N-dimensional auto_combine

[TomNicholas]

What I did

Generalised the auto_combine() function to be able to concatenate and merge datasets along any number of dimensions, instead of just one.

Provides one solution to #2159, and relevant for discussion in #2039.

Currently it cannot deduce the order in which datasets should be concatenated along any one dimension from the coordinates, so it just concatenates them in order they are supplied. This means for an N-D
concatenation the datasets have to be supplied as a list of lists, which is nested as many times as there are dimensions to be concatenated along.

How it works

In _infer_concat_order_from_nested_list() the nested list of datasets is recursively traversed in order to create a dictionary of datasets, where the keys are the corresponding "tile IDs". These tile IDs are
tuples serving as multidimensional indexes for the position of the dataset within the hypercube of all datasets which are to be combined. For example four datasets which are to be combined along two dimensions would be supplied as

datasets = [[ds0, ds1], [ds2, ds3]]

and given tile_IDs to be stored as

combined_ids = {(0, 0): ds0, (0, 1): ds1, (1, 0): ds2, (1, 1): ds3}

Using this unambiguous intermediate structure means that another method could be used to organise the datasets for concatenation (i.e. reading the values of their coordinates), and a new keyword argument
infer_order_from_coords used to choose the method. The _combine_nd() function concatenates along one dimension at a time, reducing the length of the tile_ID tuple by one each time _combine_along_first_dim() is called. After each concatenation the different variables are merged, so the new auto_combine() is essentially like calling the old one once for each dimension in concat_dims.

Still to do

I would like people's opinions on the method I've chosen to do this, and any feedback on the code quality would be appreciated. Assuming we're happy with the method I used here, then the remaining tasks include:

• More tests of the final auto_combine() function

• Add option to deduce concatenation order from coords (or this could be a separate PR)

• Integrate this all the way up to open_mfdataset().

• Unit tests for open_mfdataset()

• More tests that the user has inputted a valid structure of datasets

• Possibly parallelize the concatenation step?

• A few other small TODOs which are in combine.py

• Proper documentation showing how the input should be structured.

• Fix failing unit tests on python 2.7 (though support for 2.7 is being dropped at the end of 2018?)

• Fix failing unit tests on python 3.5

• Update what's new

This PR was intended to solve the common use case of collecting output from a simulation which was parallelized in multiple dimensions. I would like to write a tutorial about how to use xarray to do this, including examples of how to preprocess the data and discard processor ghost cells.

[pep8speaks]

Hello @TomNicholas! Thanks for updating the PR.

• In the file xarray/backends/api.py, following are the PEP8 issues :

Line 647:80: E501 line too long (81 > 79 characters)

• In the file xarray/core/combine.py, following are the PEP8 issues :

Line 381:36: E226 missing whitespace around arithmetic operator
Line 536:75: E211 whitespace before '('
Line 537:57: E126 continuation line over-indented for hanging indent

• In the file xarray/testing.py, following are the PEP8 issues :

Line 148:1: W391 blank line at end of file

• There are no PEP8 issues in the file xarray/tests/__init__.py !

• In the file xarray/tests/test_backends.py, following are the PEP8 issues :

Line 2165:36: E127 continuation line over-indented for visual indent
Line 2172:36: E126 continuation line over-indented for hanging indent

• In the file xarray/tests/test_combine.py, following are the PEP8 issues :

Line 495:65: E211 whitespace before '('
Line 523:52: E231 missing whitespace after ','
Line 523:54: E231 missing whitespace after ','
Line 523:61: E231 missing whitespace after ','
Line 524:57: E241 multiple spaces after ','
Line 525:55: E241 multiple spaces after ','
Line 526:55: E241 multiple spaces after ','
Line 527:57: E241 multiple spaces after ','
Line 597:52: E203 whitespace before ','
Line 646:80: E501 line too long (85 > 79 characters)
Line 647:80: E501 line too long (85 > 79 characters)
Line 658:80: E501 line too long (83 > 79 characters)
Line 659:80: E501 line too long (87 > 79 characters)

Comment last updated on December 12, 2018 at 00:54 Hours UTC

[shoyer]

a few high level comments

[shoyer]

I should have suggested this before, but the internal implementation of numpy.block() is a nice example of this sort of logic: https://github.com/numpy/numpy/blob/e726c045c72833d5c826e8a13f889746ee0bfdf4/numpy/core/shape_base.py#L661

[TomNicholas]

This is basically done now - I've implemented everything I wanted to, and included unit tests for the new functionality.
I'm also successfully using it in my personal code now.

@shoyer I haven't changed the way the combine function is applied repeatedly to match the implementation in numpy.block(), I think the way I've done it is relatively clear, but let me know if you think it should be refactored.

I don't understand why some of the CI tests are failing - all test which run on my machine pass, and the errors in the log seem to come from dask arrays not being loaded, i.e:

AssertionError: 
    <xarray.Dataset>
    Dimensions:  (x: 10, y: 8)
    Dimensions without coordinates: x, y
    Data variables:
        foo      (x, y) float64 -1.19 0.4476 0.6283 ... -0.09699 -0.2311 -0.6482
    <xarray.Dataset>
    Dimensions:  (x: 10, y: 8)
    Dimensions without coordinates: x, y
    Data variables:
        foo      (x, y) float64 dask.array<shape=(10, 8), chunksize=(5, 4)>

[TomNicholas]

Also if you specify dim='something' when that dimension doesn't already exist, you are explicitly telling it to add a new dimension.

[shoyer]

Specifying dim=None will specify (1) explicitly, but otherwise yes, that's exactly what it does.

Yes, this was a typo on my part :).

[TomNicholas]

I appreciate that this is pretty complicated, perhaps it should have its own section in the docs (I can do another PR?)

[shoyer]

Do you think it would make sense to try to get rid of "both merging and concatenating" in favor of requiring another dimension in the grid? I guess we do sort of need this for this for open_mfdataset but it feels more complex than desired.

[TomNicholas]

I think that would mean there might be some situations that 1D auto_combine() could deal with but nested_concat_and_merge() could not. Situations where you can only merge the variables once they have been individually concatenated. (actually I'm not sure - I would need to think about that)

I think this is a question of whether you think that:

a) auto_combine() is a well-defined operation that is still well-defined in N-D

or

b) auto_combine()'s purpose is to do it's best to combine whatever is thrown at it, and more surgical N-D operations should be the job of another function

I personally think a), but I expect users who haven't read the source code for auto_combine() might disagree.

[shoyer]

I think that would mean there might be some situations that 1D auto_combine() could deal with but nested_concat_and_merge() could not. Situations where you can only merge the variables once they have been individually concatenated. (actually I'm not sure - I would need to think about that)

I think multiple levels nested_concat_and_merge() could always replace auto_combine() as long as merging happens last, at the outer-most level. That's basically what auto_combine() already does.

There are probably some edge cases where the existing auto_combine would you be a little more sloppy about how you put together your lists of variables. But I don't think that's really a good idea for us to support in xarray, since it removes most of the freedom from the implementation.

I personally think a), but I expect users who haven't read the source code for auto_combine() might disagree.

Yes, this is my concern. The way this API handles nested lists makes sense from an implementation perspective, but not really from a user perspective. For users, I think it makes sense to have:

1. The precisely defined nested_concat_and_merge() for cases where you already data sorted. This is pretty common (e.g., with hierarchical directory structures).
2. A magical auto_combine() that inspects coordinates to figure out how to put everything together. We might even rename this to combine() (if we think it's comprehensive enough).

The original version of auto_combine() that I wrote was aiming towards this second use-case, but was obviously fragile and incomplete.

If you like, we could also merge this work (which is excellent progress towards user-facing APIs) but keep the changes internal to xarray for now until we figure out the public APIs.

[TomNicholas]

There are probably some edge cases where the existing auto_combine would you be a little more slopping about how you put together your lists of variables.

I think that's probably the case, but I also think that those edge cases will be so specific that maybe we don't have to explicitly support them. We could just say that anyone who has a combination of datasets that is that funky can just concatenate them themselves?

For users, I think it makes sense to have...

I agree, two separate functions is a lot more intuitive than having auto_combine() behave either auto-magically or totally manually depending on the value of an optional keyword argument.

The only problem with that idea is that both of these functions should be options for open_mfdataset(), so you still have multiple possible argument structures for that function. I think that's okay though - it would be a lot easier to explain if you say "if the keyword argument combine='auto' then the function auto_combine() is used, if combine='manual' then nested_concat_and_merge() is used, and the filepaths must be structured accordingly."

If you like, we could also merge this work (which is excellent progress towards user-facing APIs) but keep the changes internal to xarray for now until we figure out the public APIs.

That would be great. Then I could start using the master branch of xarray again in my code, while we redo the public API. If I set concat_dims back to concat_dim then it should be completely backwards-compatible.

[shoyer]

The only problem with that idea is that both of these functions should be options for open_mfdataset(), so you still have multiple possible argument structures for that function. I think that's okay though - it would be a lot easier to explain if you say "if the keyword argument combine='auto' then the function auto_combine() is used, if combine='manual' then nested_concat_and_merge() is used, and the filepaths must be structured accordingly."

This sounds pretty good to me.

If you like, we could also merge this work (which is excellent progress towards user-facing APIs) but keep the changes internal to xarray for now until we figure out the public APIs.

That would be great. Then I could start using the master branch of xarray again in my code, while we redo the public API. If I set concat_dims back to concat_dim then it should be completely backwards-compatible.

OK, do you want to go ahead and revert the documentation/public API for now?

I would even be OK supporting nested lists temporarily in xarray via APIs like open_mfdataset and auto_concat without official documentation/support -- using this feature would be "at your own risk". As long as we don't break any existing code, that is totally fine. We just shouldn't advertise an interface until we've happy with it.

[TomNicholas]

This sounds pretty good to me.

Okay good. So with this API then auto_combine() (and therefore the default behaviour of open_mfdatset()) would be to throw a warning if there isn't enough info in the dataset coords to order them?

OK, do you want to go ahead and revert the documentation/public API for now?

Yes, I'll do that.

I would even be OK supporting nested lists temporarily ... "at your own risk"

If I revert the documentation and you merge this PR then that's exactly what we will have, which would be useful for me until we do the public API.

(Also it seems that whatever was wrong with cftime has been fixed now, as the CI tests are passing.)

[shoyer]

Okay good. So with this API then auto_combine() (and therefore the default behaviour of open_mfdatset()) would be to throw a warning if there isn't enough info in the dataset coords to order them?

Yes, but throwing a warning should probably only be a temporary solution. Long term we should pick a default value (e.g., combine='auto') and require explicitly opting in to different behavior (e.g., combine='manual').

[TomNicholas]

Okay I've reverted the API, so it basically converts a concat_dim input to it's multidimensional version before continuing.

If we merge this then should I start a separate Pull Request for further discussion about the API?

One of the Travis CI builds failed but again I don't think that was me.

[shoyer]

you can indeed ignore the test failure on dask-dev -- that looks like a dask issue (dask/dask#4291)

[TomNicholas]

Removed the unnecessary argument.

[shoyer]

thanks @TomNicholas! I'm looking forward to the follow-ups :)