2010-06-16 6 views
9

j'ai une table de consultation qui est définie de la manière suivante:Scipy interpolation sur une matrice numpy

 | <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] ]) 
  • Les éléments de la ligne d'en-tête sont (hh) < 1,2,3,4,5+
  • La colonne d'en-tête (inc) éléments sont < 10000, 20000 20001+

l'utilisateur va entrer un exemple de valeur (1,3, 25 000), (0,2, 50 000), ainsi de suite. scipy.interpolate() doit interpoler pour déterminer la valeur correcte.

Actuellement, la seule façon de le faire est avec un groupe de if/elifs comme illustré ci-dessous. Je suis sûr qu'il ya une meilleure façon plus efficace de le faire

Voici ce que j'ai jusqu'à présent:

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

Merci pour les précisions! J'ai mis à jour ma réponse ci-dessous. Je pense qu'il fait exactement ce que tu veux, maintenant. –

Répondre

8

Modifier: les choses à jour pour refléter vos éclaircissements ci-dessus. Votre question est beaucoup plus claire maintenant, merci! Fondamentalement, vous voulez juste interpoler un tableau 2D à un point arbitraire.

scipy.ndimage.map_coordinates est ce que vous voulez ....

Si je comprends bien votre question, vous avez un tableau 2D de valeurs « z » qui va de certains xmin à xmax et ymin à ymax dans chaque direction.

N'importe quoi en dehors de ces coordonnées de délimitation vous voulez renvoyer des valeurs à partir des bords du tableau. Map_coordinates a plusieurs options pour gérer les points en dehors des limites de la grille, mais aucun d'entre eux ne fait exactement ce que vous voulez. Au lieu de cela, nous pouvons simplement définir n'importe quoi en dehors des limites pour rester sur le bord, et utiliser map_coordinates comme d'habitude.

Donc, pour utiliser map_coordinates, vous devez transformer vos coodinates réels:

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

en coordonnées index:

 | 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 

Notez que vos limites se comportent différemment dans chaque sens ... Dans la direction x, les choses se comportent bien, mais dans la direction y, vous demandez une rupture "dure", où y = 20000 -> 3.9, mais y = 20000.000001 -> 4.5.

À titre d'exemple:

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 

Cela donne:

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

Espérons que helps ...

+0

qui a fonctionné comme un charme, merci beaucoup :) – dassouki

Questions connexes