球体上的 SciPy RectSphereBivariateSpline 插值返回 ValueError

Question

我有一个非常粗糙的球体上的 3D 测量数据，我想进行插值。 我发现 scipy.interpolate 中的 RectSphereBivariateSpline 应该是最合适的。 我使用 RectSphereBivariateSpline 文档中的示例作为起点，现在有以下代码：

""" read csv input file, post process and plot 3D data """
import csv
import numpy as np
from mayavi import mlab
from scipy.interpolate import RectSphereBivariateSpline

# user input
nElevationPoints = 17 # needs to correspond with csv file
nAzimuthPoints = 40 # needs to correspond with csv file
threshold = - 40 # needs to correspond with how measurement data was captured
turnTableStepSize = 72 # needs to correspond with measurement settings
resolution = 0.125 # needs to correspond with measurement settings

# read data from file
patternData = np.empty([nElevationPoints, nAzimuthPoints]) # empty buffer
ifile  = open('ttest.csv') # need the 'b' suffix to prevent blank rows being inserted
reader = csv.reader(ifile,delimiter=',')
reader.next() # skip first line in csv file as this is only text
for nElevation in range (0,nElevationPoints):
    # azimuth
    for nAzimuth in range(0,nAzimuthPoints):  
        patternData[nElevation,nAzimuth] = reader.next()[2]
ifile.close()

# post process
def r(thetaIndex,phiIndex):
    """r(thetaIndex,phiIndex): function in 3D plotting to return positive vector length from patternData[theta,phi]"""
    radius = -threshold + patternData[thetaIndex,phiIndex]
    return radius

#phi,theta = np.mgrid[0:nAzimuthPoints,0:nElevationPoints]
theta = np.arange(0,nElevationPoints)
phi = np.arange(0,nAzimuthPoints)
thetaMesh, phiMesh = np.meshgrid(theta,phi)
stepSizeRad = turnTableStepSize * resolution * np.pi / 180
theta = theta * stepSizeRad
phi = phi * stepSizeRad

# create new grid to interpolate on
phiIndex = np.linspace(1,360,360)
phiNew = phiIndex*np.pi/180
thetaIndex = np.linspace(1,180,180)
thetaNew = thetaIndex*np.pi/180
thetaNew,phiNew = np.meshgrid(thetaNew,phiNew)
# create interpolator object and interpolate
data = r(thetaMesh,phiMesh)
lut = RectSphereBivariateSpline(theta,phi,data.T)
data_interp = lut.ev(thetaNew.ravel(),phiNew.ravel()).reshape((360,180)).T

x = (data_interp(thetaIndex,phiIndex)*np.cos(phiNew)*np.sin(thetaNew))
y = (-data_interp(thetaIndex,phiIndex)*np.sin(phiNew)*np.sin(thetaNew))
z = (data_interp(thetaIndex,phiIndex)*np.cos(thetaNew))

# plot 3D data
obj = mlab.mesh(x, y, z, colormap='jet')
obj.enable_contours = True
obj.contour.filled_contours = True
obj.contour.number_of_contours = 20
mlab.show()

文档中的示例有效，但是当我尝试使用以下测试数据运行上述代码时：testdata 我在声明 RectSphereBivariateSpline 插值器对象的代码位置收到 ValueError：

值错误： 错误：输入时，输入数据受有效性控制 必须满足以下限制。 -1= mumin（见上文），mv >= 4，nuest >=8，nvest >= 8， kwrk>=5+mu+mv+nuest+nvest, lwrk >= 12+nuest*(mv+nvest+3)+nvest*24+4*mu+8*mv+max(nuest,mv+nvest) 0=0 如果 s=0：nuest>=mu+6+iopt(2)+iopt(3)，nvest>=mv+7 如果发现违反这些条件之一，则控制是 立即重新传递给调用程序。在这种情况下，没有 返回近似值。

我已经尝试过，但我完全不知道我应该改变什么以满足 RectSphereBivariateSpline 对象。

有人对我可能做错了什么有任何暗示吗？

