Constructing a 2d interpolator given scattered input data

Question

I have three lists as follows:

x = [100.0, 75.0, 50.0, 0.0, 100.0, 75.0, 50.0, 0.0, 100.0, 75.0, 50.0, 0.0, 100.0, 75.0, 50.0, 0.0, 100.0, 75.0, 50.0, 0.0, 100.0, 75.0, 50.0, 0.0, 100.0, 75.0, 50.0, 0.0, 100.0, 75.0, 50.0, 0.0, 100.0, 75.0, 50.0, 0.0, 100.0, 75.0, 50.0, 0.0]

y = [300.0, 300.0, 300.0, 300.0, 500.0, 500.0, 500.0, 500.0, 700.0, 700.0, 700.0, 700.0, 1000.0, 1000.0, 1000.0, 1000.0, 1500.0, 1500.0, 1500.0, 1500.0, 2000.0, 2000.0, 2000.0, 2000.0, 3000.0, 3000.0, 3000.0, 3000.0, 5000.0, 5000.0, 5000.0, 5000.0, 7500.0, 7500.0, 7500.0, 75000.0, 10000.0, 10000.0, 10000.0, 10000.0]

z = [100.0, 95.0, 87.5, 77.5, 60.0, 57.0, 52.5, 46.5, 40.0, 38.0, 35.0, 31.0, 30.0, 28.5, 26.25, 23.25, 23.0, 21.85, 20.125, 17.825, 17.0, 16.15, 14.875, 13.175, 13.0, 12.35, 11.375, 10.075, 10.0, 9.5, 8.75, 7.75, 7.0, 6.65, 6.125, 5.425, 5.0, 4.75, 4.375, 3.875]

Each entry of each list is read as a point so point 0 is (100,300,100) point 1 is (75,300,95) and so on.

I am trying to do 2d interpolation, so that I can compute a z value for any given input (x0, y0) point.

I was reading that using meshgrid I can interpolate with RegularGridInterpolator from scipy but I am not sure how to set it up when I do:

x_,y_,z_ = np.meshgrid(x,y,z) # both indexing ij or xy

I don't get values for x_,y_,z_ that make sense and I am not sure how to go from there.

I am trying to use the data points I have above to find intermediate values so something similar to scipy's interp1d where

f = interp1d(x, y, kind='cubic')

where I can later call f(any (x,y) point within range) and get the corresponding z value.

Answer 1

You need 2d interpolation over scattered data. I'd default to using scipy.interpolate.griddata in this case, but you seem to want a callable interpolator, whereas griddata needs a given set of points onto which it will interpolate.

Not to worry: griddata with 2d cubic interpolation uses a CloughTocher2DInterpolator. So we can do exactly that:

import numpy as np
import scipy.interpolate as interp

x = [100.0, 75.0, 50.0, 0.0, 100.0, 75.0, 50.0, 0.0, 100.0, 75.0, 50.0, 0.0, 100.0, 75.0, 50.0, 0.0, 100.0, 75.0, 50.0, 0.0, 100.0, 75.0, 50.0, 0.0, 100.0, 75.0, 50.0, 0.0, 100.0, 75.0, 50.0, 0.0, 100.0, 75.0, 50.0, 0.0, 100.0, 75.0, 50.0, 0.0]
y = [300.0, 300.0, 300.0, 300.0, 500.0, 500.0, 500.0, 500.0, 700.0, 700.0, 700.0, 700.0, 1000.0, 1000.0, 1000.0, 1000.0, 1500.0, 1500.0, 1500.0, 1500.0, 2000.0, 2000.0, 2000.0, 2000.0, 3000.0, 3000.0, 3000.0, 3000.0, 5000.0, 5000.0, 5000.0, 5000.0, 7500.0, 7500.0, 7500.0, 75000.0, 10000.0, 10000.0, 10000.0, 10000.0]
z = [100.0, 95.0, 87.5, 77.5, 60.0, 57.0, 52.5, 46.5, 40.0, 38.0, 35.0, 31.0, 30.0, 28.5, 26.25, 23.25, 23.0, 21.85, 20.125, 17.825, 17.0, 16.15, 14.875, 13.175, 13.0, 12.35, 11.375, 10.075, 10.0, 9.5, 8.75, 7.75, 7.0, 6.65, 6.125, 5.425, 5.0, 4.75, 4.375, 3.875]

interpolator = interp.CloughTocher2DInterpolator(np.array([x,y]).T, z)

Now you can call this interpolator with 2 coordinates to give you the corresponding interpolated data point:

>>> interpolator(x[10], y[10]) == z[10]
True
>>> interpolator(2, 300)
array(77.81343)

Note that you'll have to stay inside the convex hull of the input points, otherwise you'll get nan (or whatever is passed as the fill_value keyword to the interpolator):

>>> interpolator(2, 30)
array(nan)

Extrapolation is usually meaningless anyway, and your input points are scattered in a bit erratic way:

So even if extrapolation was possible I wouldn't believe it.

---

Just to demonstrate how the resulting interpolator is constrained to the convex hull of the input points, here's a surface plot of your data on a gridded mesh we create just for plotting:

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# go linearly in the x grid
xline = np.linspace(min(x), max(x), 30)
# go logarithmically in the y grid (considering y distribution)
yline = np.logspace(np.log10(min(y)), np.log10(max(y)), 30)
# construct 2d grid from these
xgrid,ygrid = np.meshgrid(xline, yline)
# interpolate z data; same shape as xgrid and ygrid
z_interp = interpolator(xgrid, ygrid)

# create 3d Axes and plot surface and base points
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(xgrid, ygrid, z_interp, cmap='viridis',
                vmin=min(z), vmax=max(z))
ax.plot(x, y, z, 'ro')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')

plt.show()

Here's the output from two angles (it's better to rotate around interactively; such stills don't do the 3d representation justice):

There are two main features to note:

1. The surface nicely fits the red points, which is expected from interpolation. Fortunately the input points are nice and smooth so everything goes well with interpolation. (The fact that the red points are usually hidden by the surface is only due to how pyplot's renderer mishandles the relative position of complex 3d objects)
2. The surface is cut (due to nan values) along the convex hull of the input points, so even though our gridded arrays define a rectangular grid we only get a cut of the surface where interpolation makes sense.