There are many possible ways to explain the why, but I'd like to try one.
Start with some basic definitions for the Poincaré disk model. A point of the hyperbolic plane is a point inside the unit disk. A hyperbolic line is a circle arc which is perpendicular to the unit disk. From that one can eventually deduce that the first four Eulidean axiom as well as the hyperbolic version of the axiom of parallels do hold. So you have a model of hyperbolic geometry.
What you need to study next is the set of isometries of this model. An isometry has to be a bijection between hyperbolic points, so it maps the interior of the unit disk onto itself. Furthermore, an isometry will preserve angles, and since the Poincaré models are conformal, that means Euclidean angles will be preserved.
The set of angle-preserving (conformal) transformations of the plane (or rather the Riemann sphere, but don't worry about that distinction just now) is the set of all Möbius transformations. And one of the most natural ways to represent Möbius transformations is by interpreting points in the plane as complex numbers. (Actually I consider the orojective $\mathbb{CP}^2$ view even more natural, but that's again besides the point.) This is where complex numbers enter the scene. So taken together with mapping the unit disk onto itself, this means that you are after the set of all Möbius transformations which fix the unit disk.
Now one can write down a formula for these, and then do distance comparisons. Suppose you had a hyperbolic line segment from some point $P$ to some nearby point $Q$. You could move that point $P$ into the center, using an isometry which would preserve hyperbolic distances but change Euclidean ones. Applying that transformation to $Q$ would move it close to the center, and then you could measure distances close to the center. You could move your point $Q$ infinitesimally close to $P$, and do the same thing again. What you would get is the distance metric you mentioned, up to a scale factor.
A finite distance is simply an integral over infinitely many infinitesimal distances, so you'd have to integrate from the start point ($z_1$ in your question) to the end point ($z_2$) along a geodesic path. Integrating along geodesics might be tough, though. So you can avoid that by first applying a transformation which moves one of them to the center of the unit disk.
$$
f(z)=\frac{z-z_1}{-\overline{z_1}z+1}\qquad
d\bigl(z_1,z_2\bigr)=d\bigl(0,f(z_2)\bigr)
=2\tanh^{-1}\bigl\lvert f(z_2)\bigr\rvert
=2\tanh^{-1}\bigl\lvert\frac{z_2-z_1}{1-\overline{z_1}z_2}\bigr\rvert
$$
This is the formula you quoted, except for a sign change inside the absolute value which is of no relevance. One can show that all Möbius transformations which fix the unit disk have the form
$$z\mapsto\frac{az+b}{\bar bz+\bar a}$$
so that when you know the numerator for $f$ (which should be zero for $z=z_1$) then you already know the denominator as well. You might choose $a\neq 1$, but the result would simply be the same as above followed by a rotation around the origin.
So now you'd be able to compare distances, like “that hyperbolic segment is twice as long as that other one”. The constant term in that formula is a matter of convention. It is usually chosen in such a way that the Gaussian curvature of the plane is not only constant but actually $-1$. See also this question for details on that convention.
With respect to your question about alternatives, I'd like to point out that cross ratios can be a useful tool there. Take two points in the disk and connect them by a geodesic. That geodesic will intersect the unit circle in two ideal points. Treat all four points as complex numbers and compute their cross ratio. It will not change under Möbius transformations, so it can be used to define distances. Sicne cross ratios are multiplicative and lengths are additive, you need to take the logarithm. Then choose the correct constant as above. The question I referenced there uses this cross ratio approach, so it may be useful here as well.
Note that my above post so far deliberately ignored the fact that Anti-Möbius transformations which fix the unit disk are hyperbolic isometries as well. This is however irrelevant to the question at hand. You should only keep it in mind when you build up your intuition about the set of all hyperbolic isometries.
The simplest way to move between the disc and half-plane models is using complex analysis:
The Mobius transformation $z \mapsto \frac{z+i}{iz+1}$ maps the unit disk to the upper half-plane conformally, inducing an isometry between the Poincare metrics of the half-plane and the unit disk.
Another, more explicit way to think about it is as follows:
Consider the stereographic projection $\varphi: (x,y,z) \mapsto (\frac{x}{1-z}, \frac{y}{1-z})$ from the unit sphere $S^2$ to the plane.
The preimage of the upper half-plane $\{(x,y) : y > 0\}$ under $\varphi$ is the "northern" half-sphere with $y > 0$.
The preimage of the unit disc $\{(x,y) : x^2+y^2<1\}$ is the "bottom" half-sphere with $z < 0$.
Now consider the map $r: (x,y,z) \mapsto (x,-z,y)$ that rotates the unit sphere by $90^\circ$. It maps the bottom half-sphere to the northern half-sphere.
The composition $\varphi \circ r \circ \varphi^{-1}$ maps the Poincare disk to the Poincare half-plane, inducing an isometry of the Poincare metrics.
The relation between these two explanations is the Riemann sphere: The stereographic projection maps complex numbers to their corresponding points on the Riemann sphere, and the particular Mobius transformation $z \mapsto \frac{z+i}{iz+1}$ translates to a rotation of the Riemann sphere by $90^\circ$.
Best Answer
There are two relevant maps between disk and hyperboloid. The one you describe has the disk at one of the vertices, and the center of projection at the center of symmetry of the hyperboloid. The resulting disk is not the Poincaré disk model, but instead the Beltrami-Klein model:
If you want to project to Poincaré disk model instead, you should place your disk at the plane of symmetry (i.e. at $z=0$) and your center of projection at one of the vertices (i.e. $(0,0,-1)$):
One way to reason about these things is the following: A hyperbolic line is uniquely defined by the two ideal points it's incident with, and a point uniquely defined by two lines intersecting with it. So for both of the above mappings, it is sufficient to show that the set of ideal points is projected in a metric-preserving fashion, and that the set of hyperbolic lines matches. In other words, that the set of hyperbolic geodesics on the hyperboloid maps to straight lines for Beltrami-Klein and to circles orthogonal with the unit circle for Poincaré.
Figures taken from my PhD. thesis Section 2.1.5 Figure 2.5 Page 20. Higher resolution available upon request.