concat changes variable order

[kmuehlbauer]

Code Sample, a copy-pastable example if possible

A "Minimal, Complete and Verifiable Example" will make it much easier for maintainers to help you:
http://matthewrocklin.com/blog/work/2018/02/28/minimal-bug-reports

• Case 1: Creation of Dataset without Coordinates

data = np.zeros((2,3))
ds = xr.Dataset({'test': (['c', 'b'],  data)})
print(ds.dims)
ds2 = xr.concat([ds, ds], dim='c')
print(ds2.dims)

yields (assumed correct) output of:

Frozen(SortedKeysDict({'c': 2, 'b': 3}))
Frozen(SortedKeysDict({'c': 4, 'b': 3}))

• Case 2: Creation of Dataset with Coordinates

data = np.zeros((2,3))
ds = xr.Dataset({'test': (['c', 'b'],  data)}, 
                coords={'c': (['c'], np.arange(data.shape[0])),
                        'b': (['b'], np.arange(data.shape[1])),})
print(ds.dims)
ds2 = xr.concat([ds, ds], dim='c')
print(ds2.dims)

yields (assumed false) output of:

Frozen(SortedKeysDict({'c': 2, 'b': 3}))
Frozen(SortedKeysDict({'b': 3, 'c': 4}))

Problem description

xr.concat changes the dimension order for .dims as well as .sizes to an alphanumerically sorted representation.

Expected Output

xr.concat should not change the dimension order in any case.

Frozen(SortedKeysDict({'c': 2, 'b': 3}))
Frozen(SortedKeysDict({'c': 4, 'b': 3}))

Output of xr.show_versions()

INSTALLED VERSIONS ------------------ commit: None python: 3.7.1 | packaged by conda-forge | (default, Nov 13 2018, 18:33:04) [GCC 7.3.0] python-bits: 64 OS: Linux OS-release: 4.12.14-lp150.12.48-default machine: x86_64 processor: x86_64 byteorder: little LC_ALL: None LANG: de_DE.UTF-8 LOCALE: de_DE.UTF-8 libhdf5: 1.10.4 libnetcdf: 4.6.2

xarray: 0.11.3
pandas: 0.24.1
numpy: 1.16.1
scipy: 1.2.0
netCDF4: 1.4.2
pydap: None
h5netcdf: 0.6.2
h5py: 2.9.0
Nio: None
zarr: None
cftime: 1.0.3.4
PseudonetCDF: None
rasterio: None
cfgrib: None
iris: None
bottleneck: 1.2.1
cyordereddict: None
dask: None
distributed: None
matplotlib: 3.0.2
cartopy: 0.17.0
seaborn: None
setuptools: 40.8.0
pip: 19.0.2
conda: None
pytest: 4.2.0
IPython: 7.2.0
sphinx: None

[kmuehlbauer]

Xref: Gitter Chat

[kmuehlbauer]

This has also implications for the output using .to_netcdf(). If we read a netcdf dataset (same structure as above) with xr.open_dataset and then do the above xr.concat and save the resulting dataset with .to_netcdf then the dimensions of the dataset will be reversed in the resulting file.

Now, as the xr.concat operation need to change the length of the dimension ('c', which is not allowed by netCDF library), this is done by creating a new dataset. In this creation process xarray obviously uses the alphanumerically sorted representation of the source dataset dimension's and not the creation order as in the source dataset.

I did not find any hints in the docs on that topic. I need to preserve the original dimension ordering as declared in the source dataset. How can I achieve this using xarray?

[shoyer]

Your system might print dataset dimensions like Frozen(SortedKeysDict({'c': 2, 'b': 3})), but the iteration order will always be sorted (including if you write the dataset to disk as netcdf file).

When we drop support for Python 3.5, xarray might switch to dimensions matching order of insertion, since we'll get that for free with Python dictionary. But I still doubt we would make any guarantees about preserving dimension order in xarray operations, just like we don't guarantee variable order as part of xarray's API. It should be deterministic (with fixed versions of xarray and dependencies), but you shouldn't write your code in a way that breaks if changes.

