Interpolation with scipy



In [ ]:

    
%matplotlib inline

Official documentation: https://docs.scipy.org/doc/scipy/reference/interpolate.html

Import modules and initialize data



In [ ]:

    
import numpy as np
import pandas as pd
import scipy.interpolate
import matplotlib.pyplot as plt



In [ ]:

    
from mpl_toolkits.mplot3d import axes3d

Interpolate 1D functions

In the following examples, we interpolate $f(x) \mapsto \sin(x)$



In [ ]:

    
xmin, xmax = 0., 4*np.pi

x = np.linspace(xmin, xmax, 10)
y = np.sin(x)

x2 = np.linspace(xmin, xmax, 100)

Linear interpolation



In [ ]:

    
# Linear interpolation with extrapolation
f = scipy.interpolate.interp1d(x, y,
                               kind='linear',
                               fill_value="extrapolate")

y2 = f(x2)

plt.plot(x, y, "o:b", label="original")
plt.plot(x2, y2, ".-r", label="interpolated")
plt.legend();

B-Splines interpolation

https://docs.scipy.org/doc/scipy/reference/interpolate.html#d-splines



In [ ]:

    
# https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.splrep.html#scipy.interpolate.splrep
# https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.splprep.html#scipy.interpolate.splprep

spl = scipy.interpolate.splrep(x, y)

y2 = scipy.interpolate.splev(x2, spl)

plt.plot(x, y, "o:b", label="original")
plt.plot(x2, y2, ".-r", label="interpolated")
plt.legend();



In [ ]:

    
spl = scipy.interpolate.splrep(x, y,
                               xb=x[0], xe=x[-1],   # The interval to fit
                               #s=0.,               # A smoothing factor
                               k=1)                 # The degree fo the spline fit

y2 = scipy.interpolate.splev(x2, spl)

plt.plot(x, y, "o:b", label="original")
plt.plot(x2, y2, ".-r", label="interpolated")
plt.legend();

Spline linear interpolation



In [ ]:

    
# Spline linear interpolation with extrapolation (should be the same than spline1...)
f = scipy.interpolate.interp1d(x, y,
                               kind='slinear',
                               fill_value="extrapolate")

y2 = f(x2)

plt.plot(x, y, "o:b", label="original")
plt.plot(x2, y2, ".-r", label="interpolated")
plt.legend();

Interpolate 2D functions

In the following examples, we interpolate $f(x, y) \mapsto \sin(x) + \sin(y)$



In [ ]:

    
# Build data

x = np.arange(-1*np.pi, 1*np.pi, np.pi/4)
y = np.arange(-1*np.pi, 1*np.pi, np.pi/4)

xx, yy = np.meshgrid(x, y)

z = np.sin(xx) + np.sin(yy)



In [ ]:

    
# Plot data

fig = plt.figure(figsize=(12, 8))
ax = axes3d.Axes3D(fig)

#ax.plot_wireframe(xx, yy, z)

surf = ax.plot_surface(xx,
                       yy,
                       z,
                       cmap='gnuplot2', # 'jet' # 'gnuplot2'
                       rstride=1,
                       cstride=1,
                       shade=False)

plt.title("Original data")
plt.show();

Linear interpolation

Documentation: https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp2d.html#scipy.interpolate.interp2d



In [ ]:

    
f = scipy.interpolate.interp2d(x, y, z, kind='linear',
                               bounds_error=True)     # Let 'f' raise an exception when the requested point is outside the range defined by x and y



In [ ]:

    
# Build data

x_hd = np.arange(-1*np.pi, 1*np.pi-np.pi/4, np.pi/32)
y_hd = np.arange(-1*np.pi, 1*np.pi-np.pi/4, np.pi/32)

xx_hd,yy_hd = np.meshgrid(x_hd, y_hd)

z_hd = np.zeros(xx_hd.shape)

for xi in range(z_hd.shape[0]):
    for yi in range(z_hd.shape[1]):
        z_hd[xi, yi] = f(x_hd[xi], y_hd[yi])



In [ ]:

    
# Plot data

fig = plt.figure(figsize=(12, 8))
ax = axes3d.Axes3D(fig)

surf = ax.plot_surface(xx_hd,
                       yy_hd,
                       z_hd,
                       cmap='gnuplot2', # 'jet' # 'gnuplot2'
                       rstride=1,
                       cstride=1,
                       shade=False)

plt.show();

Non uniform grid data

Documentation: https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp2d.html#scipy.interpolate.interp2d



In [ ]:

    
# Build data

x_nu = np.arange(-1*np.pi, 1*np.pi, np.pi/4)
y_nu = np.arange(-1*np.pi, 1*np.pi, np.pi/4)

x_nu = x_nu.tolist()
y_nu = y_nu.tolist()

del x_nu[2]
del y_nu[2]

xx, yy = np.meshgrid(x_nu, y_nu)

z_nu = np.sin(xx) + np.sin(yy)



In [ ]:

    
f = scipy.interpolate.interp2d(x_nu, y_nu, z_nu, kind='linear',
                               bounds_error=True)     # Let 'f' raise an exception when the requested point is outside the range defined by x and y



In [ ]:

    
# Build data

x_nu_hd = np.arange(-1*np.pi, 1*np.pi-np.pi/4, np.pi/32)
y_nu_hd = np.arange(-1*np.pi, 1*np.pi-np.pi/4, np.pi/32)

xx_nu_hd,yy_nu_hd = np.meshgrid(x_nu_hd, y_nu_hd)

