SciPy：为什么二维插值的输入在内部排序？

iyfamqjs 于 2022-11-10 发布在其他

关注(0)|答案(2)|浏览(137)

我遇到了scipy.interpolate.interp2d问题。
对于已排序的输入，插值正常。
当我要求得到一个未排序数组的插值时，我得到的输出就像是按SciPy内部排序的一样。为什么呢？
方法是在循环中获取插值。
下面是我的演示代码：


# !/usr/bin/env python

# -*- coding: utf-8 -*-

'''
SciPy interp2d test.
Why is the input of 2D interpolation sorted internally?
'''
import matplotlib.pyplot as plt
import scipy.interpolate as itp
import numpy as np

def fMain():

 nx=11
 ny=21

 ax=np.linspace(0,1,nx)
 mx=np.empty((nx,ny))
 for i in range(ny):
  mx[:,i] = ax
 pass

 ay=np.linspace(0,1,ny)
 my=np.empty((nx,ny))
 for i in range(nx): # can I do this without loop?
  my[i,:] = ay
 pass   

 mz=np.empty((nx,ny))
 mz=mx**2 + my**3

 f2Di = itp.interp2d( mx, my, mz, kind='linear')
 #this provides identical results, ok
 #f2Di = itp.interp2d( ax, ay, mz.transpose(), kind='linear')

 if True :
  # just to check the interpolation   
  mzi = f2Di(ax,ay)   
  fig = plt.figure()
  axis = fig.add_subplot(projection='3d')
  axis.plot_wireframe( mx, my, mz )
  axis.scatter(mx, my, mzi.transpose(), marker="o",color="red")
  axis.set_xlabel("x")
  axis.set_ylabel("y")
  axis.set_zlabel("z")
  plt.tight_layout()
  plt.show()
  plt.close()
 pass 

 if True:
  y = 0.5
  az = f2Di( ax , y )
  axf=np.flip(ax)
  azf1 = f2Di( axf , y )
  azf2 = np.empty(nx)
  for i in range(nx):
   azf2[i] = f2Di( axf[i] , y )
  pass    
  plt.plot(ax,az,label="Normal",linewidth=3,linestyle="dashed")
  plt.plot(axf,azf1,label="Reversed")
  plt.plot(axf,azf2,label="Reversed loop")
  plt.legend()
  plt.xlabel("x")
  plt.ylabel("z")
  plt.tight_layout()
  plt.show()
  plt.close()  
 pass
pass    

if __name__ == "__main__":
 fMain()   
pass

scipy

来源：https://stackoverflow.com/questions/72447813/scipy-why-is-the-input-of-2d-interpolation-sorted-internally

2条答案

按热度按时间

0wi1tuuw1#

要回答另一个问题（来自代码中的注解）：

ax=np.linspace(0,1,nx)
mx=np.empty((nx,ny))
for i in range(ny):
    mx[:,i] = ax

(and对于ay类似）。
我可以循环执行此操作吗？
是的（技术上是不可以的，因为在引擎盖下面会有一个C循环，但实际上是可以的）。使用numpy.tile：

ax = np.linspace(0, 1, nx)
mx = np.tile(ax, (ny, 1)).T

而下面的np.empty没有意义：您分配内存，但立即将变量（重新）赋值为另一个值：


# mz = np.empty((nx, ny))   # This line is redundant

mz = mx**2 + my**3

这也是为什么np.empty从for-replacement代码中消失的原因。

赞(0）回复(0）举报 2022-11-10

qkf9rpyu2#

有一个内置的函数可以生成像mx,my这样的网格：

In [68]: I,J = np.meshgrid(ay,ax)
In [69]: I.shape
Out[69]: (11, 21)

In [71]: np.allclose(I,my)
Out[71]: True

In [72]: np.allclose(J,mx)
Out[72]: True

或者，您可以使用广播分配这些值

In [76]: my = np.empty((nx,ny)); my[:]=ay
In [78]: mx = np.empty((nx,ny)); mx[:]=ax[:,None]

interp2d文档说输入数组是扁平的，即使输入是二维的，并且x,y可以是ax,ay中的坐标;它们不必从整个网格中构建。因此，设置f2Di的两种方法是等效的。
使用f2Di(x,y)的完整文档
https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp2d.call.html#scipy.interpolate.interp2d.call
它明确指出输入x,y必须排序，否则它将为您执行此操作。
一次插值：

In [86]: mzi = f2Di(ax,ay)
In [87]: mzi.shape
Out[87]: (21, 11)

另一个输入相反：

In [89]: azf1 = f2Di(ax[::-1], ay[::-1] )
In [90]: azf1.shape
Out[90]: (21, 11)
In [91]: np.allclose(mzi, azf1)
Out[91]: True

正如您所注意到的，并试图用大量绘图代码来显示，结果是相同的-输入在用于插值之前已经排序。
如果我错误地告诉它坐标已排序：

In [94]: azf1 = f2Di(ax[::-1], ay[::-1] , assume_sorted=True)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
...

File ~\anaconda3\lib\site-packages\scipy\interpolate\_fitpack_impl.py:1054, in bisplev(x, y, tck, dx, dy)
   1052 z, ier = _fitpack._bispev(tx, ty, c, kx, ky, x, y, dx, dy)
   1053 if ier == 10:
-> 1054     raise ValueError("Invalid input data")
   1055 if ier:
   1056     raise TypeError("An error occurred")

ValueError: Invalid input data

请注意，这个错误是由_fitpack中的一个函数引发的。这个名字暗示着它使用了某种编译的插值代码，一个可能用C或Fortran编写的库。我不是一个开发人员，但我可以想象，在假设输入是排序的情况下编写这样的代码是最容易的。这样的共享库在对输入有明确和相对简单的期望时工作得最好。

赞(0）回复(0）举报 2022-11-10

我来回答

SciPy：为什么二维插值的输入在内部排序？

2条答案

相关问题

热门标签

最新问答