I have a set of scalar data corresponding to f(phi0,theta0) = c0 f(phi1,theta1) = c1 f(phi2,theta2) = c2 ... ... The phis and thetas are directions in spherical coordinates and are not regularly spaced. The c's are scalars. I want a f(phi,theta) = linear function of all c's I would imagine I would have to triangulate these points and determine which triangle encompasses an input phi,theta I provide, and then interpolate on the triangle. What are some other methods? This is similar to the 2D planar case.

You could use a weighted average by geodesic distance. This would be linear in c1,...,cN.

Radial basis function interpolation is likely to be well-behaved. Remember to use great circle distance, not euclidean distance.

You could also use spherical harmonics (as many of them as points you have) and find the linear combination of them that matches your data.

By "weighted average by geodesic distance" I guess I'd really meant "weighted average by reciprocal distance." Then my suggestion is a basis function expansion (like Sneftel and alvaro's suggestions), in which the ith basis function is,



where d(p_a, p_b) is the geodesic distance from p_a to p_b for any two points p_a and p_b on the sphere (which is just acos(dot(p_a, p_b)) if we represent points in R^3). Note that,

b_i(p_i) = 1
b_i(p_j) = 0 for all j not equal to i

so you can just use

f(p) = c_1 b_1(p) + ... + c_N b_N(p)

as your interpolating function.

[Edited to correct misinformation. ;-)]