Multivariate Dateninterpolation auf einem regulären Gitter (RegularGridInterpolator)

Angenommen, Sie haben N-dimensionale Daten auf einem regulären Gitter und möchten diese interpolieren. In diesem Fall kann RegularGridInterpolator nützlich sein. Mehrere Interpolationsstrategien werden unterstützt: Nearest-Neighbor, linear und Tensorprodukt-Splines ungeraden Grades.

Genau genommen behandelt diese Klasse effizient Daten, die auf **rechtwinkligen** Gittern gegeben sind: hyperkubische Gitter mit möglicherweise ungleicher Abständen zwischen den Punkten. Die Anzahl der Punkte pro Dimension kann für verschiedene Dimensionen unterschiedlich sein.

Das folgende Beispiel demonstriert seine Verwendung und vergleicht die Interpolationsergebnisse mit jeder Methode.

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

Angenommen, wir möchten diese 2D-Funktion interpolieren.

>>> def F(u, v):
...     return u * np.cos(u * v) + v * np.sin(u * v)

Angenommen, wir kennen nur einige Daten auf einem regulären Gitter.

>>> fit_points = [np.linspace(0, 3, 8), np.linspace(0, 3, 11)]
>>> values = F(*np.meshgrid(*fit_points, indexing='ij'))

Erstellen von Testpunkten und wahren Werten für Auswertungen.

>>> ut, vt = np.meshgrid(np.linspace(0, 3, 80), np.linspace(0, 3, 80), indexing='ij')
>>> true_values = F(ut, vt)
>>> test_points = np.array([ut.ravel(), vt.ravel()]).T

Wir können den Interpolator erstellen und Testpunkte mit jeder Methode interpolieren.

>>> interp = RegularGridInterpolator(fit_points, values)
>>> fig, axes = plt.subplots(2, 3, figsize=(10, 6))
>>> axes = axes.ravel()
>>> fig_index = 0
>>> for method in ['linear', 'nearest', 'slinear', 'cubic', 'quintic']:
...     im = interp(test_points, method=method).reshape(80, 80)
...     axes[fig_index].imshow(im)
...     axes[fig_index].set_title(method)
...     axes[fig_index].axis("off")
...     fig_index += 1
>>> axes[fig_index].imshow(true_values)
>>> axes[fig_index].set_title("True values")
>>> fig.tight_layout()
>>> fig.show()

Wie erwartet, sind die Spline-Interpolationen höheren Grades den wahren Werten am nächsten, obwohl ihre Berechnung teurer ist als bei linear oder nearest. Die slinear-Interpolation stimmt auch mit der linear-Interpolation überein.

Wenn Ihre Daten dazu führen, dass Spline-Methoden Ringing-Effekte erzeugen, können Sie method="pchip" in Betracht ziehen, das das Tensorprodukt von PCHIP-Interpolatoren verwendet, einem PchipInterpolator pro Dimension.

Wenn Sie eine funktionale Schnittstelle gegenüber der expliziten Erstellung einer Klasseninstanz bevorzugen, bietet die Komfortfunktion interpn die äquivalente Funktionalität.

Insbesondere ergeben diese beiden Formen identische Ergebnisse

>>> from scipy.interpolate import interpn
>>> rgi = RegularGridInterpolator(fit_points, values)
>>> result_rgi = rgi(test_points)

und

>>> result_interpn = interpn(fit_points, values, test_points)
>>> np.allclose(result_rgi, result_interpn, atol=1e-15)
True

Für Daten, die auf einen (N-1)-dimensionalen Unterraum des N-dimensionalen Raums beschränkt sind, d.h. wenn eine der Gitterachsen die Länge 1 hat, wird die Extrapolation entlang dieser Achse durch den Schlüsselwortparameter fill_value gesteuert

>>> x = np.array([0, 5, 10])
>>> y = np.array([0])
>>> data = np.array([[0], [5], [10]])
>>> rgi = RegularGridInterpolator((x, y), data,
...                               bounds_error=False, fill_value=None)
>>> rgi([(2, 0), (2, 1), (2, -1)])   # extrapolates the value on the axis
array([2., 2., 2.])
>>> rgi.fill_value = -101
>>> rgi([(2, 0), (2, 1), (2, -1)])
array([2., -101., -101.])

Hinweis

Wenn die Eingabedaten so sind, dass die Eingabedimensionen inkommensurable Einheiten haben und sich um viele Größenordnungen unterscheiden, kann der Interpolator numerische Artefakte aufweisen. Ziehen Sie in Erwägung, die Daten vor der Interpolation zu skalieren.

Batch-Dimensionen von values

Angenommen, Sie haben eine Vektorfunktion \(f(x) = y\), wobei \(x\) und \(y\) Vektoren sind, möglicherweise unterschiedlicher Länge, und Sie möchten die Funktion auf einem Gitter von \(x\)-Werten abtasten. Ein Weg, dies anzugehen, ist die Nutzung der Tatsache, dass RegularGridInterpolator values mit **nachfolgenden** Dimensionen zulässt.

