SciPy interpolazione su un array NumPy

Question

Ho una tabella di ricerca che è definita come segue:

 | <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] ]) 

• Gli elementi riga di intestazione sono (hh) < 1,2,3,4,5+
• La colonna intestazione (inc) elementi sono < 10000, 20000, 20001+

L'utente immetterà un esempio valore (1,3, 25,000), (0,2, 50000), ecc. scipy.interpolate() dovrebbe interpolare per determinare il valore corretto.

Attualmente, l'unico modo per farlo è con un gruppo di if/elifs come illustrato di seguito. Sono abbastanza sicuro che ci sia un modo migliore, più efficiente di fare questo

Ecco quello che ho finora:

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

Modifica: Aggiornato cose da riflettere le vostre precisazioni di cui sopra. La tua domanda è molto più chiara ora, grazie!

In pratica, si desidera semplicemente interpolare un array 2D in un punto arbitrario.

scipy.ndimage.map_coordinates è quello che vuoi ....

Come ho capito la tua domanda, si dispone di un array 2D di valori "Z" che varia da qualche xmin a Xmax e Ymin a yMax in ogni direzione.

Qualsiasi cosa al di fuori di queste coordinate di delimitazione, si desidera restituire i valori dai bordi della matrice.

map_coordinates ha diverse opzioni per gestire i punti al di fuori dei confini della griglia, ma nessuno di loro fa esattamente quello che vuoi. Invece, possiamo semplicemente impostare qualsiasi cosa al di fuori dei limiti per giacere sul bordo, e usare map_coordinates come al solito.

Quindi, per usare map_coordinates, è necessario trasformare le vostre coodinate effettivi:

 | <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 

Into indice coordinate:

 | 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 

Si noti che i vostri confini si comportano in modo diverso in ogni direzione ... Nel x-direction, le cose si comportano senza intoppi, ma nella direzione y, stai chiedendo una pausa "dura", dove y = 20000 -> 3.9, ma y = 20000.000001 -> 4.5.

Per fare un esempio:

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 

Questo produce:

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

Speriamo che aiuta ...