Interpolation unstructured data to regular grid

[XavierTolza]

Is your feature request related to a problem?

Hello all. First, thanks for this amazing library!

I'm working with unstructured data ( points that are not in a grid ), and I want to interpolate them onto a regular grid.
However, it seems it is not directly supported by xarray and I have to manually call scipy's griddata.
This is annoying because:

• It's verbose and I'm lazy
• I loose my array's metadata

It seems the problem is encoutered by some other people as a google search on xarray interp irregular data yields to some results:

• https://stackoverflow.com/questions/60495410/add-irregular-data-lat-lons-to-xarray-dataset-for-interpolation
• https://gis.stackexchange.com/questions/455149/interpolate-irregularly-sampled-data-to-a-regular-grid

Searching for issues related to this topic led me to this single issue #5281 where the user uses directly scipy's griddata to interpolate.

Here is a sample script to reproduce and illustrate the problem:

from xarray import DataArray
import numpy as np
import matplotlib.pyplot as plt

# Sample 1000 uniform points in [-1,1]
# Thus, those points are unstructured
x, y = points = np.random.uniform(-1, 1, (2, 1000))
x, y = points = DataArray(
    points,
    coords=dict(
        variable=list("XY"), point=np.arange(1000), x=("point", x), y=("point", y)
    ),
    dims=("variable", "point"),
)

# Compute a dummy "function"
z = np.sqrt(x**2 + y**2)
print(z)
# We now have a function sampled on unstructured points
fig, ax = plt.subplots(1, 1)
ax.scatter(x, y, c=z)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_aspect(1)

plt.show()

I get the following value for z:

<xarray.DataArray (point: 1000)>
array([...blablabla whatever...])
Coordinates:
  * point    (point) int64 0 1 2 3 4 5 6 7 8 ... 992 993 994 995 996 997 998 999
    x        (point) float64 -0.4088 -0.7709 -0.7345 ... 0.5191 0.6061 -0.7329
    y        (point) float64 -0.8546 -0.6544 -0.6771 ... -0.8286 -0.04476 -0.574

And it looks like so:

So we can see z is correctly annotated and we've got the x and y value for each point, but I still cannot interpolate on a regular grid:

z.interp(x=np.linspace(-1, 1, 100), y=np.linspace(-1, 1, 100))

It raises the following error:

Dimensions {'x', 'y'} do not exist. Expected one or more of Frozen({'point': 1000})

Describe the solution you'd like

It would be nice to allow a direct unstructured interpolation if variables are correctly annotated (I supposed it's done with panda's multi indexes?)

This could be done using Scipy's griddata

A check could be added in interp function to see if the user wants to interpolate on sub-dimensions (i'm not sure if the name is appropriate), and if so, using griddata instead of the other interpolation methods.

Describe alternatives you've considered

Alternative solutions could be to use only scipy grid data by hand, but that would be sad.

Additional context

No response

[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!