Im Einklang damit, wie 1D-Interpolatoren mehrdimensionale Arrays interpretieren, wird angenommen, dass die ersten \(N\) Dimensionen der values-Arrays Daten-Dimensionen sind (d.h. sie entsprechen den durch das grid-Argument definierten Punkten), und die **nachfolgenden** Dimensionen sind Batch-Achsen. Beachten Sie, dass dies von den üblichen NumPy-Broadcasting-Konventionen abweicht, wo das Broadcasting entlang der **führenden** Dimensionen erfolgt.

Zur Veranschaulichung

>>> n = 5   # the number of batch components

>>> # make a 3D grid
>>> x1 = np.linspace(-np.pi, np.pi, 10)
>>> x2 = np.linspace(0.0, np.pi, 15)
>>> x3 = np.linspace(0.0, np.pi/2, 20)
>>> points = (x1, x2, x3)
>>>
>>> # define a function and sample it on the grid
>>> def f(x1, x2, x3, n):
...     lst = [np.sin(np.pi*x1/2) * np.exp(x2/2) + x3 + i for i in range(n)]
...     return np.asarray(lst)
>>>
>>> X1, X2, X3 = np.meshgrid(x1, x2, x3, indexing="ij")
>>> values = f(X1, X2, X3, n)
>>> values.shape
(5, 10, 15, 20)
>>>
>>> # prepare the data and construct the interpolator
>>> values = np.moveaxis(values, 0, -1)
>>> values.shape
(10, 15, 20, 5)     # the batch dimension is 5
>>> rgi = RegularGridInterpolator(points, values)
>>>
>>> # Coordinates to compute the interpolation at
>>> x = np.asarray([0.2, np.pi/2.1, np.pi/4.1])
>>>
# evaluate
>>> rgi(x).shape
(1, 5)

In diesem Beispiel haben wir einen Batch von \(n=5\) Funktionen auf einem dreidimensionalen Gitter ausgewertet. Im Allgemeinen sind mehrere Batch-Dimensionen zulässig, und die Form des Ergebnisses ergibt sich durch Anhängen der Batch-Form (in diesem Beispiel (5,)) an die Form der Eingabe x (in diesem Beispiel (1,)`).

Gleichmäßig beabstandete Daten

Wenn Sie mit Daten auf kartesischen Gittern mit ganzzahligen Koordinaten arbeiten, z. B. beim Neuabtasten von Bilddaten, sind diese Routinen möglicherweise nicht die optimale Wahl. Betrachten Sie stattdessen scipy.ndimage.map_coordinates.

Für Gleitkommadaten auf Gittern mit gleichem Abstand kann map_coordinates leicht in eine RegularGridInterpolator-ähnliche Funktion verpackt werden. Das Folgende ist ein minimales Beispiel, das aus dem 'regulargrid'-Paket von Johannes Buchner stammt

class CartesianGridInterpolator:
    def __init__(self, points, values, method='linear'):
        self.limits = np.array([[min(x), max(x)] for x in points])
        self.values = np.asarray(values, dtype=float)
        self.order = {'linear': 1, 'cubic': 3, 'quintic': 5}[method]

    def __call__(self, xi):
        """
        `xi` here is an array-like (an array or a list) of points.

        Each "point" is an ndim-dimensional array_like, representing
        the coordinates of a point in ndim-dimensional space.
        """
        # transpose the xi array into the ``map_coordinates`` convention
        # which takes coordinates of a point along columns of a 2D array.
        xi = np.asarray(xi).T

        # convert from data coordinates to pixel coordinates
        ns = self.values.shape
        coords = [(n-1)*(val - lo) / (hi - lo)
                  for val, n, (lo, hi) in zip(xi, ns, self.limits)]

        # interpolate
        return map_coordinates(self.values, coords,
                               order=self.order,
                               cval=np.nan)  # fill_value

Dieser Wrapper kann als (fast) direkter Ersatz für die RegularGridInterpolator verwendet werden

>>> x, y = np.arange(5), np.arange(6)
>>> xx, yy = np.meshgrid(x, y, indexing='ij')
>>> values = xx**3 + yy**3
>>> rgi = RegularGridInterpolator((x, y), values, method='linear')
>>> rgi([[1.5, 1.5], [3.5, 2.6]])
array([ 9. , 64.9])
>>> cgi = CartesianGridInterpolator((x, y), values, method='linear')
>>> cgi([[1.5, 1.5], [3.5, 2.6]])
array([ 9. , 64.9])

Beachten Sie, dass das obige Beispiel die Randbedingungen von map_coordinates verwendet. Daher können die Ergebnisse der cubic- und quintic-Interpolationen von denen der RegularGridInterpolator abweichen. Weitere Details zu Randbedingungen und anderen zusätzlichen Argumenten finden Sie in der Dokumentation zu scipy.ndimage.map_coordinates. Schließlich stellen wir fest, dass dieses vereinfachte Beispiel davon ausgeht, dass die Eingabedaten in aufsteigender Reihenfolge gegeben sind.