python - Iterated Interpolation : First interpolate grids, 然后插值

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

这将允许您旋转表面并查看它们包含什么样的工件(具有适当的后端)。