2011-02-26 16 views
10

Estoy tratando de hacer algo de interpolación con scipy. He analizado muchos ejemplos, pero no encuentro exactamente lo que quiero.Interpolación de Python/Scipy (map_coordinates)

Digamos que tengo algunos datos donde la variable de fila y columna puede variar de 0 a 1. Los cambios delta entre cada fila y columna no siempre son los mismos (ver a continuación).

 | 0.00 0.25 0.80 1.00 
------|---------------------------- 
0.00 | 1.40 6.50 1.50 1.80 
0.60 | 8.90 7.30 1.10 1.09 
1.00 | 4.50 9.20 1.80 1.20 

Ahora quiero poder tomar un conjunto de puntos x, y y determinar los valores interpolados. Sé que puedo hacer esto con map_coordinates. Me pregunto si hay alguna manera fácil/inteligente de hacer un valor x, y para el índice apropiado en la matriz de datos.

Por ejemplo, si ingreso x, y = 0.60, 0.25, entonces debería recuperar el índice correcto para interpolar. En este caso, eso sería 1.0, 1.0 dado que 0.60, 0.25 se correlacionaría exactamente con la segunda fila y la segunda columna. x = 0.3 se correlacionaría con 0.5 ya que está a la mitad entre 0.00 y 0.60.

Sé cómo obtener el resultado que quiero, pero estoy seguro de que hay uno o dos líneas muy rápidas/claras (o una función que ya existe) que pueden hacer esto para hacer que mi código sea más claro. Básicamente, necesita interpolar por partes entre algunos arreglos.

Aquí es un ejemplo (basado en gran medida en el código de Scipy interpolation on a numpy array) - Pongo TODO en esta nueva función iría:

from scipy.ndimage import map_coordinates 
from numpy import arange 
import numpy as np 
#   0.000, 0.175, 0.817, 1.000 
z = array([ [ 3.6, 6.5, 9.1, 11.5], # 0.0000 
      [ 3.9, -7.3, 10.0, 13.1], # 0.2620 
      [ 1.9, 8.3, -15.0, -12.1], # 0.6121 
      [-4.5, 9.2, 12.2, 14.8] ]) # 1.0000 
ny, nx = z.shape 
xmin, xmax = 0., 1. 
ymin, ymax = 0., 1. 

xrange = array([0.000, 0.175, 0.817, 1.000 ]) 
yrange = array([0.0000, 0.2620, 0.6121, 1.0000]) 

# Points we want to interpolate at 
x1, y1 = 0.20, 0.45 
x2, y2 = 0.30, 0.85 
x3, y3 = 0.95, 1.00 

# 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 

yi[yi > ymax] = ymax 
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) 
yi = (ny - 1) * (yi - ymin)/(ymax - ymin) 
# TODO: Instead, xi and yi need to be mapped as described. This can only work with 
# even spacing...something like: 
#xi = SomeInterpFunction(xi, xrange) 
#yi = SomeInterpFunction(yi, yrange) 

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

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

Respuesta

15

Creo desea una bivariate spline on a rectangular structured mesh:

import numpy 
from scipy import interpolate 
x = numpy.array([0.0, 0.60, 1.0]) 
y = numpy.array([0.0, 0.25, 0.80, 1.0]) 
z = numpy.array([ 
    [ 1.4 , 6.5 , 1.5 , 1.8 ], 
    [ 8.9 , 7.3 , 1.1 , 1.09], 
    [ 4.5 , 9.2 , 1.8 , 1.2 ]]) 
# you have to set kx and ky small for this small example dataset 
# 3 is more usual and is the default 
# s=0 will ensure this interpolates. s>0 will smooth the data 
# you can also specify a bounding box outside the data limits 
# if you want to extrapolate 
sp = interpolate.RectBivariateSpline(x, y, z, kx=2, ky=2, s=0) 

sp([0.60], [0.25]) # array([[ 7.3]]) 
sp([0.25], [0.60]) # array([[ 2.66427408]]) 
+0

perfecto. Esto es exactamente lo que quiero, ¡gracias! –

+0

Publiqué una pregunta de seguimiento a esto para que también puedas ayudarme. Está aquí: . Gracias de nuevo. –