You could do a standard linear least-squares fit to a linear combination of spherical harmonics with variable coefficients:
$$f(\vec x)=\sum_{lm}c_{lm}Y_{lm}(\theta,\phi)$$
$$\sum_{\vec x}\left|f(\vec x)-f_{\vec x}\right|^2\rightarrow \min$$
$$\sum_{\vec x}\left|\sum_{lm}c_{lm}Y_{lm}(\theta,\phi)-f_{\vec x}\right|^2\rightarrow \min$$
$$\sum_{\vec x}Y^*_{lm}(\theta,\phi)\left(\sum_{l'm'}c_{l'm'}Y_{l'm'}(\theta,\phi)-f_{\vec x}\right)=0$$
$$\sum_{l'm'}\left(\sum_{\vec x}Y^*_{lm}(\theta,\phi)Y_{l'm'}(\theta,\phi)\right)c_{lm}=\sum_{\vec x}Y^*_{lm}(\theta,\phi)f_{\vec x}\;,$$
where $\theta$ and $\phi$ are the angular variables corresponding to $\vec x$, $f(\vec x)$ is the fitted function and $f_{\vec x}$ are the function values you have on the grid. This is a system of linear equations for the coefficients $c_{lm}$; you could solve it and say that your function is "nearest" to the spherical harmonic whose coefficient has the highest magnitude.
However, by using only function values within a spherical volume, you can consider the inner sum on the left-hand side as an approximation of an integral:
$$
\begin{eqnarray}
\sum_{\vec x}Y^*_{lm}(\theta,\phi)Y_{l'm'}(\theta,\phi)
&\approx&
\int Y^*_{lm}(\theta,\phi)Y_{l'm'}(\theta,\phi)\mathrm dV
\\
&=&
\iint Y^*_{lm}(\theta,\phi)Y_{l'm'}(\theta,\phi)\mathrm d\Omega r^2\mathrm dr
\\
&=&
\delta_{ll'}\delta_{mm'}\int r^2\mathrm dr
\\
&=&
\frac13\delta_{ll'}\delta_{mm'}R^3\;,
\end{eqnarray}
$$
where $R$ is the radius of the spherical volume. This may actually be quite a good approximation if you have e.g. a cubical grid, since the symmetry of the grid causes many of the coefficients approximated as zero to actually vanish.
This radically simplifies the least-squares solution, since the matrix on the left is now the identity and the coefficients are simply given by the sums on the right.