Python-4D lineare interpolation auf einem rechteckigen Gitter

In der gleichen ticket, das Sie verlinkt haben, es ist ein Beispiel für die Umsetzung dessen, was Sie nennen tensorprodukt-interpolation, zeigt den richtigen Weg zum nest rekursive Aufrufe interp1d. Dies entspricht quadrilineare interpolation wenn Sie die Standardeinstellung kind='linear' parameter für Ihre interp1d's.

Zwar mag dies gut genug sein, diese ist nicht linear interpoliert, und es wird höher sein, um Begriffe in der interpolation-Funktion, wie dieses Bild aus den wikipedia-Eintrag über Bilineare interpolation zeigt:

Diese kann sehr gut sein, gut genug für das, was Sie nach sind, aber es gibt Anwendungen, bei denen eine trianguliert, wirklich stückweise lineare, interpoaltion bevorzugt. Wenn Sie wirklich brauchen, gibt es eine einfache Art und Weise zu arbeiten, um die Langsamkeit von qhull.

Einmal LinearNDInterpolator eingerichtet, es gibt zwei Schritte, um zu kommen mit einer interpolierten Wert für einen gegebenen Punkt:

1. herauszufinden, die im inneren Dreieck (4D hypertetrahedron in Ihrem Fall) der Punkt ist, und
2. interpolieren mit der barycentric coordinates des Punktes relativ zu den vertices als GEWICHTE.

Werden Sie wahrscheinlich nicht wollen, um Durcheinander mit barycentric coordinates, also besser verlassen, dass LinearNDInterpolator. Aber Sie weiß einige Dinge über die triangulation. Meist, weil Sie ein regelmäßiges raster, in jeder hypercube der triangulation ist die gleiche. So zu interpolieren, um einen einzelnen Wert, könnten Sie zuerst bestimmen, in dem Teilcube, der Ihre Stelle ist, bauen Sie ein LinearNDInterpolator mit der 16 Eckpunkte des Würfels, und es verwenden, um zu interpolieren Wert:

from itertools import product

def interpolator(coords, data, point) :
    dims = len(point)
    indices = []
    sub_coords = []
    for j in xrange(dims) :
        idx = np.digitize([point[j]], coords[j])[0]
        indices += [[idx - 1, idx]]
        sub_coords += [coords[j][indices[-1]]]
    indices = np.array([j for j in product(*indices)])
    sub_coords = np.array([j for j in product(*sub_coords)])
    sub_data = data[list(np.swapaxes(indices, 0, 1))]
    li = LinearNDInterpolator(sub_coords, sub_data)
    return li([point])[0]

>>> point = np.array([12.3,-4.2, 500.5, 2.5])
>>> interpolator((lats, lons, alts, time), data, point)
0.386082399091

Kann das nicht funktionieren auf vektorisierte Daten, denn das würde erfordern das speichern einer LinearNDInterpolator für jeden möglich Teilcube, und obwohl es wäre wahrscheinlich schneller als triangulation wird die ganze Sache, wäre es immer noch sehr langsam.

Für sehr ähnliche Dinge, die ich verwenden Wissenschaftliche.Funktionen.Interpolation.InterpolatingFunction.

    import numpy as np
    from Scientific.Functions.Interpolation import InterpolatingFunction

    lats = np.arange(-90,90.5,0.5)
    lons = np.arange(-180,180,0.5)
    alts = np.arange(1,1000,21.717)
    time = np.arange(8)
    data = np.random.rand(len(lats)*len(lons)*len(alts)*len(time)).reshape((len(lats),len(lons),len(alts),len(time)))

    axes = (lats, lons, alts, time)
    f = InterpolatingFunction(axes, data)

Können Sie jetzt lassen Sie es an die Benutzer zu nennen, die InterpolatingFunction mit den Koordinaten:

>>> f(0,0,10,3)
0.7085675631375401

InterpolatingFunction hat nette zusätzliche features, wie die integration und aufschneiden.

Allerdings weiß ich nicht sicher, ob die interpolation linear ist. Sie würden schauen müssen, die im Modul Quelle zu finden, aus.