2010-06-16 37 views
9

tengo una tabla de búsqueda que se define de la siguiente manera:Scipy interpolación en una matriz de 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] ]) 
  • Los elementos de fila de encabezado son (hh) < 1,2,3,4,5+
  • La columna de cabecera (inc) elementos son < 10000, 20000, 20001+

El usuario introducirá un ejemplo de valor (1,3, 25 000), (0,2, 50 000), así sucesivamente. scipy.interpolate() debe interpolar para determinar el valor correcto.

Actualmente, la única forma en que puedo hacer esto es con un grupo de if/elifs como se ejemplifica a continuación. Estoy bastante seguro de que hay una forma mejor y más eficiente de hacerlo

Esto es lo que tengo hasta ahora:

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

Gracias por las aclaraciones! He actualizado mi respuesta a continuación. Creo que hace exactamente lo que quieres, ahora. –

Respuesta

8

Editar: Actualización cosas que reflejan sus aclaraciones anteriormente. Tu pregunta es mucho más clara ahora, ¡gracias!

Básicamente, solo quiere interpolar una matriz 2D en un punto arbitrario.

scipy.ndimage.map_coordinates es lo que quiere ....

Como entiendo su pregunta, tiene una matriz 2D de los valores "z" que va desde algunos xmin a xmax, ymin y ymax que en cada dirección.

Cualquier cosa fuera de esas coordenadas de límite desea devolver valores desde los bordes de la matriz.

map_coordinates tiene varias opciones para manejar puntos fuera de los límites de la grilla, pero ninguno de ellos hace exactamente lo que usted desea. En su lugar, podemos establecer cualquier cosa fuera de los límites para ubicarnos en el borde, y usar map_coordinates como de costumbre.

lo tanto, para utilizar map_coordinates, que necesita para convertir su coodinates reales:

 | <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 el índice Coordenadas:

 | 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 

Nota que sus límites se comportan de manera diferente en cada dirección ... En el x-direction, las cosas se comportan sin problemas, pero en la dirección y, estás pidiendo un break "duro", donde y = 20000 -> 3.9, pero y = 20000.000001 -> 4.5.

A modo de ejemplo:

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 

Esto produce:

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

Esperemos que ayuda ...

+0

que funcionó como un encanto, muchas gracias :) – dassouki

Cuestiones relacionadas