SciPy RectSphereBivariateSpline interpolation over sphere returning ValueError

Question

I have 3D measurement data on a sphere that is very coarse and I want to interpolate. I found that RectSphereBivariateSpline from scipy.interpolate should be most suitable. I used the example in the RectSphereBivariateSpline documentation as a starting point and now have the following code:

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

The example from the documentation works, but when I try to run the above code with the following test data: testdata I get a ValueError at the code position where the RectSphereBivariateSpline interpolator object is declared:

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.

I have tried and tried, but I am absolutely clueless what I should change in order to satisfy the RectSphereBivariateSpline object.

Does anyone have any hint as to what I may be doing wrong?

-- EDIT -- With the suggestions from #HYRY, I now have the following code that runs without runtime errors:

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

However, the plot is much different than the non-interpolated data, see picture here as reference.

Also, when running the interactive session, data_interp is much larger in value (>3e5) than the original data (this is around 20 max).

Any further tips?

Answer 1

It looks like that theta[0] can't be 0, if you change it a litte before call RectSphereBivariateSpline:

theta[0] += 1e-6