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?