迭代插值：首先插值网格，然后插值

Question

我想从x 插值到z。但有一个警告：

根据状态 y，我有不同的 xGrid - 我需要对其进行插值。

我有一个y、yGrid 的网格。说yGrid=[0,1]。而xGrid 由

给出

1    10
2    20
3    30 

对应的zGrid，是

100  1000
200  2000
300  3000

这意味着对于y=0，[1,2,3] 是 x 的正确网格，对于 y=1，[10,20,30] 是正确的网格。和 z 类似。

为了演示问题，一切都是线性且均匀分布的，但它不在实际数据中。

一句话，

• 如果y=0、x=1.5、z 是[1,2,3] 到[100, 200, 300] 在1.5 处的插值，即150。
• 如果y=1、x=10、z=1000

问题来了：如果是y=0.5 呢？在这个简单的例子中，我希望插值网格是 [5.5, 11, 33/2] 和 [550, 1100, 1650]，所以 x=10 将接近 1000。

---

在我看来，我需要插值 3 次：

• 两次得到正确的xGrid和zGrid，并且
• 一次插值xGrid-> xGrid

这是瓶颈的一部分，效率至关重要。我如何最有效地编码？

---

以下是我可以非常低效地编写代码的方法：

xGrid = np.array([[1, 10], [2, 20], [3, 30]])
zGrid = np.array([[100, 1000], [200, 2000], [300, 3000]])
yGrid = np.array([0, 1])
yValue = 0.5
xInterpolated = np.zeros(xGrid.shape[0])
zInterpolated = np.zeros(zGrid.shape[0])
for i in np.arange(xGrid.shape[0]): 
    f1 = interpolate.interp1d(pGrid, xGrid[i,:])
    f2 = interpolate.interp1d(pGrid, zGrid[i,:])
    xInterpolated[i] = f1(yValue)
    zInterpolated[i] = f2(yValue)
f3 = interpolate.interp1d(xInterpolated, zInterpolated)

输出是

In[73]: xInterpolated, zInterpolated
Out[73]: (array([  5.5,  11. ,  16.5]), array([  550.,  1100.,  1650.]))
In[75]: f3(10)
Out[75]: array(1000.0)

---

实际用例数据

xGrid:

array([[  0.30213582,   0.42091889,   0.48596506,   0.55045007,
          0.61479495,   0.67906768,   0.74328653,   0.8074641 ,
          0.8716093 ,   0.93572867,   0.99982708,   1.06390825,
          1.12797508,   1.19202984,   1.25607435,   1.32011008,
          1.38413823,   1.44815978,   1.51217558,   1.57618631],
       [  1.09945362,   1.17100971,   1.23588956,   1.30034354,
          1.36467675,   1.42894086,   1.49315319,   1.55732567,
          1.62146685,   1.68558297,   1.74967873,   1.8137577 ,
          1.87782269,   1.94187589,   2.00591907,   2.06995365,
          2.1339808 ,   2.1980015 ,   2.26201653,   2.32602659],
       [  1.96474476,   2.03281806,   2.09757883,   2.16200519,
          2.22632562,   2.29058026,   2.35478537,   2.41895223,
          2.48308893,   2.54720144,   2.61129424,   2.67537076,
          2.73943368,   2.80348513,   2.86752681,   2.93156011,
          2.99558615,   3.05960586,   3.12362004,   3.18762935],
       [  2.97271432,   3.03917779,   3.10382629,   3.16822546,
          3.23253177,   3.29677589,   3.36097295,   3.42513351,
          3.48926519,   3.55337363,   3.61746308,   3.68153682,
          3.74559741,   3.80964688,   3.87368686,   3.93771869,
          4.00174345,   4.06576206,   4.12977526,   4.1937837 ],
       [  4.17324037,   4.23880534,   4.30336811,   4.36773934,
          4.43202986,   4.49626215,   4.56045011,   4.62460351,
          4.68872947,   4.75283326,   4.81691888,   4.88098942,
          4.94504732,   5.0090945 ,   5.07313252,   5.13716266,
          5.20118595,   5.26520326,   5.32921533,   5.39322276],
       [  5.64337535,   5.70841895,   5.77290336,   5.83724805,
          5.90152063,   5.96573939,   6.02991687,   6.094062  ,
          6.15818132,   6.22227969,   6.28636083,   6.35042763,
          6.41448236,   6.47852685,   6.54256256,   6.60659069,
          6.67061223,   6.73462802,   6.79863874,   6.86264497],
       [  7.51378714,   7.57851747,   7.6429358 ,   7.70725236,
          7.77150412,   7.83570702,   7.89987216,   7.9640075 ,
          8.0281189 ,   8.09221078,   8.15628654,   8.22034883,
          8.28439974,   8.34844097,   8.41247386,   8.47649955,
          8.54051897,   8.60453289,   8.66854195,   8.73254673],
       [ 10.03324294,  10.09777483,  10.162134  ,  10.22641722,
         10.29064401,  10.35482771,  10.41897777,  10.48310105,
         10.54720264,  10.61128646,  10.67535549,  10.73941211,
         10.80345821,  10.8674953 ,  10.93152463,  10.99554722,
         11.05956392,  11.12357544,  11.1875824 ,  11.25158529],
       [ 13.77079831,  13.83519161,  13.89949459,  13.96373623,
         14.02793138,  14.09209044,  14.15622093,  14.2203284 ,
         14.28441705,  14.34849012,  14.41255015,  14.47659914,
         14.54063872,  14.6046702 ,  14.66869465,  14.73271299,
         14.79672596,  14.86073419,  14.92473821,  14.9887385 ],
       [ 20.60440125,  20.66868421,  20.7329108 ,  20.79709436,
         20.8612443 ,  20.92536747,  20.98946899,  21.05355274,
         21.11762172,  21.1816783 ,  21.24572435,  21.30976141,
         21.37379071,  21.43781328,  21.50182995,  21.56584146,
         21.6298484 ,  21.69385127,  21.75785053,  21.82184654]])

