Let us consider a slightly modified Laplacian $L$ on the unit sphere $\mathbb{S}^2$:
\begin{align*} L&=\frac {1}{\sin\theta}\partial_\theta\sin\theta\partial_\theta+\frac{k(\theta)}{\sin^2\theta}\partial_\varphi. \end{align*}
Above $(x:=\sin\theta\cos\varphi,y:=\sin\theta\sin\varphi)$ with $(\theta,\varphi)\in(0,\pi)\times(0,2\pi)$ are local coordinates of $\mathbb{S}^2$ in which we expressed $L$, and $k:[0,\pi]\mapsto[c,C]$ is a smooth positive function bounded below by $c>0$ and bounded above by $C>c$. When $k\equiv1$ then $L=\Delta_{\mathbb{S}^2}$.
I would like to use oblate spheroidal coordinates to see $L$ as a Laplacian on an oblate spheroid (why? because it would explicitly make $L$ smooth near the poles $\{\theta=0\}$ and $\{\theta=\pi\}$). Using the Wikipedia page, in the coordinates
\begin{align*} x&=a\cosh\mu\cos\nu\cos\varphi,\\ y&=a\cosh\mu\cos\nu\sin\varphi,\\ z&=a\sinh\mu\sin\nu \end{align*}
where $(\mu,\nu)\in\mathbb{R}\times[-\pi/2,\pi/2]$, the Laplacian reads
\begin{align*} \Delta&=\frac{1}{a^2(\sinh^2\mu+\sin^2\nu)} \left[\frac{1}{\cosh \mu}\partial_\mu\cosh\mu\partial_\mu+\frac{1}{\cos\nu}\partial_\nu\cos\nu\partial_\nu\right]+\frac{1}{a^2\cosh^2\mu\cos^2\nu}\partial_\varphi^2. \end{align*}
Taking $\mu$ constant and letting $\theta:=\pi/2-\nu$, we obtain:
\begin{align*} \Delta_{\vert\mu=\text{const}}&=\frac{1}{a^2(\sinh^2\mu+\cos^2\theta)} \left[\frac{1}{\sin\theta}\partial_\theta\sin\theta\partial_\theta+\frac{\sinh^2\mu+\cos^2\theta}{\cosh^2\mu\sin^2\theta}\partial_\varphi^2\right]. \end{align*}
It is not clear to me how to connect with $L$, or if there are better coordinates to explicitly make $L$ smooth at the poles of $\mathbb{S}^2$ (maybe not using oblate spheroidal coordinates?). Any help would be appreciated.
I answer my question in the case where $k(\theta)-k(\alpha)=\mathcal{O}(\theta^2)$ for $\alpha\in\{0,\pi\}$; I do not use oblate spheroidal coordinates.
Near the north pole characterized by $\{\theta=0\}$, we introduce
\begin{align*} \bar{\varphi}&:=\varphi/\sqrt{k(0)} \end{align*}
then
\begin{align*} x&:=\sin\theta\cos\bar{\varphi},\qquad y:=\sin\theta\sin\bar{\varphi}. \end{align*}
We have:
\begin{align*} \begin{pmatrix} \mathrm{d}x\\[1mm]\mathrm{d}y \end{pmatrix}&=\begin{pmatrix} \cos\theta\cos\bar{\varphi}&-y\\[1mm]\cos\theta\sin\bar{\varphi}&x \end{pmatrix}\begin{pmatrix} \mathrm{d}\theta\\[1mm]\mathrm{d}\bar{\varphi} \end{pmatrix},\qquad\begin{pmatrix} \mathrm{d}\theta\\[1mm]\mathrm{d}\bar{\varphi} \end{pmatrix}=\begin{pmatrix} \frac{\cos\bar{\varphi}}{\cos\theta}&\frac{\sin\bar{\varphi}}{\cos\theta}\\[1mm]-\frac{\sin\bar{\varphi}}{\sin\theta}&\frac{\cos\bar{\varphi}}{\sin\theta} \end{pmatrix}\begin{pmatrix} \mathrm{d}x\\[1mm]\mathrm{d}y \end{pmatrix}. \end{align*}
Now $L$ is the Laplace-Beltrami operator associated to the Riemannian metric $g$ defined in polar coordinates by:
\begin{align*} g&=\mathrm{d}\theta^2+\frac{\sin^2\theta}{k(\theta)}\mathrm{d}\varphi^2,\qquad g^{-1}=\partial_\theta^2+\frac{k(\theta)}{\sin^2\theta}\partial\varphi^2. \end{align*}
In the new $(x,y)$ coordinates, we have $\det g = \frac{1}{k(0)\cos^2\theta}$ (with $\cos^2\theta=1-x^2-y^2$ smooth near the north pole $\{x=y=0\}$) and
\begin{align*} g^{xx}&=g^{\theta\theta}\left(\frac{\partial x}{\partial \theta}\right)^2+g^{\varphi\varphi}\left(\frac{\partial\bar{\varphi}}{\partial\varphi}\right)^2\left(\frac{\partial x}{\partial\bar{\varphi}}\right)^2=\cos^2\theta\cos^2\bar{\varphi}+\frac{k(\theta)y^2}{k(0)\sin^2\theta},\\ % g^{yy}&=g^{\theta\theta}\left(\frac{\partial y}{\partial \theta}\right)^2+g^{\varphi\varphi}\left(\frac{\partial\bar{\varphi}}{\partial\varphi}\right)^2\left(\frac{\partial y}{\partial\bar{\varphi}}\right)^2=\cos^2\theta\sin^2\bar{\varphi}+\frac{k(\theta)x^2}{k(0)\sin^2\theta},\\ % g^{xy}&=g^{\theta\theta}\left(\frac{\partial x}{\partial \theta}\right)\left(\frac{\partial y}{\partial \theta}\right)+g^{\varphi\varphi}\left(\frac{\partial\bar{\varphi}}{\partial\varphi}\right)^2\left(\frac{\partial x}{\partial\bar{\varphi}}\right)\left(\frac{\partial y}{\partial\bar{\varphi}}\right)=\cos^2\theta\sin\bar{\varphi}\cos\bar{\varphi}-\frac{k(\theta)xy}{k(0)\sin^2\theta}. \end{align*}
These coefficients are smooth with respect to the $x$ and $y$ variables: for example, writing $k(\theta)=k(0)+\widetilde{k}(\theta)\theta^2$ with $\widetilde{k}\in\mathcal{C}^\infty([0,\pi],\mathbb{R})$,
\begin{align*} g^{xx}&=\cos^2\bar{\varphi}-\sin^2\theta\cos^2\bar{\varphi}+\frac{(k(0)+\widetilde{k}(\theta)\theta^2)y^2}{k(0)\sin^2\theta}\\ &=\cos^2\bar{\varphi}-\sin^2\theta\cos^2\bar{\varphi}+k(0)\sin^2\bar{\varphi}+\frac{\widetilde{k}(\theta)}{k(0)\mathrm{sinc}^2\theta}y^2\\ &=1-x^2+\frac{\widetilde{k}(\theta)}{k(0)\mathrm{sinc}^2\theta}y^2 \end{align*}
which is obviously smooth near the north pole (note that $\frac{\partial\widetilde{k}}{\partial x}=\frac{\cos\bar{\varphi}}{\cos\theta}\frac{\partial\widetilde{k}}{\partial\theta}$ and $\frac{\partial\widetilde{k}}{\partial y}=-\frac{\sin\bar{\varphi}}{\cos\theta}\frac{\partial\widetilde{k}}{\partial\theta}$ are smooth there). Similar computations show that $g^{yy}$ and $g^{xy}$ are also smooth. Then
\begin{align*} L&=\frac{1}{\sqrt{\det g}}\partial_x\sqrt{\det g}g^{xx}\sqrt{\det g}\partial_x+\frac{1}{\sqrt{\det g}}\partial_x\sqrt{\det g}g^{xy}\sqrt{\det g}\partial_y+\frac{1}{\sqrt{\det g}}\partial_y\sqrt{\det g}g^{xy}\sqrt{\det g}\partial_x+\frac{1}{\sqrt{\det g}}\partial_y\sqrt{\det g}g^{y}\sqrt{\det g}\partial_y \end{align*}
is obviously smooth near the north pole. Similar computations show that $L$ is also smooth near the south pole, where we define this time $\bar{\varphi}:=\varphi/\sqrt{k(\pi)}$.