SciPy LinearNDInterpolator RegularGridInterpolator produces different result for the same data

I am experimenting with the LinearNDInterpolator and RegularGridInterpolator functions from SciPy to interpolate 1D, 2D and 3D data linearly.

I understand that the main difference between the two methods is that the RegularGridInterpolator fis designed to handle data on a rectilinear (or "regular") grid with either even or uneven spacing. This function is also more efficient compared to LinearNDInterpolator.

However, when I generate random data on a rectilinear grid with even spacings, the results obtained from these two methods don't match. I was expecting them to be similar since the grid is regular and I'm using linear interpolation in both cases.

Here is the code I used to compare these two methods:

import numpy as np
from scipy.interpolate import LinearNDInterpolator
from scipy.interpolate import RegularGridInterpolator
import matplotlib.pyplot as plt

#original data
x = np.linspace(0, 1, num=20)
y = np.linspace(1, 2, num=10)
X, Y = np.meshgrid(x, y)
values = np.random.rand(20, 10)
points = np.column_stack((X.flatten(), Y.flatten()))
values_flat = values.flatten()

#LinearNDInterpolation
interfunc = LinearNDInterpolator(points, values_flat)
x1 = np.linspace(0, 1, num=3000)
y1 = np.linspace(1, 2, num=3000)
X1, Y1 = np.meshgrid(x1, y1)
interpolated_values = interfunc(np.column_stack((X1.flatten(), Y1.flatten())))
interpolated_values = interpolated_values.reshape(X1.shape)
fig, ax = plt.subplots()
linear = ax.contourf(X1, Y1, interpolated_values.T)
fig.colorbar(linear, ax=ax)

#RegularGridInterpolation
fig2, ax2 = plt.subplots()
x2 = np.linspace(0, 1, num=3000)
y2 = np.linspace(1, 2, num=3000)
X2, Y2 = np.meshgrid(x2, y2)
points_grid = (x, y)
interfunc_grid = RegularGridInterpolator(points_grid, values, method="linear")
interpolated_values_grid = interfunc_grid(np.column_stack((X2.flatten(), Y2.flatten())))
interpolated_values_grid = interpolated_values_grid.reshape(X2.shape)
d = ax2.contourf(X2, Y2, interpolated_values_grid.T)
fig2.colorbar(d, ax=ax2)
plt.show()

The code generates two contour plots: one for the LinearNDInterpolator and the other for the RegularGridInterpolator. The plots show the interpolated values over a finer grid.

Am I making a mistake somewhere in my code or is there another reason why the results from these two methods of interpolation would differ? Any insights would be greatly appreciated.

LinearNDIbterpolator is based on triangulation (via QHull). Each interpolated value is a linear combination of exactly three points. And the choice of these three is degenerate for data on a grid.

Bottom line : do not use LinearNDinterpolator for data on a grid.