zGrid:

array([[ 0.30213582,  0.42091889,  0.48596506,  0.55045007,  0.61479495,
         0.67906768,  0.74328653,  0.8074641 ,  0.8716093 ,  0.93572867,
         0.99982708,  1.06390825,  1.12797508,  1.19202984,  1.25607435,
         1.32011008,  1.38413823,  1.44815978,  1.51217558,  1.57618631],
       [ 0.35871288,  0.43026897,  0.49514882,  0.5596028 ,  0.62393601,
         0.68820012,  0.75241245,  0.81658493,  0.88072611,  0.94484223,
         1.00893799,  1.07301696,  1.13708195,  1.20113515,  1.26517833,
         1.32921291,  1.39324006,  1.45726076,  1.52127579,  1.58528585],
       [ 0.37285697,  0.44093027,  0.50569104,  0.5701174 ,  0.63443782,
         0.69869247,  0.76289758,  0.82706444,  0.89120114,  0.95531365,
         1.01940644,  1.08348296,  1.14754589,  1.21159734,  1.27563902,
         1.33967232,  1.40369835,  1.46771807,  1.53173225,  1.59574155],
       [ 0.38688189,  0.45334537,  0.51799386,  0.58239303,  0.64669934,
         0.71094347,  0.77514053,  0.83930108,  0.90343277,  0.96754121,
         1.03163066,  1.0957044 ,  1.15976498,  1.22381445,  1.28785443,
         1.35188626,  1.41591103,  1.47992963,  1.54394284,  1.60795127],
       [ 0.40252392,  0.46808889,  0.53265166,  0.59702289,  0.66131341,
         0.7255457 ,  0.78973366,  0.85388706,  0.91801302,  0.98211681,
         1.04620243,  1.11027297,  1.17433087,  1.23837805,  1.30241607,
         1.36644621,  1.4304695 ,  1.49448681,  1.55849888,  1.62250631],
       [ 0.42106765,  0.48611125,  0.55059566,  0.61494035,  0.67921293,
         0.74343169,  0.80760917,  0.87175431,  0.93587362,  0.99997199,
         1.06405313,  1.12811993,  1.19217466,  1.25621915,  1.32025486,
         1.38428299,  1.44830454,  1.51232032,  1.57633104,  1.64033728],
       [ 0.4442679 ,  0.50899823,  0.57341657,  0.63773312,  0.70198488,
         0.76618779,  0.83035293,  0.89448826,  0.95859966,  1.02269154,
         1.08676731,  1.15082959,  1.21488051,  1.27892173,  1.34295463,
         1.40698032,  1.47099973,  1.53501365,  1.59902272,  1.66302749],
       [ 0.47525152,  0.53978341,  0.60414258,  0.6684258 ,  0.73265259,
         0.79683629,  0.86098635,  0.92510963,  0.98921122,  1.05329504,
         1.11736407,  1.18142069,  1.24546679,  1.30950388,  1.37353321,
         1.4375558 ,  1.5015725 ,  1.56558403,  1.62959098,  1.69359387],
       [ 0.52099935,  0.58539265,  0.64969564,  0.71393728,  0.77813242,
         0.84229149,  0.90642197,  0.97052944,  1.03461809,  1.09869116,
         1.16275119,  1.22680018,  1.29083976,  1.35487124,  1.4188957 ,
         1.48291403,  1.546927  ,  1.61093523,  1.67493926,  1.73893954],
       [ 0.60440125,  0.66868421,  0.7329108 ,  0.79709436,  0.8612443 ,
         0.92536747,  0.98946899,  1.05355274,  1.11762172,  1.1816783 ,
         1.24572435,  1.30976141,  1.37379071,  1.43781328,  1.50182995,
         1.56584146,  1.6298484 ,  1.69385127,  1.75785053,  1.82184654]])

