1D-Interpolation… von 2D-Gittern

Question

Ich weiß, dass dies ziemlich verwirrend sein kann. Bitte lassen Sie mich wissen, wenn diese Erklärung bearbeitet werden muss.

Angenommen, ich habe Eingabedaten in diesem Format:

für einen gegebenen Druck p_0-> 2x2 Gitter von Temperaturen ( T_0), die sich auf diesen Druckwert beziehen

für einen gegebenen Druck p_1-> 2x2 Gitter von Temperaturen ( T_1), die sich auf diesen Druckwert beziehen

p_0 = 0
T_0 = np.array([[1, 4], [3, 2]])

p_1 = 1
T_1 = np.array([[1, 6], [4, 4]])

p = np.array([p_0, p_1])
T = np.array([T_0, T_1])

Jetzt bekomme ich ein 2x2-Raster mit neuen Druckwerten

p_target = np.array([[0.1, 0.4], [0.3, 0.2]])

und ich möchte ein 2x2-Raster interpolierter Temperaturwerte unter Verwendung der Eingabedaten erhalten.

Ich mache das für jeden Punkt des Gitters. Ich erstelle eine Interpolationsfunktion und verwende sie dann, um den neuen interpolierten Temperaturwert für diesen Gitterpunkt zu erhalten:

from scipy.interpolate import interp1d

T_new = np.empty(p_target.shape)

for ix,iy in np.ndindex(p_target.shape):
    f = interp1d(p, T[:,ix,iy])
    T_new[ix,iy] = f(p_target[ix,iy])

T_new

array([[1. , 4.8],
       [3.3, 2.4]])

Wie man leicht erraten kann, ist dies für große Arrays ziemlich langsam, und es scheint ziemlich gegen die numpy Art und Weise zu sein, Dinge zu tun.

BEARBEITEN: Ich verwende es interp1dauch, weil es auch eine Extrapolation ermöglicht, eine Option, die ich behalten möchte.

Answer 1

Sie können die Interpolation einfach selbst berechnen. Hier gehe ich davon aus, dass Sie mehr als zwei TWerte haben und das pnicht unbedingt gleichmäßig verteilt ist. Außerdem geht der Code davon aus, dass Sie mehrere p_targetWerte haben, funktioniert aber offensichtlich nur für einen einzigen.

import numpy as np

p_0 = 0
T_0 = np.array([[1., 4.], [3., 2.]])
p_1 = 1
T_1 = np.array([[1., 6.], [4., 4.]])
p = np.array([p_0, p_1])
T = np.array([T_0, T_1])
p_target = np.array([[0.1, 0.4], [0.3, 0.2]])
# Assume you may have several of p_target values
p_target = np.expand_dims(p_target, 0)

# Find the base index for each interpolated value (assume p is sorted)
idx_0 = (np.searchsorted(p, p_target) - 1).clip(0, len(p) - 2)
# And the next index
idx_1 = idx_0 + 1
# Get p values for each interpolated value
a = p[idx_0]
b = p[idx_1]
# Compute interpolation factor
alpha = ((p_target - a) / (b - a)).clip(0, 1)
# Get interpolation values
v_0 = np.take_along_axis(T, idx_0, axis=0)
v_1 = np.take_along_axis(T, idx_1, axis=0)
# Compute interpolation
out = (1 - alpha) * v_0 + alpha * v_1
print(out)
# [[[1.  4.8]
#   [3.3 2.4]]]

BEARBEITEN: Wenn Sie eine lineare Extrapolation wünschen, schneiden Sie die alphaWerte einfach nicht ab :

alpha = ((p_target - a) / (b - a))

Answer 2

Ich habe einige Parameter für die Dimensionen hinzugefügt. Nach Ihrer Wahl n_x = n_y = n_p = 2waren die Abhängigkeiten nicht so klar.

from scipy.interpolate import interp1d, interp2d, dfitpack

