Find a point along line on Earth given initial point and distance

Question

I have a line that has a series of vertices and segments. From iterating through the vertices, I can calculate the length of the line (approximately) by determine the great circle distance using the Haversine formula.

My next challenge to calculate a point at a given distance. Initial research shows that in a Cartesian space, I can treat the line segment as a vector, and thus normalize it as shown here. However, this doesn't work due to the curvature of the Earth.

So how would one find a point for my scenario? Alter the Haversine to reverse engineer it?

Answer 1

We can find the coordinates of an intermediate point through repeated application of the law of cosines of spherical trigonometry: $$\cos a = \cos b \cos c + \sin b \sin c \cos A$$
or equivalently $$\cos A = \frac{\cos a - \cos b \cos c}{\sin b \sin c}$$

In the above figure, A is the north pole, and arc BC is the great circle route from B to C. We are given the spherical coordinates of B and C, and the angle a'. We want to compute the spherical coordinates of point C'.

In spherical trigonometry, a side of a triangle is measured by the angle subtended at the center of the sphere. To convert an angle to a distance on the surface of the earth, which is assumed to be a sphere, we would multiply the angle in radians by the radius of the earth. The convention in spherical trigonometry is that a point is identified by its longitude, $\theta$, and its co-latitude, $\phi$. The co-latitude is zero at the north pole, unlike geographic latitude which is $90$ degrees at the north pole.

Let's say the spherical coordinates of B are $(\theta_1, \phi_1)$, those of C are $(\theta_2, \phi_2)$, and we know the angle a'. We want to compute the spherical coordinates of C, $(\theta_3, \phi_3)$. The steps of the computation are as follows: $$\begin{align} A &= \theta_2 - \theta_1 \\ b &= \phi_2 \\ c &= \phi_1 \\ a &= \cos^{-1} (\cos b \cos c + \sin b \sin c \cos A) &&(1) \\ B &= \cos^{-1} \left( \frac{\cos b - \cos a \cos c}{\sin a \sin c} \right) &&(2) \\ b' &= \cos^{-1} (\cos a' \cos c + \sin a' \sin c \cos B) &&(3) \\ A' &= \cos^{-1} \left( \frac{\cos a' - \cos b' \cos c}{\sin b' \sin c} \right) &&(4) \\ \theta_3 &= \theta_1 + A' \\ \phi_3 &= b' \end{align}$$

Equations (1), (2), (3), and (4) are all by applications of the spherical law of cosines, applied at various vertices of the spherical triangles.

A possible computational problem arises at equations (2) and (4) since there is a potential for division by zero. We can get around this by using the atan2 function available in many libraries of mathematical functions and the identity $$\cos^{-1} \frac{a}{b} = \tan^{-1} \frac{\sqrt{b^2-a^2}}{a}$$ The atan2 function is defined by $\text{atan2(y,x)} = \tan^{-1} (y/x)$, with the difference that the case $x=0$ does not cause an error. Using the identity, (2) becomes $$B = \tan^{-1} \left( \frac{\sqrt{(\sin a \sin c)^2 - (\cos b - \cos a \cos c)^2}} {\cos b - \cos a \cos c} \right) $$ which we might code as

x = cos(b) - cos(a) * cos(c)
y = sqrt(max(0, (sin(a) * sin(c))^2 - x^2))
B = atan2(y, x)

The purpose of the max function in the second line of code above is to guard against the possibility that the argument of the sqrt function might be a small negative number due to round-off in floating point computation in cases where the argument is theoretically zero.

Similarly, for computational purposes we would rewrite (4) as $$A' = \tan^{-1} \left( \frac{\sqrt{(\sin b' \sin c)^2 - (\cos a' - \cos b' \cos c)^2}}{\cos a' - \cos b' \cos c} \right) $$ with computer code similar to the code for (2).

Reference: Wikipedia article on spherical trigonometry