Order of coordinates not preserved when calling `interp`

[schmoelder]

What happened?

When calling the interp method on DataArrays or Datasets, the order of the coordinates is not preserved.

What did you expect to happen?

Keep order of coordinates.

Minimal Complete Verifiable Example

import numpy as np
import xarray as xr

# %% DataArray
coordinates={
    "x": np.linspace(0, 1, 10),
    "y": np.linspace(0, 1, 5),
    "time": np.linspace(0, 1, 6),
}
grid_shape = tuple(len(v) for v in coordinates.values())

data = np.zeros(grid_shape)

data_array = xr.DataArray(data=data, coords=coordinates)

print(data_array.to_numpy().shape)
print(data_array.coords)

data_array_interp = data_array.interp(x=np.linspace(0, 1, 20))
print(data_array_interp.to_numpy().shape)
print(data_array_interp.coords)

MVCE confirmation

• Minimal example — the example is as focused as reasonably possible to demonstrate the underlying issue in xarray.
• Complete example — the example is self-contained, including all data and the text of any traceback.
• Verifiable example — the example copy & pastes into an IPython prompt or Binder notebook, returning the result.
• New issue — a search of GitHub Issues suggests this is not a duplicate.
• Recent environment — the issue occurs with the latest version of xarray and its dependencies.

Relevant log output

(10, 5, 6)
Coordinates:
  * x        (x) float64 80B 0.0 0.1111 0.2222 0.3333 ... 0.7778 0.8889 1.0
  * y        (y) float64 40B 0.0 0.25 0.5 0.75 1.0
  * time     (time) float64 48B 0.0 0.2 0.4 0.6 0.8 1.0
(20, 5, 6)
Coordinates:
  * y        (y) float64 40B 0.0 0.25 0.5 0.75 1.0
  * time     (time) float64 48B 0.0 0.2 0.4 0.6 0.8 1.0
  * x        (x) float64 160B 0.0 0.05263 0.1053 0.1579 ... 0.8947 0.9474 1.0

Anything else we need to know?

Note, the shape is as expected, but the order of the coordinates has been modified.

The same holds for datasets:

# %% Dataset
components = ["A", "B"]

coordinates={
    "x": np.linspace(0, 1, 10),
    "y": np.linspace(0, 1, 5),
    "time": np.linspace(0, 1, 6),
}
grid_shape = tuple(len(v) for v in coordinates.values())

data = np.zeros((len(components), *grid_shape))

data_vars = {
    comp: xr.DataArray(data=data[i, :], coords=coordinates)
    for i, comp in enumerate(components)
}
data_set = xr.Dataset(data_vars=data_vars, coords=coordinates)
print(data_set.to_array().shape)
print(data_set.coords)

data_set_interp = data_set.interp(x=np.linspace(0, 1, 20))
print(data_set_interp.to_array().shape)
print(data_set_interp.coords)

Output:

(2, 10, 5, 6)
Coordinates:
  * x        (x) float64 80B 0.0 0.1111 0.2222 0.3333 ... 0.7778 0.8889 1.0
  * y        (y) float64 40B 0.0 0.25 0.5 0.75 1.0
  * time     (time) float64 48B 0.0 0.2 0.4 0.6 0.8 1.0
(2, 20, 5, 6)
Coordinates:
  * y        (y) float64 40B 0.0 0.25 0.5 0.75 1.0
  * time     (time) float64 48B 0.0 0.2 0.4 0.6 0.8 1.0
  * x        (x) float64 160B 0.0 0.05263 0.1053 0.1579 ... 0.8947 0.9474 1.0
(10, 5, 6)
Coordinates:
  * x        (x) float64 80B 0.0 0.1111 0.2222 0.3333 ... 0.7778 0.8889 1.0
  * y        (y) float64 40B 0.0 0.25 0.5 0.75 1.0
  * time     (time) float64 48B 0.0 0.2 0.4 0.6 0.8 1.0
(20, 5, 6)
Coordinates:
  * y        (y) float64 40B 0.0 0.25 0.5 0.75 1.0
  * time     (time) float64 48B 0.0 0.2 0.4 0.6 0.8 1.0
  * x        (x) float64 160B 0.0 0.05263 0.1053 0.1579 ... 0.8947 0.9474 1.0

Environment

INSTALLED VERSIONS