What's your actual use-case here? What are you trying to do that needs preserving of dimension order?

[kmuehlbauer]

Thanks for looking into this @shoyer.

Your system might print dataset dimensions like Frozen(SortedKeysDict({'c': 2, 'b': 3})), but the iteration order will always be sorted (including if you write the dataset to disk as netcdf file).

This isn't true for my system. If we consider this example:

data = np.zeros((2,3))
ds = xr.Dataset({'test': (['c', 'b'],  data)}, 
                coords={'c': (['c'], np.arange(data.shape[0])),
                        'b': (['b'], np.arange(data.shape[1])),})
ds.to_netcdf('test_dims.nc')
ds2 = xr.concat([ds, ds], dim='c')
ds2.to_netcdf('test_dims2.nc')

Dumping the created files gives the following:

netcdf test_dims {
dimensions:
        c = 2 ;
        b = 3 ;
variables:
        double test(c, b) ;
                test:_FillValue = NaN ;
        int64 c(c) ;
        int64 b(b) ;
data:

 test =
  0, 0, 0,
  0, 0, 0 ;

 c = 0, 1 ;

 b = 0, 1, 2 ;
}

netcdf test_dims2 {
dimensions:
        b = 3 ;
        c = 4 ;
variables:
        int64 b(b) ;
        double test(c, b) ;
                test:_FillValue = NaN ;
        int64 c(c) ;
data:

 b = 0, 1, 2 ;

 test =
  0, 0, 0,
  0, 0, 0,
  0, 0, 0,
  0, 0, 0 ;

 c = 0, 1, 0, 1 ;
}

My use case is, well, I have to use some legacy code.

Concerning my code, yes I'm trying to write it as robust as possible. Finally I wan't to replace the legacy code with the implementation relying completely on xarray, but that's a long way to go.

[shoyer]

Dimensions are written to netCDF files in the order in which they appear on variables in the Dataset:

xarray/xarray/backends/common.py

Lines 325 to 329 in f382fd8

	dims = OrderedDict()
	for v in unlimited_dims: # put unlimited_dims first
	dims[v] = None
	for v in variables.values():
	dims.update(dict(zip(v.dims, v.shape)))

It sounds like your use-case is writing netCDF files to disk with a desired dimension order? We could conceivably add an "encoding" option to datasets for specifying dimension order, like how we support controlling unlimited dimensions.

[kmuehlbauer]

Dimensions are written to netCDF files in the order in which they appear on variables in the Dataset:

I was assuming something along that lines. But in my variable test the 'c' dim is first. And it is written correctly, if there are no coordinates in that dataset (test_dim0). If there are coordinates (test_dims2) the dimensions are written in wrong order. So there is something working in one config and not in the other.

My use case is, that the dimensions should appear in the same order as in the source files.

[shoyer]

I was assuming something along that lines. But in my variable test the 'c' dim is first. And it is written correctly, if there are no coordinates in that dataset (test_dim0). If there are coordinates (test_dims2) the dimensions are written in wrong order. So there is something working in one config and not in the other.

The order of dimensions in the netCDF file matches the order of their appearance on variables in the netCDF files. In your first file, it's (c, b) on variable test. In the second file, it's b on variable b and c from variable test.

My use case is, that the dimensions should appear in the same order as in the source files.

Sorry, xarray is not going to satisfy this use. If you want this guarantee in all cases, you should pick a different tool.

[kmuehlbauer]

@shoyer I'm sorry if I did not explain well enough and if my intentions were vague. So let me first clarify, I really appreciate all your hard work to make xarray better. I've adapted many of my workflows to use xarray and I'm happy that such a library exist.

Let's consider just one more example where I hopefully get better to the point of my problems in understanding.

Two files are created, same dimensions, same data, but one without coordinates the other with coordinates.

data = np.zeros((2,3))
src_dim0 = xr.Dataset({'test': (['c', 'b'],  data)})
src_dim0.to_netcdf('src_dim0.nc')

