Xarray: Is it possible to perform this interpolation with xarray?

[aragong]

I'm trying to interpolate information from a 3D dataset (lon,lat,time) ussing directly xarray.

When I made a simply interpolation with only one point I have no problem at all.

lat = [44.25]
lon = [-4.5]
t = datetime.strptime('2000-02-28 01:00:00', '%Y-%m-%d %H:%M:%S')

ds = xr.open_dataset('file.nc')
vx = ds['uo_surface'].interp(longitude=lon, latitude=lat, time=t)
But now I'm trying to interpolate in the same way several points and the result of this operation following the same syntax shows more results of what I will expected.

lat = [44.25, 45.25]
lon = [-4.5, -5]
t = datetime.strptime('2000-02-28 01:00:00', '%Y-%m-%d %H:%M:%S')

ds = xr.open_dataset('Currents\oceanTESEO.nc')
vx = ds['uo_surface'].interp(longitude=lon, latitude=lat, time=[t, t])
The result is this array:

array([[[0.01750018, 0.05349977],
[0.03699994, 0.11299999]],
[[0.01750018, 0.05349977],
[0.03699994, 0.11299999]]])

However, I expect only 2 values, one for each (lon,lat,t) point. Do I have to implement a loop to do that? I suposse this feature is already included in xarray. Do you know other way to calculate this sort of point interpolation faster and with 4D datarrays (lon,lat,z,time)?

Thank you in advance!!!
https://stackoverflow.com/questions/56144678/interpolation-syntax

[fmaussion]

Yes, it is possible. It is a bit "less intuitive" at first sight, but powerful and documented here: http://xarray.pydata.org/en/stable/interpolation.html#advanced-interpolation

The call you need to make is:

blah.interp(longitude=('z', lon), latitude=('z', lat))

[aragong]

Yes, it is possible. It is a bit "less intuitive" at first sight, but powerful and documented here: http://xarray.pydata.org/en/stable/interpolation.html#advanced-interpolation

The call you need to make is:

blah.interp(longitude=('z', lon), latitude=('z', lat))

Thank you so much! It works fine! I guess that we are creating a new one common dimension with only this points to interpolate the data. I did this:

lat = [44.25, 45.25]
lon = [-4.5, -5]
t = datetime.strptime('2000-02-28 01:00:00', '%Y-%m-%d %H:%M:%S')
vx = da.interp(longitude=('p', lon), latitude=('p', lat), time=('p', [t, t]))

But I can't figure out what was the code doing without create this new one "common dimension". ¿Do you have any clue about that?¿is the code making a subset interpolation between the coordinates?

[fmaussion]

But I can't figure out what was the code doing without create this new one "common dimension"

It is doing what is called "orthogonal indexing", but with interpolation. The resulting shape of the output is then (2, 2, 2) in your case, but could be any (t, y, x) as given by the size of each dimension indexer.

Maybe this helps a little: http://xarray.pydata.org/en/stable/indexing.html#vectorized-indexing

[aragong]

Thank you @fmaussion, I will review that link!

[fmaussion]

Thanks! If you want you can also accept the stackoverflow answer for later reference