If you think of your surface as the upper half plane modulo a group of Moebius transformations $G$, start by representing each of your Moebius transformations $ z \longmapsto \frac{az+b}{cz+d}$ by a Matrix.
$$A = \pmatrix{ a & b \\\ c & d}$$
And since only the representative in $PGL_2(\mathbb R)$ matters, people usually normalize to have $Det(A) = \pm 1$.
The standard classification of Moebius transformations as elliptic / parabolic / hyperbolic (loxodromic) is in terms of the determinant and trace squared. You're hyperbolic if and only if the trace squared is larger than $4$. Hyperbolic transformations are the ones with no fixed points in the interior of the Poincare disc, and two fixed points on the boundary, and they are rather explicitly "translation along a geodesic".
Elliptic transformations fix a point in the interior of the disc so they can't be covering transformations. Parabolics you only get as covering transformations if the surface is non-compact, because parabolics have one fixed point and its on the boundary -- if you had such a covering transformation it would tell you your surface has non-trivial closed curves such that the length functional has no lower bound in its homotopy class.
So your covering tranformations are only hyperbolic. That happens only when $tr(A)^2 > 4$. So how do you find your axis? It's the geodesic between the two fixed points on the boundary, so you're looking for solutions to the equation:
$$ t = \frac{at+b}{ct+d}$$
for $t$ real, this is a quadratic equation in the real variable $t$. If I remember the quadratic equation those two points are:
$$ \frac{tr(A) \pm \sqrt{tr(A)^2 - 4Det(A)}}{2c}$$
Is this what you're after?
I don't know that you can make a good distinction between the Riemann uniformization theorem, and the theorem that a surface has a constant-curvature metric in the same conformal class. Each of these results is a relatively easy corollary of the other one, and you can view them as the same theorem.
The old methods of Riemann and Weierstrass give you one difficult proof of the uniformization theorem. Ricci flow is another, beautiful proof of the uniformization theorem. I find it conceptually easier, but it's still a formidable proof. The Ricci flow proof requires various non-trivial results in differential geometry and parabolic PDEs.
Your proof #2 is indeed a proof of a much simpler fact. I think that the main lesson of this proof is the great distance between the existence of a constant-curvature metric on a surface, and the bijection between all such metrics and conformal structures.
Best Answer
If you just want examples for which it's not hard to figure out how the geodesics behave, here's a class of examples with negative and non-constant curvauture in the plane where the geodesics are relatively easy to understand:
Let $a$ and $b$ be smooth functions on $\mathbb{R}$ such that $a(x)+b(y)>0$ for all $x,y\in\mathbb{R}$ and consider the metric $$ g = \bigl(a(x) + b(y)\bigr)(\mathrm{d}x^2 + \mathrm{d}y^2) $$ on $\mathbb{R}^2$. The curvature of this metric is $$ K = \frac{a'(x)^2+b'(y)^2-\bigl(a''(x)+b''(y)\bigr)\bigl(a(x)+b(y)\bigr)} {2\,\bigl(a(x)+b(y)\bigr)^3} $$ It's easy to choose $a$ and $b$ so that $K<0$. For example, $a(x) = x^2+1$ and $b(y) = y^2+1$ gives a complete metric on $\mathbb{R}^2$ that has non-constant negative curvature $K = -4/(x^2{+}y^2{+}2)^{3}<0$.
Note that, taking $a$ (respectively, $b$) to be a constant gives a metric $g$ that has a Killing vector field, namely $\partial/\partial x$ (respectively, $\partial/\partial y$), but, for generic choices of $a$ and $b$, the metric $g$ will have no Killing vector field.
As for geodesics, the good thing about these metrics (called Liouville metrics in the literature) is that their geodesic flows are integrable: Any unit speed geodesic $(x(t),y(t))$ satisfies $$ \bigl(a(x)+b(y)\bigr)\bigl(\dot x^2+\dot y^2\bigr) = 1 \quad\text{and}\quad \bigl(a(x)+b(y)\bigr)\bigl(b(y)\,\dot x^2- a(x)\,\dot y^2\bigr) = c $$ for some constant $c$. (Note that, when either $a$ or $b$ is constant, this second 'first integral' of the geodesic equations specializes to the well-known 'Clairaut integral' for surfaces of revolution.)
In particular, $$ \bigl(b(y)-c)\bigr)\,\dot x^2 - \bigl(a(x)+c)\bigr)\,\dot y^2 = 0, $$ and, assuming that you are in a region when $a(x){+}c$ and $b(y){-}c$ are both positive, $$ \frac{\mathrm{d}x}{\sqrt{a(x)+c}} \pm \frac{\mathrm{d}y}{\sqrt{b(y)-c}}=0, $$ which gives two foliations of this region by geodesics, which can be found by quadrature.
In any case, you will have good qualitative control over these geodesics and can draw some nice pictures.
Added remark (12 May 2020): As an example of what one can do with this more explicit information, you might be interested in this answer of mine to an old question about Riemannian surfaces for which one can compute an explicit distance function.