multinterp

A Unified Interface for Multivariate Interpolation in the Scientific Python Ecosystem

from __future__ import annotations

from time import time

import matplotlib.pyplot as plt
import numpy as np
from multinterp.curvilinear import Curvilinear2DInterp, Warped2DInterp

Suppose we have a collection of values for an unknown function along with their respective coordinate points. For illustration, assume the values come from the following function:

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

The points are randomly scattered in the unit square and therefore have no regular structure. This is achieved by randomly shifting a well structured grid at every point.

rng = np.random.default_rng(0)
warp_factor = 0.01
x_list = np.linspace(0, 1, 20)
y_list = np.linspace(0, 1, 20)
x_temp, y_temp = np.meshgrid(x_list, y_list, indexing="ij")
rand_x = x_temp + warp_factor * (rng.random((x_list.size, y_list.size)) - 0.5)
rand_y = y_temp + warp_factor * (rng.random((x_list.size, y_list.size)) - 0.5)
values = function_1(rand_x, rand_y)

Now suppose we would like to interpolate this function on a rectilinear grid, which is known as “regridding”.

grid_x, grid_y = np.meshgrid(
    np.linspace(0, 1, 100),
    np.linspace(0, 1, 100),
    indexing="ij",
)

To do this, we use multinterp’s Warped2DInterp and Curvilinear2DInterp classes. The class takes the following arguments:

  • values: an ND-array of values for the function at the points
  • grids: a list of ND-arrays of coordinates for the points
  • backend: the backend to use for interpolation, currently only scipy and numba are supported for Warped2DInterp and only scipy is supported for Curvilinear2DInterp
warped_interp = Warped2DInterp(values, (rand_x, rand_y), backend="numba")
warped_interp.warmup()

Once we create the interpolator objects, we can evaluate the functions on the query grids and compare their time performance.

start = time()
warped_grid = warped_interp(grid_x, grid_y)
print(f"Warped interpolation took {time() - start:.5f} seconds")
Warped interpolation took 0.00240 seconds
curvilinear_interp = Curvilinear2DInterp(values, (rand_x, rand_y))
start = time()
curvilinear_grid = curvilinear_interp(grid_x, grid_y)
print(f"Curvilinear interpolation took {time() - start:.5f} seconds")
Curvilinear interpolation took 0.00802 seconds

Now we can compare the results of the interpolation with the original function. Below we plot the original function and the sample points that are known. Notice that the points are almost rectilinear, but have been randomly shifted to create a more challenging interpolation problem.

plt.imshow(function_1(grid_x, grid_y).T, extent=(0, 1, 0, 1), origin="lower")
plt.plot(rand_x.flat, rand_y.flat, "ok", ms=2, label="input points")
plt.title("Original")
plt.legend(loc="lower right")
<Figure size 640x480 with 1 Axes>

Then, we can look at the result for each method of interpolation and compare it to the original function.

fig, axs = plt.subplots(1, 3, figsize=(9, 6))
titles = ["Original", "WarpedInterp", "CurvilinearInterp"]
grids = [function_1(grid_x, grid_y), warped_grid, curvilinear_grid]

for ax, title, grid in zip(axs.flat, titles, grids):
    im = ax.imshow(grid.T, extent=(0, 1, 0, 1), origin="lower")
    ax.set_title(title)

plt.tight_layout()
plt.show()
<Figure size 900x600 with 3 Axes>

In short, multinterp’s Warped2DInterp and Curvilinear2DInterp classes are useful for interpolating functions on curvilinear grids which have a quadrilateral structure but are not perfectly rectangular.