n_x = 30
n_y = 40
n_p = 50
T = np.random.random((n_p, n_x, n_y)) * 100
p = np.random.random(n_p)
p[np.argmin(p)] = 0
p[np.argmax(p)] = 1
p_target = np.random.random((n_x, n_y))

T_new = np.empty(p_target.shape)

for ix, iy in np.ndindex(p_target.shape):
    f = interp1d(p, T[:, ix, iy])
    T_new[ix, iy] = f(p_target[ix, iy])

Als ein Wort zu Ihrer Modellierung. Wenn ich richtig verstanden habe, möchten Sie temperature_xy = fun_xy(pressure)eine separate Funktion für jede Koordinate in Ihrem räumlichen Raster. Eine andere Möglichkeit könnte darin bestehen, die räumlichen Komponenten in eine kombinierte Funktion aufzunehmen temperature_xy = fun(pressure, x, y). Für den zweiten Ansatz schauen Sie sich scipy.interpolate.griddata an .

Sie können den ersten Ansatz neu anordnen, damit er funktioniert interp2d(). Hierzu ist die erste Dimension der Druck x=pressureund die zweite Dimension die kombinierten räumlichen Dimensionen y=product(x, y). Damit sich dies als n_x * n_yunabhängige Interpolation der Druckwerte verhält , verwende ich nur die gleichen Dummy-Werte 0, 1, 2 ... für die räumlichen Komponenten, sowohl beim Erstellen der Interpolation als auch bei der Auswertung. Da die Auswertung der interp2d()Normalität nur für Gitterkoordinaten funktioniert, habe ich die von user6655984 bereitgestellte Methode verwendet , um die Funktion nur für einen bestimmten Satz von Punkten auszuwerten.

def evaluate_interp2d(f, x, y):
    """https://stackoverflow.com/a/47233198/7570817"""
    return dfitpack.bispeu(f.tck[0], f.tck[1], f.tck[2], f.tck[3], f.tck[4], x, y)[0]

f2 = interp2d(x=p, y=np.arange(n_x*n_y), z=T.reshape(n_p, n_x*n_y).T)

T_new2 = evaluate_interp2d(f=f2, x=p_target.ravel(), y=np.arange(n_x*n_y))
T_new2 = T_new2.reshape(n_x, n_y)

print(np.allclose(T_new, T_new2))
# True

Mit diesen Einstellungen bekomme ich eine Zeitverbesserung von fast 10x. Wenn Sie jedoch noch größere Werte verwenden, wird n_x=n_y=1000die Speichernutzung dieses benutzerdefinierten interp2d- Ansatzes zu groß und Ihr iterativer Ansatz gewinnt.

# np=50
#    nx*ny      1e2      1e4      1e5      1e6
# interp1d  0.0056s  0.3420s  3.4133s  33.390s
# interp2d  0.0004s  0.0388s  2.0954s  191.66s

Mit diesem Wissen könnten Sie ein großes 1000x1000Gitter durchlaufen und 100x100Teile nacheinander verarbeiten, dann würden Sie bei etwa 3 Sekunden statt 30 Sekunden enden.

def interpolate2d_flat(p, p_target_flat, T_flat):
    n_p, n_xy = T_flat.shape
    f2 = interp2d(x=p, y=np.arange(n_xy), z=T_flat.T)
    return evaluate_interp2d(f=f2, x=p_target_flat, y=np.arange(n_xy))


n_splits = n_x * n_y // 1000  # So each patch has size n_p*1000, can be changed 

# Flatten and split the spatial dimensions
T_flat_s = np.array_split(T.reshape(n_p, n_x*n_y), n_splits, axis=1)
p_target_flat_s = np.array_split(p_target.ravel(), n_splits, axis=0)

# Loop over the patches
T_new_flat = np.concatenate([interpolate2d_flat(p=p, p_target_flat=ptf, T_flat=Tf)
                             for (ptf, Tf) in zip(p_target_flat_s, T_flat_s)])
T_new2 = T_new_flat.reshape(n_x, n_y)