Interpolating non-uniformly distributed points on a 3D sphere

Question

I have several points on the unit sphere that are distributed according to the algorithm described in https://www.cmu.edu/biolphys/deserno/pdf/sphere_equi.pdf (and implemented in the code below). On each of these points, I have a value that in my particular case represents 1 minus a small error. The errors are in [0, 0.1] if this is important, so my values are in [0.9, 1].

Sadly, computing the errors is a costly process and I cannot do this for as many points as I want. Still, I want my plots to look like I am plotting something "continuous". So I want to fit an interpolation function to my data, to be able to sample as many points as I want.

After a little bit of research I found scipy.interpolate.SmoothSphereBivariateSpline which seems to do exactly what I want. But I cannot make it work properly.

Question: what can I use to interpolate (spline, linear interpolation, anything would be fine for the moment) my data on the unit sphere? An answer can be either "you misused scipy.interpolation, here is the correct way to do this" or "this other function is better suited to your problem".

Sample code that should be executable with numpy and scipy installed:

import typing as ty

import numpy
import scipy.interpolate


def get_equidistant_points(N: int) -> ty.List[numpy.ndarray]:
    """Generate approximately n points evenly distributed accros the 3-d sphere.

    This function tries to find approximately n points (might be a little less
    or more) that are evenly distributed accros the 3-dimensional unit sphere.

    The algorithm used is described in
    https://www.cmu.edu/biolphys/deserno/pdf/sphere_equi.pdf.
    """
    # Unit sphere
    r = 1

    points: ty.List[numpy.ndarray] = list()

    a = 4 * numpy.pi * r ** 2 / N
    d = numpy.sqrt(a)
    m_v = int(numpy.round(numpy.pi / d))
    d_v = numpy.pi / m_v
    d_phi = a / d_v

    for m in range(m_v):
        v = numpy.pi * (m + 0.5) / m_v
        m_phi = int(numpy.round(2 * numpy.pi * numpy.sin(v) / d_phi))
        for n in range(m_phi):
            phi = 2 * numpy.pi * n / m_phi
            points.append(
                numpy.array(
                    [
                        numpy.sin(v) * numpy.cos(phi),
                        numpy.sin(v) * numpy.sin(phi),
                        numpy.cos(v),
                    ]
                )
            )
    return points


def cartesian2spherical(x: float, y: float, z: float) -> numpy.ndarray:
    r = numpy.linalg.norm([x, y, z])
    theta = numpy.arccos(z / r)
    phi = numpy.arctan2(y, x)
    return numpy.array([r, theta, phi])


n = 100
points = get_equidistant_points(n)
# Random here, but costly in real life.
errors = numpy.random.rand(len(points)) / 10

# Change everything to spherical to use the interpolator from scipy.
ideal_spherical_points = numpy.array([cartesian2spherical(*point) for point in points])
r_interp = 1 - errors
theta_interp = ideal_spherical_points[:, 1]
phi_interp = ideal_spherical_points[:, 2]
# Change phi coordinate from [-pi, pi] to [0, 2pi] to please scipy.
phi_interp[phi_interp < 0] += 2 * numpy.pi

# Create the interpolator.
interpolator = scipy.interpolate.SmoothSphereBivariateSpline(
    theta_interp, phi_interp, r_interp
)

# Creating the finer theta and phi values for the final plot
theta = numpy.linspace(0, numpy.pi, 100, endpoint=True)
phi = numpy.linspace(0, numpy.pi * 2, 100, endpoint=True)

# Creating the coordinate grid for the unit sphere.
X = numpy.outer(numpy.sin(theta), numpy.cos(phi))
Y = numpy.outer(numpy.sin(theta), numpy.sin(phi))
Z = numpy.outer(numpy.cos(theta), numpy.ones(100))

thetas, phis = numpy.meshgrid(theta, phi)
heatmap = interpolator(thetas, phis)

Issue with the code above:

• With the code as-is, I have a

  ValueError: The required storage space exceeds the available storage space: nxest or nyest too small, or s too small. The weighted least-squares spline corresponds to the current set of knots.
  

  that is raised when initialising the interpolator instance.
