Numerical evaluation of the Laplace operator near the singularity points in spherical coordinates

Question

If one had to evaluate $\Delta Y_l^m$ numerically everywhere on the unit sphere, including the singularity points $\theta = 0,\pi$, how would they do it? Let's say $Y_l^m$ is a spherical harmonic. I'm wondering if there is an "automatic" procedure that avoids this sort of problems. My operator is a bit different, but it also has similar singularities, and in general I don't know, where they're located. Thanks!