How to get a non-smoothing 2D spline interpolation with scipy
I want a 2D cubic spline fit to some irregullary spaced data - i.e. a function that exactly fits the data at the given points - but can also return values in between.
All I can find (for irregural spaced data) is scipy.interpolate.SmoothBivariateSpline. I can't figure out how to turn 'smoothing' off (no matter what value I put in the s parameter.
I did, however, find I can get mostly what I want with scipy.interpolate.griddata - though this has to recalculate it every time (i.e. doesn't just generate a function). Is there any difference, fundamentally between these two - i.e. is griddata doing something different than a 'spline'? Is there anyway to turn off smoothing in the SmoothBivariateSpline or an equivalent function that doesn't smooth?
The following is a script I'm using to test fitting of a spline vs a polynomial
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import scipy.optimize
import scipy.interpolate
import matplotlib.pyplot as plt
import numpy.polynomial.polynomial as poly
# Grid and test function
N = 9;
x,y = np.linspace(-1,1, N), np.linspace(-1,1, N)
X,Y = np.meshgrid(x,y)
F = lambda X,Y : X+Y-1*X*Y-(X*Y)**2 -2*X*Y**2 + X**2*Y + 3*np.exp(-((X+1)**2+(Y+1)**2)*5)
Z = F(X,Y)
noise = 0.4
Z *= 1+(np.random.random(Z.shape)*2-1)*noise # noise
# Finer Grid and test function
N2 = 19;
x2,y2 = np.linspace(-1,1, N2), np.linspace(-1,1, N2)
X2,Y2 = np.meshgrid(x2,y2)
Z2 = F(X2,Y2)
# Make data into lists
Xl = X.reshape(X.size)
Yl = Y.reshape(Y.size)
Zl = Z.reshape(Z.size)
# Polynomial fit
# polyval(x,y,p) = p[0,0]+p[0,1]y+p[1,0]x+p[1,1]xy+p[1,2]xy^2 ..., etc
# I use a flat (1D) array for p, so it needs to be reshaped into a 2D array before
# passing to polyval
order = 3
p0 = np.zeros(order**2) # guess parameters (all 0 for now)
f_poly = lambda x,y,p : poly.polyval2d(x,y,p.reshape((order,order))) # Wrapper for our polynomial
errf = lambda p : np.mean((f_poly(Xl,Yl,p.reshape((order,order)))-Zl)**2) # error function to find least square error
sol = scipy.optimize.minimize(errf, p0)
psol = sol['x']
# Spline interpolation
# Bivariate (2D), Smoothed (doesn't fit points *exactly*) cubic (3rd order - i.e. kx=ky=3) spline
spl = scipy.interpolate.SmoothBivariateSpline(Xl, Yl, Zl, kx=3,ky=3)
f_spline = spl.ev
# regular Interpolate
f_interp = lambda x,y : scipy.interpolate.griddata((Xl, Yl), Zl, (x,y), method='cubic')
# Plot
fig = plt.figure(1, figsize=(7,8))
plt.clf()
# poly fit
ax = fig.add_subplot(311, projection='3d')
ax.scatter3D(X2,Y2,Z2,s=3, color='red', label='actual data')
fit = f_poly(X2,Y2, psol)
l = 'order {} poly fit'.format(order)
ax.plot_wireframe(X2,Y2, fit, color='black', label=l)
ax.scatter3D(X,Y,Z, color='blue', label='noisy data')
plt.legend()
print("Average {} error: {}".format(l, np.sqrt(np.mean((fit-Z2)**2))))
# spline fit
ax = fig.add_subplot(312, projection='3d')
ax.scatter3D(X2,Y2,Z2,s=3, color='red', label='actual data')
l = 'smoothed spline'
fit = f_spline(X2,Y2)
ax.plot_wireframe(X2,Y2, fit, color='black', label=l)
ax.scatter3D(X,Y,Z, color='blue', label='noisy data')
plt.legend()
print("Average {} error: {}".format(l, np.sqrt(np.mean((fit-Z2)**2))))
# interp fit
ax = fig.add_subplot(313, projection='3d')
ax.scatter3D(X2,Y2,Z2,s=3, color='red', label='actual data')
l='3rd order interp '
fit=f_interp(X2,Y2)
ax.plot_wireframe(X2,Y2, fit, color='black', label=l)
ax.scatter3D(X,Y,Z, color='blue', label='noisy data')
plt.legend()
print("Average {} error: {}".format(l, np.sqrt(np.mean((fit-Z2)**2))))
plt.show(False)
plt.pause(1)
raw_input('press key to continue') # Change to input() if using python3
For unstructured mesh, griddata is the right interpolation tool. However, the triangulation (Delaunay) and the interpolation is performed each time. One workaround is to use either CloughTocher2DInterpolator for a C1 smooth interpolation or LinearNDInterpolator for a linear interpolation. These are the functions actually used by griddata. The difference is that it is possible to use as input a Delaunay object and it returns an interpolation function.
Here is an example based on your code:
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from scipy.interpolate import CloughTocher2DInterpolator
from scipy.spatial import Delaunay
# Example unstructured mesh:
nodes = np.array([[-1. , -1. ],
[ 1. , -1. ],
[ 1. , 1. ],
[-1. , 1. ],
[ 0. , 0. ],
[-1. , 0. ],
[ 0. , -1. ],
[-0.5 , 0. ],
[ 0. , 1. ],
[-0.75 , 0.4 ],
[-0.5 , 1. ],
[-1. , -0.6 ],
[-0.25 , -0.5 ],
[-0.5 , -1. ],
[-0.20833333, 0.5 ],
[ 1. , 0. ],
[ 0.5 , 1. ],
[ 0.36174242, 0.44412879],
[ 0.5 , -0.03786566],
[ 0.2927264 , -0.5411368 ],
[ 0.5 , -1. ],
[ 1. , 0.5 ],
[ 1. , -0.5 ]])
# Theoretical function:
def F(x, y):
return x + y - x*y - (x*y)**2 - 2*x*y**2 + x**2*y + 3*np.exp( -((x+1)**2 + (y+1)**2)*5 )
z = F(nodes[:, 0], nodes[:, 1])
# Finer regular grid:
N2 = 19
x2, y2 = np.linspace(-1, 1, N2), np.linspace(-1, 1, N2)
X2, Y2 = np.meshgrid(x2, y2)
# Interpolation:
tri = Delaunay(nodes)
CT_interpolator = CloughTocher2DInterpolator(tri, z)
z_interpolated = CT_interpolator(X2, Y2)
# Plot
fig = plt.figure(1, figsize=(8,14))
ax = fig.add_subplot(311, projection='3d')
ax.scatter3D(nodes[:, 0], nodes[:, 1], z, s=15, color='red', label='points')
ax.plot_wireframe(X2, Y2, z_interpolated, color='black', label='interpolated')
plt.legend();
The obtained graph is:
Both the spline method and Clough-Tocher interpolation are based on constructing a piecewise polynomial function on mesh elements. The difference is that for the spline the mesh is regular and given by the algorithm (see .get_knots()). And the coefficients are fitted so that the function is close as possible to the points and smooth (fit). For the Clough-Tocher interpolation, the mesh elements are the ones given as input. The resulting function is therefore guaranteed to go through the points.