2011-02-26 104 views
10

我试图用scipy做一些插值。我已经通过了很多例子,但我并没有找到我想要的东西。Python/Scipy插值(map_coordinates)

比方说,我有一些数据,行和列变量可以从0到1变化。每一行和列之间的增量变化并不总是相同的(见下文)。

 | 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 

现在我想能够取一组x,y点并确定插值。我知道我可以用map_coordinates做到这一点。我想知道是否有任何简单/巧妙的方法来为数据数组中的适当索引创建x,y值。例如,如果我输入x,y = 0.60,0.25,那么我应该找回要插入的正确索引。在这种情况下,这将是1.0,1.0,因为0.60,0.25将完全映射到第二行和第二列。 x = 0.3会映射到0.5,因为它在0.00和0.60之间。

我知道如何得到我想要的结果,但我确定有一个非常快速/清晰的单行或双行(或已存在的函数)可以做到这一点,使我的代码更清晰。基本上它需要在一些数组之间进行分段插值。

下面是一个例子(主要基于代码从Scipy interpolation on a numpy array) - 我把TODO哪里这个新功能会去:

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 

回答

15

我想你想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

完美。这正是我想要的,谢谢! – 2011-02-28 16:55:16

+0

我发布了一个后续问题,以便您能够帮助我。它在这里:。再次感谢。 – 2011-02-28 18:31:40