SciPy RectSphereBivariate在球体返回值上的Spill插值ValueError

Question

我有一个非常粗糙的球体上的三维测量数据，我想进行插值。 我发现从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我得到的代码位置一个ValueError其中RectSphereBivariateSpline插值对象被声明：

ValueError: ERROR: on entry, the input data are controlled on validity the following restrictions must be satisfied. -1<=iopt(1)<=1, 0<=iopt(2)<=1, 0<=iopt(3)<=1, -1<=ider(1)<=1, 0<=ider(2)<=1, ider(2)=0 if iopt(2)=0. -1<=ider(3)<=1, 0<=ider(4)<=1, ider(4)=0 if iopt(3)=0. mu >= mumin (see above), 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< u(i-1)=0: s>=0 if s=0: nuest>=mu+6+iopt(2)+iopt(3), nvest>=mv+7 if one of these conditions is found to be violated,control is immediately repassed to the calling program. in that case there is no approximation returned.

我已经试了又试，但我绝对没有线索，我应该为了满足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 max）大得多。

还有什么小窍门？

Answer 1

它看起来像theta[0]不能为0，如果调用RectSphereBivariateSpline之前更改豆蔻：

theta[0] += 1e-6