yGrid:

array([   1.        ,    6.21052632,   11.42105263,   16.63157895,
         21.84210526,   27.05263158,   32.26315789,   37.47368421,
         42.68421053,   47.89473684,   53.10526316,   58.31578947,
         63.52631579,   68.73684211,   73.94736842,   79.15789474,
         84.36842105,   89.57894737,   94.78947368,  100.        ])

我按照给定的答案创建了插值器，然后插值了一些点：

yGrid = yGrid + np.zeros(xGrid.shape)
f3 = interpolate.interp2d(xGrid,yGrid,zGrid,kind='linear')
import matplotlib.pyplot as plt
plt.plot(np.linspace(0.001, 5, 100), [f3(y, 2) for y in np.linspace(0.001, 5, 100)])
plt.plot(xGrid[:, 1], zGrid[:, 1])
plt.plot(xGrid[:, 0], zGrid[:, 0])

这是输出：

蓝线是插值线。我担心对于非常小的 x 值，它应该稍微向下倾斜（遵循两个函数的加权平均值），但根本不是。

Answer 1

您实际上是在查看 2d 插值：您需要 z(x,y) 的插值值为 x 和 y。唯一的微妙之处是您需要广播yGrid 以与x 和z 数据具有相同的形状：

import scipy.interpolate as interpolate

xGrid = np.array([[1, 10], [2, 20], [3, 30]])
zGrid = np.array([[100, 1000], [200, 2000], [300, 3000]])
yGrid = np.array([0, 1]) + np.zeros(xGrid.shape)

yValue = 0.5

f3 = interpolate.interp2d(xGrid,yGrid,zGrid,kind='linear')

这是一个二元函数，你可以这样称呼它

In [372]: f3(10,yValue)
Out[372]: array([ 1000.])

您可以使用 lambda 将其转换为返回标量的单变量函数：

f4 = lambda x,y=yValue: f3(x,y)[0]

这将为您的（假定的）单个 y 值返回单个值，该值在 lambda 定义时设置为 yValue。像这样使用它：

In [376]: f4(10)
Out[376]: 1000.0

不过，通用的f3函数可能更适合您的问题，因为您可以根据需要动态更改y的值，并且可以使用数组输入来获取z的数组输出。

更新

对于形状奇特的x,y 数据，interp2d 可能会产生不令人满意的结果，尤其是在网格的边界处。所以另一种方法是改用interpolate.LinearNDInterpolator，它基于输入数据的三角剖分，本质上给出了局部分段线性近似

f4 = interpolate.LinearNDInterpolator((xGrid.flatten(),yGrid.flatten()),zGrid.flatten())

使用您的更新数据集：

plt.figure()
plt.plot(np.linspace(0.001, 5, 100), f4(np.linspace(0.001, 5, 100), 2))
plt.plot(xGrid[:, 0], zGrid[:, 0])
plt.plot(xGrid[:, 1], zGrid[:, 1])

请注意，这种插值也有其缺点。我建议将两个插值函数都绘制为一个曲面，并查看它们与原始数据相比的失真情况：

from mpl_toolkits.mplot3d import Axes3D

xx,yy=(np.linspace(0,10,20),np.linspace(0,20,40))
xxs,yys=np.meshgrid(xx,yy)
zz3=f3(xx,yy) #from interp2d
zz4=f4(xxs,yys) #from LinearNDInterpolator

#plot raw data
hf=plt.figure()
ax=hf.add_subplot(111,projection='3d')
ax.plot_surface(xGrid,yGrid,zGrid,rstride=1,cstride=1)
plt.draw()

#plot interp2d case
hf=plt.figure()
ax=hf.add_subplot(111,projection='3d')
ax.plot_surface(xxs,yys,zz3,rstride=1,cstride=1)
plt.draw()

#plot LinearNDInterpolator case
hf=plt.figure()
ax=hf.add_subplot(111,projection='3d')
ax.plot_surface(xxs,yys,f4(xxs,yys),rstride=1,cstride=1)
plt.draw()

这将允许您旋转表面并查看它们包含哪些类型的工件（使用适当的后端）。