This question extends the approach mentioned in Expressing spherical functions with zonal harmonics, which I will briefly restate for ease of reading:
I would like to express any function on the sphere (e.g. $x^2 + xy$ over the 2-sphere) as a sum of zonal harmonics. In the case of functions of one variable, this can be done by taking any function $f(x)$, and noting that it can be expressed as $$ f(x) = \sum_{l=0}^\infty a_l P_l(x) $$ where each $P_l(x)$ is the Legendre Polynomial of order l. We then find coefficients $a_l$ with the fact that $$ a_l = \frac{2l+1}{2} \int_{-1}^1 f(x) P_l(x) dx $$
Numerically, we can build an approximation $g(x)$ by finding: $$ a_l = \frac{2l+1}{2} \int_0^\pi f(\cos(x))P_l(\cos(x)) sin(x) dx $$ and then taking: $$ f(x) \approx g(x) = \sum_{l=0}^\infty a_l P_l(x) $$
So far this is all fine. The issue now is that I would like to express a function on the sphere $f(\theta, \phi)$, and thus the same integral does not work. The ADDITIONAL stipulation is that I cannot take the spherical harmonics equation so that: $$ f(\theta, \phi) = \sum_{l=0}^\infty \sum_{m = -l}^l a_{lm} Y^m_l(\phi, \theta) $$ Instead, I have been advised to form an orthonormal basis of Legendre polynomials by rotating the zonal harmonics. Numerically, this is where I am lost. Thus the question:
How do you approximate a function on the sphere by a linear combination of a rotation of the zonal harmonics $P_l(x)$? Thanks!