Tutorials

Interpolation with SciPy

Interpolation estimates values between known data points. It's different from fitting a model: you're not assuming a particular functional form — you just want a smooth curve that passes exactly through your measurements. You'd use it to resample a signal at a finer resolution, reconstruct a continuous surface from a coarse grid of measurements, or fill in time-series data at regular intervals. SciPy's `interpolate` module offers a range of methods from simple linear interpolation to smooth cubic splines and 2D gridded surfaces.

### Linear vs Cubic Interpolation

`interp1d` builds a callable interpolant from known (x, y) pairs. The `kind` parameter selects the method: `"linear"` connects dots with straight lines, while `"cubic"` fits smooth cubic polynomials between each pair of points.

from scipy.interpolate import interp1d
import numpy as np
import matplotlib.pyplot as plt

x_known = np.array([0, 1, 2, 3, 4, 5], dtype=float)
y_known = np.array([0, 0.8, 0.9, 0.1, -0.8, -1.0])

x_new = np.linspace(0, 5, 300)
f_linear = interp1d(x_known, y_known, kind="linear")
f_cubic  = interp1d(x_known, y_known, kind="cubic")

plt.figure(figsize=(8, 4))
plt.plot(x_known, y_known, "o", color="black", markersize=9, zorder=5, label="Known points")
plt.plot(x_new, f_linear(x_new), linestyle="--", color="orange", linewidth=2, label="Linear")
plt.plot(x_new, f_cubic(x_new), color="crimson", linewidth=2, label="Cubic")
plt.title("1D Interpolation: Linear vs Cubic")
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.tight_layout()
plt.show()
- `interp1d` returns a callable — call it like a function `f(x_new)` to get interpolated values.
- Linear interpolation (orange dashed) has visible kinks at the data points; cubic produces a smooth curve but can oscillate between points if the data is irregular.
- `interp1d` will raise a `ValueError` if you ask it to extrapolate outside the original x range. Pass `fill_value="extrapolate"` to allow it, or `bounds_error=False` to get NaN instead.

### CubicSpline for Smooth Curves

`CubicSpline` is the preferred tool when you need smooth, accurate 1D interpolation. Unlike `interp1d`'s `"cubic"` mode, it also gives you access to the first and second derivatives of the interpolant — useful for computing slopes and curvatures.

from scipy.interpolate import CubicSpline
import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(0, 2 * np.pi, 8)
y = np.sin(x)

cs = CubicSpline(x, y)
x_dense = np.linspace(0, 2 * np.pi, 400)

plt.figure(figsize=(9, 4))
plt.plot(x, y, "o", color="black", markersize=8, zorder=5, label="Sample points")
plt.plot(x_dense, np.sin(x_dense), color="steelblue", linewidth=2, label="True sin(x)")
plt.plot(x_dense, cs(x_dense), color="crimson", linewidth=2, linestyle="--", label="CubicSpline")
plt.plot(x_dense, cs(x_dense, 1), color="orange", linewidth=1.5, linestyle=":", label="Derivative (≈ cos)")
plt.title("CubicSpline: Interpolation and Derivative")
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.tight_layout()
plt.show()
- `cs(x_dense)` evaluates the spline; `cs(x_dense, 1)` evaluates its first derivative — the `1` is the derivative order.
- With only 8 sample points, the spline (red dashed) closely tracks the true sine (blue), demonstrating the accuracy of cubic interpolation for smooth functions.
- The derivative (orange dotted) approximates cos(x) — useful for computing local slopes without finite differencing.

### Monotone Interpolation with PchipInterpolator

`CubicSpline` can overshoot between data points if the data has rapid changes. `PchipInterpolator` (Piecewise Cubic Hermite Interpolating Polynomial) avoids overshooting by preserving monotonicity — if the data increases, the interpolant never dips below the data.

from scipy.interpolate import CubicSpline, PchipInterpolator
import numpy as np
import matplotlib.pyplot as plt

x = np.array([0, 1, 2, 3, 4, 5, 6], dtype=float)
y = np.array([0, 0, 0, 1, 1, 1, 1], dtype=float)  # step-like data

x_new = np.linspace(0, 6, 300)
cs  = CubicSpline(x, y)
pch = PchipInterpolator(x, y)

plt.figure(figsize=(8, 4))
plt.plot(x, y, "o", color="black", markersize=9, zorder=5, label="Data")
plt.plot(x_new, cs(x_new),  color="orange", linewidth=2, linestyle="--", label="CubicSpline (overshoots)")
plt.plot(x_new, pch(x_new), color="crimson", linewidth=2, label="PCHIP (monotone)")
plt.title("CubicSpline vs PCHIP on Step Data")
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.tight_layout()
plt.show()
- `CubicSpline` (orange) dips below 0 and overshoots above 1 around the step — physically wrong if y represents a quantity that can't be negative.
- `PchipInterpolator` (red) stays within the data range, making it the right choice for monotone or bounded data like probabilities or concentrations.
- Use `CubicSpline` for smooth natural data; use `PchipInterpolator` when the data has abrupt changes and you need the interpolant to stay well-behaved.

### 2D Gridded Interpolation

`RegularGridInterpolator` interpolates over a regular 2D (or higher-dimensional) grid — for example, a coarse elevation map that you want to sample at higher resolution.

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

# Coarse 6×6 grid
x = np.linspace(0, 2 * np.pi, 6)
y = np.linspace(0, 2 * np.pi, 6)
X, Y = np.meshgrid(x, y, indexing="ij")
Z = np.sin(X) * np.cos(Y)

interp = RegularGridInterpolator((x, y), Z, method="cubic")

# Evaluate on a fine 100×100 grid
x_fine = np.linspace(0, 2 * np.pi, 100)
y_fine = np.linspace(0, 2 * np.pi, 100)
Xf, Yf = np.meshgrid(x_fine, y_fine, indexing="ij")
pts = np.stack([Xf.ravel(), Yf.ravel()], axis=-1)
Z_interp = interp(pts).reshape(100, 100)

fig, axes = plt.subplots(1, 2, figsize=(10, 4))
kw = dict(extent=[0, 2*np.pi, 0, 2*np.pi], origin="lower", cmap="RdBu_r", vmin=-1, vmax=1)
axes[0].imshow(Z.T, **kw)
axes[0].set_title("Original (6×6 grid)")
im = axes[1].imshow(Z_interp.T, **kw)
axes[1].set_title("Interpolated (100×100)")
for ax in axes:
    ax.set_xlabel("x")
    ax.set_ylabel("y")
plt.colorbar(im, ax=axes[1])
plt.tight_layout()
plt.show()
- `RegularGridInterpolator((x, y), Z)` takes the 1D grid coordinates and the 2D array of values. The axes of `Z` must align with the order of coordinates.
- `np.stack([Xf.ravel(), Yf.ravel()], axis=-1)` packages the query points as an (N, 2) array — the format `RegularGridInterpolator` expects.
- The coarse 6×6 surface (left) becomes a smooth 100×100 surface (right) while preserving the values at the original grid points exactly.

Interpolation and smoothing are complementary: interpolation passes through the data exactly, while smoothing deliberately deviates to remove noise. For the latter, see [signal filtering and smoothing](/tutorials/signal-filtering-and-smoothing-with-scipy).