commit: None
python: 3.12.11 | packaged by conda-forge | (main, Jun 4 2025, 14:45:31) [GCC 13.3.0]
python-bits: 64
OS: Linux
OS-release: 6.14.0-29-generic
machine: x86_64
processor: x86_64
byteorder: little
LC_ALL: None
LANG: de_DE.UTF-8
LOCALE: ('C', 'UTF-8')
libhdf5: 1.14.3
libnetcdf: None

xarray: 2025.8.0
pandas: 2.3.2
numpy: 2.2.6
scipy: 1.16.1
netCDF4: None
pydap: None
h5netcdf: 1.6.4
h5py: 3.13.0
zarr: None
cftime: None
nc_time_axis: None
iris: None
bottleneck: None
dask: None
distributed: None
matplotlib: 3.10.5
cartopy: None
seaborn: None
numbagg: None
fsspec: 2025.5.1
cupy: None
pint: None
sparse: None
flox: None
numpy_groupies: None
setuptools: 80.9.0
pip: 25.2
conda: None
pytest: 8.4.1
mypy: None
IPython: 8.37.0
sphinx: 8.3.0

[welcome]

Thanks for opening your first issue here at xarray! Be sure to follow the issue template!
If you have an idea for a solution, we would really welcome a Pull Request with proposed changes.
See the Contributing Guide for more.
It may take us a while to respond here, but we really value your contribution. Contributors like you help make xarray better.
Thank you!

[dcherian]

we'd happily take a PR to transpose back here:

xarray/xarray/core/missing.py

Line 765 in 7893348

return result

We'll need to account for the case where result.dims != variable.dims

[schmoelder]

Thanks for looking into it.

Another thing I noticed, for DataArrays, the order of the dims attribute is preserved:

print(data_array_interp.dims)
Out[49]: ('x', 'y', 'time')

but not for Datasets:

print(data_set_interp.dims)
Out[51]: FrozenMappingWarningOnValuesAccess({'y': 5, 'time': 6, 'x': 20})

So it's really hard for a dataset to actually check which array dimension belongs to which coordinate. Is there another way as a workaround?

Edit: for datasets, the _coord_names even has another order 😕

print(data_set_interp._coord_names)
Out[54]: {'time', 'x', 'y'}

[dcherian]

for DataArrays, the order of the dims attribute is preserved:

oh hmm.. great then.

but not for Datasets:

there is no defined dimension order on a dataset because it can contain arrays with different dimensions.

and yes _coord_names is a set so there is no order

[schmoelder]

and yes _coord_names is a set so there is no order

Ah, I missed that.

there is no defined dimension order on a dataset because it can contain arrays with different dimensions.

I understand. I'm still pretty new to xarray so I'm still missing a lot of intuition on the internals.

But my question remains: How can I programmatically infer which coordinate corresponds to which data dimension?

[jsignell]

I think the confusion is coming from the fact that coords just do not have an order. Dims have an order and some dims have a related coordinate. If a dimension has a coordinate then you can get that coordinate using dict-like indexing:

import xarray as xr

ds = xr.tutorial.open_dataset("air_temperature")

for dim in ds.air.dims:
    print(ds.coords[dim])

Or in your original example:

for data_array_interp.dims:
    print(data_array_interp[dim])

<xarray.DataArray 'x' (x: 20)> Size: 160B
array([0.      , 0.052632, 0.105263, 0.157895, 0.210526, 0.263158, 0.315789,
       0.368421, 0.421053, 0.473684, 0.526316, 0.578947, 0.631579, 0.684211,
       0.736842, 0.789474, 0.842105, 0.894737, 0.947368, 1.      ])
Coordinates:
  * x        (x) float64 160B 0.0 0.05263 0.1053 0.1579 ... 0.8947 0.9474 1.0
<xarray.DataArray 'y' (y: 5)> Size: 40B
array([0.  , 0.25, 0.5 , 0.75, 1.  ])
Coordinates:
  * y        (y) float64 40B 0.0 0.25 0.5 0.75 1.0
<xarray.DataArray 'time' (time: 6)> Size: 48B
array([0. , 0.2, 0.4, 0.6, 0.8, 1. ])
Coordinates:
  * time     (time) float64 48B 0.0 0.2 0.4 0.6 0.8 1.0

[jsignell]

I'm going to close this but if you have more questions or ideas for how to improve the docs please reopen or submit a pull request!

[jsignell]

Just came across #712 as well that I think this is essentially a duplicate of