• The issue above seems to say that I should change the value of s that is one on the parameters of scipy.interpolate.SmoothSphereBivariateSpline. I tested different values of s ranging from 0.0001 to 100000, the code above always raise, either the exception described above or:

  ValueError: Error code returned by bispev: 10
  

---

Edit: I am including my findings here. They can't really be considered as a solution, that is why I am editing and not posting as an answer.

With more research I found this question Using Radial Basis Functions to Interpolate a Function on a Sphere. The author has exactly the same problem as me and use a different interpolator: scipy.interpolate.Rbf. I changed the above code by replacing the interpolator and plotting:

# Create the interpolator.
interpolator = scipy.interpolate.Rbf(theta_interp, phi_interp, r_interp)

# Creating the finer theta and phi values for the final plot
plot_points = 100
theta = numpy.linspace(0, numpy.pi, plot_points, endpoint=True)
phi = numpy.linspace(0, numpy.pi * 2, plot_points, endpoint=True)

# Creating the coordinate grid for the unit sphere.
X = numpy.outer(numpy.sin(theta), numpy.cos(phi))
Y = numpy.outer(numpy.sin(theta), numpy.sin(phi))
Z = numpy.outer(numpy.cos(theta), numpy.ones(plot_points))

thetas, phis = numpy.meshgrid(theta, phi)
heatmap = interpolator(thetas, phis)


import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib import cm

colormap = cm.inferno
normaliser = mpl.colors.Normalize(vmin=numpy.min(heatmap), vmax=1)
scalar_mappable = cm.ScalarMappable(cmap=colormap, norm=normaliser)
scalar_mappable.set_array([])

fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.plot_surface(
    X,
    Y,
    Z,
    facecolors=colormap(normaliser(heatmap)),
    alpha=0.7,
    cmap=colormap,
)
plt.colorbar(scalar_mappable)
plt.show()

This code runs smoothly and gives the following result:

The interpolation seems OK except on one line that is discontinuous, just like in the question that led me to this class. One of the answer give the idea of using a different distance, more adapted the the spherical coordinates: the Haversine distance.

def haversine(x1, x2):
    theta1, phi1 = x1
    theta2, phi2 = x2
    return 2 * numpy.arcsin(
        numpy.sqrt(
            numpy.sin((theta2 - theta1) / 2) ** 2
            + numpy.cos(theta1) * numpy.cos(theta2) * numpy.sin((phi2 - phi1) / 2) ** 2
        )
    )


# Create the interpolator.
interpolator = scipy.interpolate.Rbf(theta_interp, phi_interp, r_interp, norm=haversine)

which, when executed, gives a warning:

LinAlgWarning: Ill-conditioned matrix (rcond=1.33262e-19): result may not be accurate.
  self.nodes = linalg.solve(self.A, self.di)

and a result that is not at all the one expected: the interpolated function have values that may go up to -1 which is clearly wrong.

Answer 1

You can use Cartesian coordinate instead of Spherical coordinate.

The default norm parameter ('euclidean') used by Rbf is sufficient

# interpolation
x, y, z = numpy.array(points).T
interpolator = scipy.interpolate.Rbf(x, y, z, r_interp)

# predict
heatmap = interpolator(X, Y, Z)

Here the result:

ax.plot_surface(
    X, Y, Z,
    rstride=1, cstride=1, 
    # or rcount=50, ccount=50,
    facecolors=colormap(normaliser(heatmap)),
    cmap=colormap,
    alpha=0.7, shade=False
)
ax.set_xlabel('x axis')
ax.set_ylabel('y axis')
ax.set_zlabel('z axis')

You can also use a cosine distance if you want (norm parameter):

def cosine(XA, XB):
    if XA.ndim == 1:
        XA = numpy.expand_dims(XA, axis=0)
    if XB.ndim == 1:
        XB = numpy.expand_dims(XB, axis=0)
    return scipy.spatial.distance.cosine(XA, XB)

In order to better see the differences, I stacked the two images, substracted them and inverted the layer.