z_nu_hd = np.zeros(xx_nu_hd.shape)

for xi in range(z_nu_hd.shape[0]):
    for yi in range(z_nu_hd.shape[1]):
        z_nu_hd[xi, yi] = f(x_nu_hd[xi], y_nu_hd[yi])



In [ ]:

    
# Plot data

fig = plt.figure(figsize=(12, 8))
ax = axes3d.Axes3D(fig)

surf = ax.plot_surface(xx_nu_hd,
                       yy_nu_hd,
                       z_nu_hd,
                       cmap='gnuplot2', # 'jet' # 'gnuplot2'
                       rstride=1,
                       cstride=1,
                       shade=False)

surf = ax.plot_surface(xx_hd,
                       yy_hd,
                       z_hd,
                       cmap='gnuplot2', # 'jet' # 'gnuplot2'
                       rstride=1,
                       cstride=1,
                       alpha=0.5,
                       shade=False)

plt.show();

Cubic splines



In [ ]:

    
f = scipy.interpolate.interp2d(x, y, z, kind='cubic',
                               bounds_error=True)     # Let 'f' raise an exception when the requested point is outside the range defined by x and y



In [ ]:

    
# Build data

x_hd = np.arange(-1*np.pi, 1*np.pi-np.pi/4, np.pi/32)
y_hd = np.arange(-1*np.pi, 1*np.pi-np.pi/4, np.pi/32)

xx_hd,yy_hd = np.meshgrid(x_hd, y_hd)

z_hd = np.zeros(xx_hd.shape)

for xi in range(z_hd.shape[0]):
    for yi in range(z_hd.shape[1]):
        z_hd[xi, yi] = f(x_hd[xi], y_hd[yi])



In [ ]:

    
# Plot data

fig = plt.figure(figsize=(12, 8))
ax = axes3d.Axes3D(fig)

surf = ax.plot_surface(xx_hd,
                       yy_hd,
                       z_hd,
                       cmap='gnuplot2', # 'jet' # 'gnuplot2'
                       rstride=1,
                       cstride=1,
                       shade=False)

plt.show();

Interpolate unstructured D-dimensional data

Documentation: https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.griddata.html

Scipy official documentation example

Example taken from https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.griddata.html



In [ ]:

    
def func(x, y):
    return x*(1-x)*np.cos(4*np.pi*x) * np.sin(4*np.pi*y**2)**2

grid_x, grid_y = np.mgrid[0:1:100j, 0:1:200j]

points = np.random.rand(1000, 2)
values = func(points[:,0], points[:,1])

grid_z0 = scipy.interpolate.griddata(points, values, (grid_x, grid_y), method='nearest')
grid_z1 = scipy.interpolate.griddata(points, values, (grid_x, grid_y), method='linear')
grid_z2 = scipy.interpolate.griddata(points, values, (grid_x, grid_y), method='cubic')

plt.subplot(221)
plt.imshow(func(grid_x, grid_y).T, extent=(0,1,0,1), origin='lower')
plt.plot(points[:,0], points[:,1], 'k.', ms=1)
plt.title('Original')
plt.subplot(222)
plt.imshow(grid_z0.T, extent=(0,1,0,1), origin='lower')
plt.title('Nearest')
plt.subplot(223)
plt.imshow(grid_z1.T, extent=(0,1,0,1), origin='lower')
plt.title('Linear')
plt.subplot(224)
plt.imshow(grid_z2.T, extent=(0,1,0,1), origin='lower')
plt.title('Cubic')
plt.gcf().set_size_inches(6, 6)
plt.tight_layout()
plt.show()

Callable version



In [ ]:

    
class InterpoledGridData:
    
    def __init__(self, x, y, z, interpolation_method='linear', fill_value=float('nan')):
        self.x = x
        self.y = y
        self.z = z
        self.fill_value = fill_value
        self.interpolation_method = interpolation_method
    
    def __call__(self, x1_mesh, x2_mesh):
        z = scipy.interpolate.griddata(points = (self.x, self.y),
                                       values = self.z,
                                       xi = (x1_mesh, x2_mesh),
                                       fill_value=self.fill_value,
                                       method = self.interpolation_method)
        
        if z.ndim == 0:
            z = float(z)
            
        return z



In [ ]:

    
x = np.random.rand(1000)
y = np.random.rand(1000)
z = func(x, y)

f = InterpoledGridData(x, y, z, interpolation_method='cubic')



In [ ]:

    
f(0.5, 0.5)



In [ ]:

    
x_hd = np.linspace(x.min(),
                   x.max(),
                   100)
y_hd = np.linspace(y.min(),
                   y.max(),
                   100)

xx_hd, yy_hd = np.meshgrid(x_hd, y_hd)

z_hd = f(xx_hd, yy_hd)



In [ ]:

    
# Plot data

fig = plt.figure(figsize=(12, 8))
ax = axes3d.Axes3D(fig, azim=150, elev=30)

surf = ax.plot_surface(xx_hd,
                       yy_hd,
                       z_hd,
                       cmap='gnuplot2',
                       rstride=1,
                       cstride=1,
                       vmin=np.nanmin(z_hd),
                       vmax=np.nanmax(z_hd),
                       shade=False)

plt.show();



In [ ]:

    
fig, ax = plt.subplots(figsize=(12, 8))

im = ax.pcolormesh(xx_hd, yy_hd, z_hd,
                   #shading='gouraud',
                   cmap='gnuplot2')

plt.colorbar(im, ax=ax)

plt.show();