After checking a bit more in older issues, this seems related: https://github.com/pydata/xarray/pull/1049, ping @fmaussion.

And also @shoyer's comment suggest that those two iterations/loops I mentioned above need to be addressed correctly.

@shoyer Attached the description of the issue source and kind of workaround.

During concat a result_vars = OrderedDict() is created. After that it is iterated over the first dataset datasets[0].variables.items() and those variables which are not affected by the concat are added to the result_vars : https://github.com/pydata/xarray/blob/a5ca64ac5988f0c9c9c6b741a5de16e81b90cad5/xarray/core/combine.py#L244-L246

After several checks the affected variables are treated and added to result_vars: https://github.com/pydata/xarray/blob/a5ca64ac5988f0c9c9c6b741a5de16e81b90cad5/xarray/core/combine.py#L301-L306

The comment indicates what you already mentioned, that the reorder might be unintentional. But due to the handling in two separate iterations over datasets[0].variables, the source variable order is not preserved (and with that in some cases the order of the dimensions).

This can be worked around by changing the second iteration to: ```python

re-initialize result_vars to write in correct order

result_vars = OrderedDict()

stack up each variable to fill-out the dataset (in order)

for k in datasets[0].variables: if k in concat_over: vars = ensure_common_dims([ds.variables[k] for ds in datasets]) combined = concat_vars(vars, dim, positions) insert_result_variable(k, combined) else: insert_result_variable(k, datasets[0].variables[k])
`` With this workaround applied, theconcat` works as expected and the variable/coordinate order (and with that the dimension order) is preserved. I'm thinking about a better solution but wanted to get some feedback from you first, if I#m on the right track. Thanks!

@shoyer I think I found the relevant lines of code in combine.py. It might be possible to preserve the (correct) order of dimensions at least with regard to their occurrence in variables (and not in coordinates). But that would mean to treat variables and coordinates consecutively in some way..

In the docs there is a Warning: We are changing the behavior of iterating over a Dataset the next major release of xarray, to only include data variables instead of both data variables and coordinates. In the meantime, prefer iterating over ds.data_vars or ds.coords. below here.

Does that mean that this also affects internal machinery (like in concat)? If so, could you point me to some code where this is taken care of or give some explanation or links where this is discussed?

Update: I' working with latest 0.12.0 release.

Just as note for me, to not have to reiterate:

• It seems, that variables are handled in creation order. Means that concat reads them, handles them (even if the variable remains unchanged) and writes them to the new dataset in that order.
• This does not happen for coordinate for some reason. There only the changed dimension is handled and written after the variables. The unaffected coordinates are written at the beginning before the variables.

Example (dst concat over 'x'): The ordering of the src_dim1 is because the dimensions in the variables/coordinates are x,y,z in that order. The ordering of the dst_dim1 is because the dimensions in the variables/coordinates are z, y, x. netcdf src_dim1 { dimensions: x = 2 ; y = 3 ; z = 4 ; variables: double test2(x, y) ; test2:_FillValue = NaN ; double test3(x, z) ; test3:_FillValue = NaN ; double test1(y, z) ; test1:_FillValue = NaN ; int64 z(z) ; int64 y(y) ; int64 x(x) ; netcdf dst_dim1 { dimensions: z = 4 ; y = 3 ; x = 4 ; variables: int64 z(z) ; int64 y(y) ; double test2(x, y) ; test2:_FillValue = NaN ; double test3(x, z) ; test3:_FillValue = NaN ; double test1(x, y, z) ; test1:_FillValue = NaN ; int64 x(x) ;

It seems, that the two coordinates (z and y) are written first, then the variables, and then the changed coordinate. Now trying to find, where this happens. If the two coordinates would be written in the same way as the variables (and after them), then the ordering would be x,y,z as in the source.

@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.

Sorry, fat fingers...

@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?

@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.

```python 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:

```python 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!

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.

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:

python 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.

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?

Xref: Gitter Chat