text: Title: Advanced interpolation returning unexpected shape
  Body: ### What happened?
  
  Hi !
  
  Iam using Xarray since quite some time and recently I started using fresh NC files from https://cds.climate.copernicus.eu/ (up to now, using NC4 files that I have from quite some time now). I have an issue regarding Xarray' interpolation with advanced indexing: the returned shape by interp() is not as expected.
  
  basically Iam just doing advanced interpolation as mentioned here https://docs.xarray.dev/en/stable/user-guide/interpolation.html#advanced-interpolation.
  It used to work perfectly with my older datasets.
  
  thanks
  
  Vianney
  
  
  
  ### What did you expect to happen?
  
  I would expect res.shape to be (100,) or similar (and actually I get this result with older datasets. but I don't understand why I get (705, 715, 2) ! this does not corrrespond to anything, as we can see from the below description of the dataset:
  
  
  ### Minimal Complete Verifiable Example
  [dataset.zip](https://github.com/user-attachments/files/17960240/dataset.zip)
  ```Python
  # -*- coding: utf-8 -*-
  
  import dask
  import numpy as np
  import xarray as xa
  import datetime as dt
  
  # open dataset (unzip dataset.zip which is attached to bug report)
  dataset = xa.open_mfdataset("dataset.nc", decode_times=True, parallel=True)
  start_date = dt.datetime(2023, 1, 1)
  
  np.random.seed(0)
  
  # points generation
  N = 100
  access_date = np.linspace(0, 5000, N) # seconds since start_date
  longitudes = -180 + np.random.rand(N)*360
  latitudes = -90 + np.random.rand(N)*180
  
  dates_folded = np.datetime64(start_date) + \
      access_date.astype("timedelta64[s]").astype("timedelta64[ns]")
  
  arguments = {
      "valid_time": xa.DataArray(dates_folded, dims="z"),
      "longitude": xa.DataArray(longitudes, dims="z"),
      "latitude": xa.DataArray(latitudes, dims="z"),
      "method": "linear",
      "kwargs": {"bounds_error": False},
  }
  with dask.config.set(**{"array.slicing.split_large_chunks": True}):
      data = dataset['tcc'].interp(**arguments)
      
  res = np.array(data).T
  print(res.shape)
  ```
  
  Output: 
  ```
  (705, 715, 2)
  ```
  
  
  ### MVCE confirmation
  
  - [X] Minimal example — the example is as focused as reasonably possible to demonstrate the underlying issue in xarray.
  - [X] Complete example — the example is self-contained, including all data and the text of any traceback.
  - [X] Verifiable example — the example copy & pastes into an IPython prompt or [Binder notebook](https://mybinder.org/v2/gh/pydata/xarray/main?urlpath=lab/tree/doc/examples/blank_template.ipynb), returning the result.
  - [X] New issue — a search of GitHub Issues suggests this is not a duplicate.
  - [X] Recent environment — the issue occurs with the latest version of xarray and its dependencies.
  
  ### Relevant log output
  
  ```Python
  (705, 715, 2)
  ```
  
  
  ### Anything else we need to know?
  
  the expected log is (100,)
  
  ### Environment
  
  <details>
  
  INSTALLED VERSIONS
  ------------------
  commit: None
  python: 3.11.9 | packaged by conda-forge | (main, Apr 19 2024, 18:27:10) [MSC v.1938 64 bit (AMD64)]
  python-bits: 64
  OS: Windows
  OS-release: 10
  machine: AMD64
  processor: Intel64 Family 6 Model 140 Stepping 1, GenuineIntel
  byteorder: little
  LC_ALL: None
  LANG: fr
  LOCALE: ('fr_FR', 'cp1252')
  libhdf5: 1.14.3
  libnetcdf: 4.9.2
  
  xarray: 2024.11.0
  pandas: 2.2.2
  numpy: 1.26.4
  scipy: 1.14.0
  netCDF4: 1.6.5
  pydap: None
  h5netcdf: None
  h5py: 3.11.0
  zarr: None
  cftime: 1.6.4
  nc_time_axis: None
  iris: None
  bottleneck: None
  dask: 2024.7.0
  distributed: 2024.7.0
  matplotlib: 3.9.1
  cartopy: 0.23.0
  seaborn: None
  numbagg: None
  fsspec: 2024.6.1
  cupy: 13.2.0
  pint: None
  sparse: None
  flox: None
  numpy_groupies: None
  setuptools: 70.3.0
  pip: 24.0
  conda: 24.5.0
  pytest: 8.2.2
  mypy: None
  IPython: 8.26.0
  sphinx: 7.4.0
  
  </details>
  