Spherical equivalent of Nyquist frequency
I'm a little confused about what you are asking and how it relates to the set of spherical basis functions.  Those are the same one that electron orbitals are based on, right?  I have never understood them to be constrained to the surface of a sphere.  Then again, that is not my usual stomping grounds.

Nevertheless, I think I can generalize your question and provide an answer.

$\phi$ is a poor choice for function name in this context.  Yes, I know it is the standard for "potential function" in many contexts, but within spherical coordinates it has a different meaning.  I will use $lon$ and $lat$ for your point location, and retain $\phi$ for the value so there is no confusion.

You have a large collection of points that fit the model:

$$ \phi_n = \phi(lat_n,lon_n) $$

You would like to find $\phi$ as a linear combination of pre-defined basis functions.

$$ \phi(lat,lon) = \sum_{k=0}^{K-1} c_k B_k(lat,lon) $$

Presumably, $B_0(lat,lon) = 1$, but that is not necessary, but I wanted to throw in the "DC" term to help orient folks in relation to a DFT.

Anyway, for each $k$, you have a computable function definition of $B_k(lon,lat)$.

This is now a standard best fit problem.  You want to find the values of the $c_k$s for the best fit of $\phi$ to your data.

Set up the following arrays:

$$
B = 
\begin{bmatrix}
B_0(lat_0,lon_0) &amp; B_1(lat_0,lon_0) &amp; B_2(lat_0,lon_0) &amp; ...  &amp; B_{K-1}(lat_0,lon_0) \\
B_0(lat_1,lon_1) &amp; B_1(lat_1,lon_1) &amp; B_2(lat_1,lon_1) &amp; ...  &amp; B_{K-1}(lat_1,lon_1) \\ 
: &amp; : &amp; : &amp; ... &amp; : \\
B_0(lat_{N-1},lon_{N-1}) &amp; B_1(lat_{N-1},lon_{N-1}) &amp; B_2(lat_{N-1},lon_{N-1}) &amp; ...  &amp; B_{K-1}(lat_{N-1},lon_{N-1}) \\ 
\end{bmatrix}
$$

$$
C = 
\begin{bmatrix}
c_0 \\
c_1 \\
c_2 \\
: \\
c_{K-1} \\
\end{bmatrix}
$$

$$
P = 
\begin{bmatrix}
\phi_0 \\
\phi_1 \\
: \\
\phi_{N-1} \\
\end{bmatrix}
$$

$B$ (NxK) and $P$ (Nx1) are known from your data points and basis function evaluations.  $C$ (Kx1) contains the unknowns you are solving for.

Ideally,

$$ BC = P $$

If $N=K$, simply multiply both sides by $B^{-1}$ and you have a solution.

$$ C = B^{-1}P $$

Generally, you want more points than that, but that is the minimum requirement unless you want to be underfitted.

If $N&gt;K$, there is a standard approach to find the best fit $C$:

$$ B^* B C = B^* P $$

$$ C = (B^* B)^{-1} B^* P $$

This is basically an interpolation technique.  Interestingly, on the domain of a sphere surface, interpolation vs extrapolation seems to be meaningless.  You will want your points to be as spread out as possible as to not to leave any large areas of your surface unsampled.  The size of $B^* B$ is KxK, and its inversion is probably your biggest computational load depending on how intensive your basis functions are to evaluate at the surface points.

In a conventional DFT, $N=K$ and $B^{-1}$ is $\frac{1}{N}B^*$ so the inverse is implicitly solved.

The concept of a "Nyquist frequency" doesn't have a direct correlation in this scenario in the general case.  It may in terms of your choice of $B_k$ functions, but it would still be irrelevant because in your situation you have randomly distributed points instead of regularly spaced ones.

Your ultimate answer is then you are going to want to have at least as many sample points as you have basis functions, but are likely going to want significantly more.  The number of points you select will not impact the size of the inverse that needs to be taken.

Hope this helps.



This is in response to Fat32's comment.

Looking at the reference given by the OP, it seems the $B_k$s are defined by equations (18) through (33), though it looks like some others are listed as well.  Of course, the list can be extended.  Since the OP's points are randomly distributed, I don't see any other feasible method than the one I provided.  Note that the first one is indeed a "DC" function.

It does raise the question though if there is a set of definable and spread out locations for which matrix $B^{-1}$ is a multiple of $B^*$ when $N=K$, or $ B^* B $ is a multiple of the identity matrix when $N&gt;K$ .  (Calling Olli) Just because continuous functions are orthogonal does not necessarily mean that discrete sampled versions of them are.  With the trig functions (the DFT) they are, with Legendre polynomials, for instance, they aren't. [Correction:  Being diagonal is sufficient for a trivial inverse.]

Regardless, for a fixed set of sample locations $(B^* B)^{-1} B^*$ can always be pre-calculated.  Then finding the coefficients is merely a simple matrix multiplication, just like with the DFT.  This is not what the OP has.  [Correction:  It is possible to select a set of points with $N \ge K$ so the matrix is not invertible.  With a random selection this is not likely, but the columns of $B$ may end up not being very orhogonal which would lead to poor numerical behavior.  Thanks for the paper reference, Olli.  That looks spot on and the OP should definitely read it.]

The concept of a "Nyquist frequency" is strictly dependent on the spacing of evenly spaced sampling locations on a 1D scale, and completely independent of the number of points in the frame.  We tend to think of it as "the last valid basis function before we get into aliases" (assuming the upper half is considered negative frequencies and are used first).  The OP's situation is similar to unevenly sampled points of a signal, so I really don't see the Nyquist concept as being meaningfully translatable.

I interpreted his question as being more of a the necessary number of points for a good answer and just speculation that it was analogous to the Nyquist.

I do think that it is instructive to understand that the DFT is indeed an example of a best fit solution just like this.