Scipy Interpolation auf einer numpy Array

Question

I eine Nachschlagetabelle aufweisen, die folgendermaßen definiert ist:

 | <1 2 3 4 5+ 
-------|---------------------------- 
<10000 | 3.6 6.5 9.1 11.5 13.8 
20000 | 3.9 7.3 10.0 13.1 15.9 
20000+ | 4.5 9.2 12.2 14.8 18.2 


TR_ua1 = np.array([ [3.6, 6.5, 9.1, 11.5, 13.8], 
        [3.9, 7.3, 10.0, 13.1, 15.9], 
        [4.5, 9.2, 12.2, 14.8, 18.2] ]) 

• Die Kopfzeile Elemente (hh) < 1,2,3,4,5+
• Die Header-Säule (inc) Elemente sind < 10000, 20000, 20001+

die Benutzer einen Wert eingeben Beispiel (1.3, 25.000), (0,2, 50,000), so weiter. scipy.interpolate() sollte interpolieren, um den korrekten Wert zu bestimmen.

Derzeit ist die einzige Möglichkeit, die ich tun kann, mit einer Reihe von if/elifs wie unten beispielhaft dargestellt. Ich bin mir ziemlich sicher, dass es eine bessere, effizientere Art und Weise, dies zu tun

Hier ist, was ich bisher habe:

import numpy as np 
from scipy import interpolate 

if (ua == 1): 
    if (inc <= low_inc): # low_inc = 10,000 
     if (hh <= 1): 
     return TR_ua1[0][0] 
     elif (hh >= 1 & hh < 2): 
     return interpolate((1, 2), (TR_ua1[0][1], TR_ua1[0][2])) 

Answer 1

Edit: Aktualisiert Dinge Ihre Erläuterungen zu reflektieren oben. Ihre Frage ist jetzt viel klarer, danke!

Grundsätzlich Sie wollen nur ein 2D-Array an einem beliebigen Punkt zu interpolieren.

scipy.ndimage.map_coordinates ist das, was Sie wollen ....

Wie ich verstehe Ihre Frage, können Sie ein 2D-Array von „z“ Werte haben, die von einigen xmin reicht bis xmax und ymin in jeder Richtung ymax.

Alles, was außerhalb dieser Begrenzungskoordinaten Sie Werte von den Rändern des Arrays zurückgeben möchten.

map_coordinates hat mehrere Optionen Punkte außerhalb der Grenzen des Gitters zu handhaben, aber keiner von ihnen genau das tun, was Sie wollen. Stattdessen können wir nur etwas außerhalb der Grenzen am Rande liegen, und map_coordinates wie gewohnt verwenden.

So verwendet map_coordinates, müssen Sie Ihren tatsächlichen von Koordinaten drehen:

 | <1 2 3 4 5+ 
-------|---------------------------- 
<10000 | 3.6 6.5 9.1 11.5 13.8 
20000 | 3.9 7.3 10.0 13.1 15.9 
20000+ | 4.5 9.2 12.2 14.8 18.2 

In Index Koordinaten:

 | 0  1 2 3 4 
-------|---------------------------- 
    0 | 3.6 6.5 9.1 11.5 13.8 
    1 | 3.9 7.3 10.0 13.1 15.9 
    2 | 4.5 9.2 12.2 14.8 18.2 

Beachten Sie, dass Ihre Grenzen unterschiedlich in jede Richtung verhalten ... In der x-Richtung, verhalten sich die Dinge glatt, aber in der y-Richtung, sind Sie für eine "harte" Pause zu fragen, wo y = 20000 -> 3,9, aber y = 20000,000001 -> 4.5.

Als Beispiel:

import numpy as np 
from scipy.ndimage import map_coordinates 

#-- Setup --------------------------- 
z = np.array([ [3.6, 6.5, 9.1, 11.5, 13.8], 
       [3.9, 7.3, 10.0, 13.1, 15.9], 
       [4.5, 9.2, 12.2, 14.8, 18.2] ]) 
ny, nx = z.shape 
xmin, xmax = 1, 5 
ymin, ymax = 10000, 20000 

# Points we want to interpolate at 
x1, y1 = 1.3, 25000 
x2, y2 = 0.2, 50000 
x3, y3 = 2.5, 15000 

# To make our lives easier down the road, let's 
# turn these into arrays of x & y coords 
xi = np.array([x1, x2, x3], dtype=np.float) 
yi = np.array([y1, y2, y3], dtype=np.float) 

# Now, we'll set points outside the boundaries to lie along an edge 
xi[xi > xmax] = xmax 
xi[xi < xmin] = xmin 

# To deal with the "hard" break, we'll have to treat y differently, 
# so we're ust setting the min here... 
yi[yi < ymin] = ymin 

# We need to convert these to (float) indicies 
# (xi should range from 0 to (nx - 1), etc) 
xi = (nx - 1) * (xi - xmin)/(xmax - xmin) 

# Deal with the "hard" break in the y-direction 
yi = (ny - 2) * (yi - ymin)/(ymax - ymin) 
yi[yi > 1] = 2.0 

# Now we actually interpolate 
# map_coordinates does cubic interpolation by default, 
# use "order=1" to preform bilinear interpolation instead... 
z1, z2, z3 = map_coordinates(z, [yi, xi]) 

# Display the results 
for X, Y, Z in zip((x1, x2, x3), (y1, y2, y3), (z1, z2, z3)): 
    print X, ',', Y, '-->', Z 

Dies ergibt:

1.3 , 25000 --> 5.1807375 
0.2 , 50000 --> 4.5 
2.5 , 15000 --> 8.12252371652 

Hoffentlich ...