Evaluate a RegularGridInterpolator on an another regular grid

Question

I'm having trouble understanding how to shape data to evaluate an interpolated view of an nD-array, using scipy.interpolate.RegularGridInterpolator

Considering A a (n1,n2,n3)-shaped numpy array, indexed along the following coordinates :

x = np.linspace(0, 10, 5) # n1 = 5
y = np.linspace(-1, 1, 10) # n2 = 10
z = np.linspace(0, 500, 1000) # n3 = 1000

For this example, you can generate A = ex_array with this bit of code from the documentation :

def f(x,y,z):
    return 2 * x**3 + 3 * y**2 - z

ex_array = f(*np.meshgrid(x, y, z, indexing='ij', sparse=True))

Let's imagine I want to interpolate the entire array along each axis. This is done with :

from scipy.interpolate import RegularGridInterpolator

interpolated = RegularGridInterpolator((x,y,z), ex_array)

Now, the part where my brain starts to hurt hard : In order to evaluate this interpolator object at any given coordinates, you have to __call__ it on said point like so :

evaluated_pts = interpolated((0,1,0)) # evaluate at (x,y,z) = (5,0.5,300)
print(evaluated_pts)

In order to evaluate it on several points, you can iterate like this :

pts = ((5,0.5,_z) for _z in np.linspace(100,200,50))
evaluated_pts = interpolated(pts)

Now, what if I want to use the same logic as above, and evaluate on an entire new grid, such as :

new_x = np.linspace(2, 3, 128)
new_y = np.linspace(-0.1, 0.1, 100)
new_z = np.linspace(350, 400, 256)

As you can see now, it's not as straightforward as interpolated(new_x, new_y, new_z), and I tried to use np.meshgrid but could not figure it out.

Ideally, I'd want to output a new (128, 100, 256) array in this example.

I have the same question. As I see it, `scipy.griddata` is useful to map *unstructured* data, i.e. a cloud of points to another cloud of points. `RGD` can map *structured* to *unstructured*. **HOWEVER** there isn't any method to map from *structured* to *structured*. The best you can do is pass a cloud of points into RGD, which seems very inefficient to me. But maybe it's no big deal :) — John Henckel, Dec 02 '21 at 17:34

score 2 · Answer 1 · answered Aug 27 '19 at 12:33

RegularGridInterpolator input values are located on a grid. The grid points are defined using a tuple of "ticks" along each axis, for instance ((x0, x1, ..., xn), (y0, y1, ..., xm), (z0, z1, ..., zk) ) in 3D. The values are given as an nd-array of shape (n, m, k) in this case.

To evaluate the interpolated function, the assumption that the points are on a grid is no more required. Then, the asked points are defined as a list of points (actually an array of coordinates): ((x1, y1, z1), (x2, y2, z2), ... (xP, yP, zP)) i.e. a nd-array of shape (Number of points, Number of dimension).

To evaluate the interpolation on a new grid, it must be constructed using meshgrid.

reshape and transpose are used to transform arrays from one shape to another (see this question).

For example:

x = [0, 1, 2]
y = [3, 4]
z = [5, 6, 7, 8]

xyz_grid = np.meshgrid(x, y, z, indexing='ij')

xyz_list = np.reshape(xyz_grid, (3, -1), order='C').T

xyz_list

̀xyz_list could be used to call the interpolation function and it looks like that:

array([[0, 3, 5],
       [0, 3, 6],
       [0, 3, 7],
       [0, 3, 8],
       [0, 4, 5],
       [0, 4, 6],
       [0, 4, 7],
       [0, 4, 8],
       [1, 3, 5],
       [1, 3, 6],
       [1, 3, 7],
       [1, 3, 8],
       [1, 4, 5],
       [1, 4, 6],
       [1, 4, 7],
       [1, 4, 8],
       [2, 3, 5],
       [2, 3, 6],
       [2, 3, 7],
       [2, 3, 8],
       [2, 4, 5],
       [2, 4, 6],
       [2, 4, 7],
       [2, 4, 8]])

Thank you for this answer. I ended up implementing a solution which is similar to yours in its logic but not in implementation. I'm going to post it later here as own-answer. I'll accept yours if you match it to my example data and finally output the required array ? Thank you — MaximGi, Aug 27 '19 at 14:47

Evaluate a RegularGridInterpolator on an another regular grid

1 Answers1