src_dim1 = xr.Dataset({'test': (['c', 'b'],  data)}, 
                      coords={'c': (['c'], np.arange(data.shape[0])),
                              'b': (['b'], np.arange(data.shape[1])),})
src_dim1.to_netcdf('src_dim1.nc')

The dump of both:

netcdf src_dim0 {
dimensions:
	c = 2 ;
	b = 3 ;
variables:
	double test(c, b) ;
		test:_FillValue = NaN ;
data:

 test =
  0, 0, 0,
  0, 0, 0 ;
}

netcdf src_dim1 {
dimensions:
	c = 2 ;
	b = 3 ;
variables:
	double test(c, b) ;
		test:_FillValue = NaN ;
	int64 c(c) ;
	int64 b(b) ;
data:

 test =
  0, 0, 0,
  0, 0, 0 ;

 c = 0, 1 ;

 b = 0, 1, 2 ;
}

Now, from the dump, the 'c' dimension is first in both. Lets read those files again and concat them along the c-dimension:

dst_dim0 = xr.open_dataset('src_dim0.nc')
dst_dim0 = xr.concat([dst_dim0, dst_dim0], dim='c')
dst_dim0.to_netcdf('dst_dim0.nc')

dst_dim1 = xr.open_dataset('src_dim1.nc')
dst_dim1 = xr.concat([dst_dim1, dst_dim1], dim='c')
dst_dim1.to_netcdf('dst_dim1.nc')

Now, and this is what confuses me, the file without coordinates has 'c' dimension first and the file with coordinates has 'b' dimension first.:

netcdf dst_dim0 {
dimensions:
	c = 4 ;
	b = 3 ;
variables:
	double test(c, b) ;
		test:_FillValue = NaN ;
data:

 test =
  0, 0, 0,
  0, 0, 0,
  0, 0, 0,
  0, 0, 0 ;
}

netcdf dst_dim1 {
dimensions:
	b = 3 ;
	c = 4 ;
variables:
	int64 b(b) ;
	double test(c, b) ;
		test:_FillValue = NaN ;
	int64 c(c) ;
data:

 b = 0, 1, 2 ;

 test =
  0, 0, 0,
  0, 0, 0,
  0, 0, 0,
  0, 0, 0 ;

 c = 0, 1, 0, 1 ;
}

I really like to understand why there is this difference. Thanks for your patience!

[shoyer]

This is due to the internal implementation of xarray.concat, which sometimes reorders variables. I doubt the reordering was intentional. It's probably just a side effect of how concat takes multiple passes over different types of variables to figure out how to combine them.

You are welcome to take a look at improving this, though I doubt this would be particularly easy to fix. Certainly the code in concat could use some clean-up, and if we can preserve the order of variables on outputs that would be an improvement in usability. But it's still not something I would consider a "bug" per se.

[kmuehlbauer]

@shoyer Yes, that was what I was assuming. But was a bit confused too, as the concat docs say, that dimension order is not affected. But maybe I get this wrong and the order of dimensions is not affected only for DataArrays.

IIUC xarray creates a new dataset during concat, because the dimensions cannot be expanded (due to netCDF4 limitations). So I would need to look at that specific part, where this creation process takes place.

I would also not speak of "bug" here, but if such reordering happens only in certain conditions users (I mean at least me) can get confused.

I'll try to find out under what conditions this happens and try to come up with some workaround. Will also try ti find my way through the concat-mechanism. Again, I really appreciate your help in this issue.
I can rename the issue to a somewhat better name, do you have a suggestion? Ambiguous dimension reorder in Dataset concat maybe?

[kmuehlbauer]

Sorry, fat fingers...

[kmuehlbauer]

@shoyer I'm working on a notebook with all testing inside. Just found that if I have 3 dimensions ('c', 'd', 'b') the ordering is preserved in any case (with and without coords) for concat along any dimension. Will link the notebook next day.

Update: Need to be more thorough...with coordinates it reorders also with 3 dims.