For the case where $|S|=2$ you can assume that $S=\{0,1\}$. It is then well-known that you can take $\pi$ to be the elliptic modular function $\lambda$. It follows easily that it works equally well to take $\pi(z)=1/\lambda(z)$. This induces a conformal isomorphism $\mathbb{H}/G\to\mathbb{C}\setminus S$, where $G$ is the group of matrices $g=\left[\begin{array}{cc}a&b \\ c &d\end{array}\right]\in SL_2(\mathbb{Z})$ with $g=1\pmod{2}$. The action on $\mathbb{H}$ is by $z\mapsto (az+b)/(cz+d)$. If we let $l$ denote the ray $(1,\infty)$ then $\pi^{-1}(l)=\lambda^{-1}((0,1))$, and one can check that this is the union of the $G$-translates of the positive imaginary axis.
If you have Maple you can enter
pi := (z) -> EllipticModulus(exp(Pi*I*z))^(-2);
phi := (w) -> log(EllipticNome(1/sqrt(w)))/(Pi*I);
Then $\pi$ is as discussed above, and $\phi$ is a local inverse for $\pi$. The group $G$ is generated by the transformations $z\mapsto z+2$ and $z\mapsto z/(2z+1)$. Using these you can plot lots of other things if you want.