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.