Python/Scipy 插值（地图坐标）

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 点并确定插值。我知道我可以用地图坐标来做到这一点。我想知道是否有任何简单/聪明的方法可以将 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]])