Python/Scipy 插值 (map_coordinates)

Question

我正在尝试使用scipy进行一些插值操作。我看了很多例子，但还是没找到完全符合我需求的。

假设我有一些数据，行和列的变量可以从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在numpy数组上的插值的代码） - 我在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

Answer 1

我觉得你想要的是一个在矩形网格上使用的双变量样条插值：

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