-- 编辑-- 根据#HYRY 的建议，我现在可以运行以下代码，并且不会出现运行时错误：

""" read csv input file, post process and plot 3D data """
import csv
import numpy as np
from mayavi import mlab
from scipy.interpolate import RectSphereBivariateSpline

# user input
nElevationPoints = 17 # needs to correspond with csv file
nAzimuthPoints = 40 # needs to correspond with csv file
threshold = - 40 # needs to correspond with how measurement data was captured
turnTableStepSize = 72 # needs to correspond with measurement settings
resolution = 0.125 # needs to correspond with measurement settings

# read data from file
patternData = np.empty([nElevationPoints, nAzimuthPoints]) # empty buffer
ifile  = open('ttest.csv') # need the 'b' suffix to prevent blank rows being inserted
reader = csv.reader(ifile,delimiter=',')
reader.next() # skip first line in csv file as this is only text
for nElevation in range (0,nElevationPoints):
    # azimuth
    for nAzimuth in range(0,nAzimuthPoints):  
        patternData[nElevation,nAzimuth] = reader.next()[2]
ifile.close()

# post process
def r(thetaIndex,phiIndex):
    """r(thetaIndex,phiIndex): function in 3D plotting to return positive vector length from patternData[theta,phi]"""
    radius = -threshold + patternData[thetaIndex,phiIndex]
    return radius

#phi,theta = np.mgrid[0:nAzimuthPoints,0:nElevationPoints]
theta = np.arange(0,nElevationPoints)
phi = np.arange(0,nAzimuthPoints)
thetaMesh, phiMesh = np.meshgrid(theta,phi)
stepSizeRad = turnTableStepSize * resolution * np.pi / 180
theta = theta * stepSizeRad
phi = phi * stepSizeRad

# create new grid to interpolate on
phiIndex = np.arange(1,361)
phiNew = phiIndex*np.pi/180
thetaIndex = np.arange(1,181)
thetaNew = thetaIndex*np.pi/180
thetaNew,phiNew = np.meshgrid(thetaNew,phiNew)
# create interpolator object and interpolate
data = r(thetaMesh,phiMesh)
theta[0] += 1e-6 # zero values for theta cause program to halt; phi makes no sense at theta=0
lut = RectSphereBivariateSpline(theta,phi,data.T)
data_interp = lut.ev(thetaNew.ravel(),phiNew.ravel()).reshape((360,180)).T

def rInterp(theta,phi):
    """rInterp(theta,phi): function in 3D plotting to return positive vector length from interpolated patternData[theta,phi]"""
    thetaIndex = theta/(np.pi/180)
    thetaIndex = thetaIndex.astype(int)
    phiIndex = phi/(np.pi/180)
    phiIndex = phiIndex.astype(int)
    radius = data_interp[thetaIndex,phiIndex]
    return radius
# recreate mesh minus one, needed otherwise the below gives index error, but why??
phiIndex = np.arange(0,360)
phiNew = phiIndex*np.pi/180
thetaIndex = np.arange(0,180)
thetaNew = thetaIndex*np.pi/180
thetaNew,phiNew = np.meshgrid(thetaNew,phiNew)

x = (rInterp(thetaNew,phiNew)*np.cos(phiNew)*np.sin(thetaNew))
y = (-rInterp(thetaNew,phiNew)*np.sin(phiNew)*np.sin(thetaNew))
z = (rInterp(thetaNew,phiNew)*np.cos(thetaNew))

# plot 3D data
obj = mlab.mesh(x, y, z, colormap='jet')
obj.enable_contours = True
obj.contour.filled_contours = True
obj.contour.number_of_contours = 20
mlab.show()

但是，该图与非插值数据有很大不同，请参阅图片here作为参考。

此外，在运行交互式会话时，data_interp 的值 (>3e5) 比原始数据大得多（最大值约为 20）。

还有什么建议吗？

Answer 1

看起来theta[0]不能为0，如果你在调用RectSphereBivariateSpline之前稍微改变一下：

theta[0] += 1e-6