Why this approximation for $\pi$ is so accurate?

This isn't actually so mysterious after all; it's just a high precision approximation of the identity $$\int_{-\infty}^\infty e^{-x^2}\,dx = \sqrt{\pi},$$ by using very small rectangles to approximate area. By considering these rectangles, you can also see that the areas cannot be equal, since some area will be left over between the rectangles and the graph of $e^{-x^2}$.

It comes from a theta function identity. In particular, if we let

$$\theta(\tau)=\sum_{n=-\infty}^\infty e^{i\pi\tau n^2}$$

for $\tau$ in the upper half plane of complex numbers, then

$$\theta(-1/\tau)=(-i\tau)^{1/2}\theta(\tau)$$

Letting $\tau=i\pi N^2$ with $N=10^5$, we see that

$$\sum_{n=-\infty}^\infty e^{-n^2/N^2}=\theta(-1/(i\pi N^2))=(\pi N^2)^{1/2}\theta(i\pi N^2)=\sqrt\pi N\sum_{n=-\infty}^\infty e^{-\pi^2N^2n^2}\approx\sqrt\pi N$$

since all terms in the final sum other than $n=0$ contribute extremely small amounts. To be somewhat more precise, we have

$$\begin{align} {1\over N}\sum_{n=-\infty}^\infty e^{-n^2/N^2} &=\sqrt\pi\left(1+2\sum_{n=1}^\infty e^{-\pi^2N^2n^2}\right)\\ &=\sqrt\pi+2\sqrt\pi\left(e^{-\pi^2N^2}+e^{-4\pi^2N^2}+e^{-9\pi^2N^2}+\cdots \right)\\ &\lt\sqrt\pi+2\sqrt\pi\left(e^{-\pi^2N^2}+e^{-2\pi^2N^2}+e^{-3\pi^2N^2}+\cdots \right)\\ &=\sqrt\pi+{2\sqrt\pi e^{-\pi^2N^2}\over1-e^{-\pi^2N^2}} \end{align}$$

and for $N=10^5$, we have $e^{-\pi^2N^2}\approx10^{-42{,}863{,}147{,}299.6}$, from which one can see how we get (more than) $42$ billion digits of $\pi$, answering question (1). We also see that the left hand side of the OP's approximation is slightly smaller than $\pi$, so the formula is not an identity, answering question (2).