Integrating over particular grids to obtain Spherical Harmonic coefficients

Theoretically the spherical harmonic expansion coefficients of a function $f$ should be calculated via a continuous integration: $$F_{lm} = \int_{0}^{2\pi}\int_{0}^{\pi} f(\theta,\phi)Y_{lm}^*(\theta,\phi)\sin(\theta) d\theta d\phi.$$ When only the numerical accuracy of up to a certain order is needed, the coefficients are accessible by numerical integration methods, e.g. using the Gaussian-Legendre quadrature: $$F_{lm}=\frac{\sqrt{2 \pi}}{2N}\sum_{i=0}^{2N-1}\sum_{j=0}^{2N-1}f(\theta_{i},\phi_{j})Y_{lm}^*(\theta_{i},\phi_{j})\omega_{N}(i)$$ In these two equations $\theta$ and $\phi$ are the polar and azimuthal angles, and $\omega$ is the gaussian weight function. Another way to perform the integration is to use equiangular quadratures. However, in my setup, the signal are emitted from the same point and received at equally spaced points, as shown in the following figure: This makes the angular coverage of the data awkwardly spaced, as shown by the blue lines in the following figure for a fixed azimuthal angle, with comparisons with the 6th order Gaussian nodes shown by the black lines:

My question is, what is the most effective way to make use of the data? what I'm thinking of is to interpolate the data to obtain the values at the Gaussian nodes, integrate using the Gaussian scheme to obtain a coarse result; then use the least square method over all the available data points to reduce the error and obtain an optimised final set of result. I would appreciate if you could kindly advise me whether this is the best way forward please, thanks!