Python SciPy RectSphereBivariateSpline interpolation generating wrong data?

Question

I seem to have solved it!

For on thing, I tried to extrapolate whereas I could only interpolate this scattered data. SO the new interpolation grid should only go up to theta = 140 degrees or so.

But the most important change is the addition of the parameter s=900 in the RectSphereBivariateSpline call.

So I 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.arange(1,361)
phiNew = phiIndex*np.pi/180
thetaIndex = np.arange(1,141)
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,s=900)
data_interp = lut.ev(thetaNew.ravel(),phiNew.ravel()).reshape((360,140)).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,140)
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
intensity = rInterp(thetaNew,phiNew)
obj = mlab.mesh(x, y, z, scalars = intensity, colormap='jet')
obj.enable_contours = True
obj.contour.filled_contours = True
obj.contour.number_of_contours = 20
mlab.show()

The resulting plot compares nicely to the original non-interpolated data:

enter image description here

I don't fully understand why s should be set at 900, since the RectSphereBivariateSpline documentation says that s=0 for interpolation. However, when reading the documentation a little further I gain some insight:

Chosing the optimal value of s can be a delicate task. Recommended values for s depend on the accuracy of the data values. If the user has an idea of the statistical errors on the data, she can also find a proper estimate for s. By assuming that, if she specifies the right s, the interpolator will use a spline f(u,v) which exactly reproduces the function underlying the data, she can evaluate sum((r(i,j)-s(u(i),v(j)))**2) to find a good estimate for this s. For example, if she knows that the statistical errors on her r(i,j)-values are not greater than 0.1, she may expect that a good s should have a value not larger than u.size * v.size * (0.1)**2. If nothing is known about the statistical error in r(i,j), s must be determined by trial and error. The best is then to start with a very large value of s (to determine the least-squares polynomial and the corresponding upper bound fp0 for s) and then to progressively decrease the value of s (say by a factor 10 in the beginning, i.e. s = fp0 / 10, fp0 / 100, ... and more carefully as the approximation shows more detail) to obtain closer fits.