How do I get geographic coordinates of a point if I know the distance from two other points?

Question

The keyword here is geographic.

I am assuming the solution has something to do with a spherical triangle.

I know that this problem has one, two, infinite or no solution at all. My specific problem has two solutions.

---

Here is the approach I am going to use:

1. Get the approximation of the third point by picking it by hand.
2. Calculate distances on ellipsoid
3. Subtract from the original distances
4. Calculate azimuths to the original points
5. Use the subtractions and the azimuths to get a better approximation of the third point.
6. Go back to 2. until the two consecuent approximations don't differ more than the permitted deviation

Answer 1

In this answer I first give a solution to the great circle distance triangulation problem of the OP using spherical trigonometry. I then present a numerical method to solve the geodesic distance triangulation problem on the WGS84 reference ellipsoid (the ellipsoid used in the Global Positioning System). This numerical method requires an initial approximate solution which is then refined using Newton's method as applied to the partial derivatives of a pair of simultaneous equations; that initial estimate can be given by hand or by using the solution found by spherical trigonometry.

In what follows, I'll (mostly) ignore degenerate solutions where 2 (or more) of the points are identical or antipodes.

The great circle / spherical problem

Those simultaneous equations in Emilio Novati's answer are not easy to solve analytically. It is possible to solve those equations using Newton's method, but fortunately this problem can be solved using straight-forward spherical trigonometry.

In what follows I use $\phi$ for latitude and $\lambda$ for longitude; I also use $\overline\phi$ for colatitude, the complement of the latitude i.e., $\overline\phi = 90° - \phi$.

In spherical trigonometry we generally work with a unit sphere and lengths are given in angular measure. The sides of a spherical triangle are arcs of great circles and the length of such an arc is simply the angle that the arc subtends at the centre of the sphere.

In any spherical triangle, given two sides and the included angle the 3rd side can be found using the spherical law of cosines, aka the cos rule.

For a spherical triangle ABC, with angles A, B, and C, and sides $a$, $b$, and $c$ (with side $a$ opposite angle A, etc): $$ \cos a = \cos b \cos c + \sin b \sin c \cos A\\ \cos b = \cos a \cos c + \sin a \sin c \cos B\\ \cos c = \cos a \cos b + \sin a \sin b \cos C\\ $$

The cos rule can also be used to determine angles, given 3 sides: $$\cos A = \frac{\cos a -\cos b \cos c}{\sin b \sin c}$$ etc

Note that in this form $\cos A$ is indeterminate / undefined if $b$ or $c$ equal 0° or 180°. In practice, this is generally not a problem, since those cases arise in simple degenerate triangles.

There's also a spherical law of sines

$$\frac{\sin A}{\sin a}=\frac{\sin B}{\sin b}=\frac{\sin C}{\sin c}$$

The great circle distance triangulation problem

Given the locations of two points A and B find point(s) X such that the great circle distance from A to X is $d_{AX}$ and the great circle distance from B to X is $d_{BX}$, (where the distances are given in angular measure).

The location of A is $\phi_A, \lambda_A$ and the location of B is $\phi_B, \lambda_B$.

Method

First, construct the polar triangle APB, where P is the North pole. Using the cos rule determine $d_{AB}$, the length of side AB, and $A_{azi}$, the angle at A (i.e., $\angle$ PAB) which is the azimuth of side AB at A.
Then construct the triangle AXB and use the cos rule to determine $\alpha$ the angle at A (i.e., $\angle$ XAB). There are (generally) two solutions, one with X "towards" the pole, the other with X "away" from the pole. The azimuth of AX at A = $A_{azi} - \alpha$ for X towards the pole and $A_{azi} + \alpha$ for X away from the pole.
Using that azimuth we can solve triangle APX to get side XP = $\overline{\phi_X}$ and $\angle APX = \lambda_X - \lambda_A$.

Details

In the polar triangle APB the angle at P is $\lambda_B - \lambda_A$. The side AP is $90° - \phi_A = \overline{\phi_A}$, BP = $\overline{\phi_B}$. $d_{AB}$, the length of AB, is found using the cos rule.

$$ \cos d_{AB} = \cos \overline{\phi_A} \cos \overline{\phi_B} + \sin \overline{\phi_A} \sin \overline{\phi_B} \cos(\lambda_B - \lambda_A)\\ \cos d_{AB} = \sin \phi_A \sin \phi_B + \cos \phi_A \cos \phi_B \cos(\lambda_B - \lambda_A) $$

Let $A_{azi}$ be the azimuth of side AB at A (i.e., $\angle$ PAB). Using the cos rule,

$$ \cos A_{azi} = (\cos \overline{\phi_B} - \cos \overline{\phi_A} \cos d_{AB}) / (\sin \overline{\phi_A} \sin d_{AB})\\ \cos A_{azi} = (\sin \phi_B - \sin \phi_B \cos d_{AB}) / (\cos \phi_A \sin d_{AB}) $$

Note that $\cos A_{azi}$ is indeterminate if $d_{AB} = 0$ i.e., if A = B, or if A and B are antipodes (diametrically opposite). $\cos A_{azi}$ is undefined if $\phi_A = 90°$, i.e., if A is at the pole.

We could also use the sin rule to find this angle, since $$\sin P / \sin d_{AB} = \sin A_{azi} / \sin \overline{\phi_B}$$ but cosines (and inverse cosines) are preferred over sines (and inverse sines) if angles greater than 90° may be encountered.

In the triangle AXB, let $\alpha$ be the angle at A (i.e., $\angle$ XAB).

$$ \cos \alpha = (\cos d_{BX} - \cos d_{AX} \cos d_{AB}) / (\sin d_{AX} \sin d_{AB}) $$

$\alpha$ is well-defined as long as neither A=X nor A=B.

The triangle towards the pole.

In the triangle APX, let $X_{azi}$ be the azimuth of AX at A (i.e., $\angle$ PAX).
$X_{azi} = A_{azi} - \alpha$
Side PX is the colatitude of X, $\overline{\phi_X}$.
The polar angle $\angle APX = \lambda_X - \lambda_A$

The latitude of X $$ \cos \overline{\phi_X} = \cos \overline{\phi_A} \cos d_{AX} + \sin \overline{\phi_A} \sin d_{AX} \cos X_{azi}\\ \sin \phi_X = \sin \phi_A \cos d_{AX} + \cos \phi_A \sin d_{AX} \cos X_{azi} $$

The longitude of X $$ \cos (\lambda_X - \lambda_A) = (\cos d_{AX} - \cos \overline{\phi_A} \cos \overline{\phi_X}) / (\sin \overline{\phi_A} \sin \overline{\phi_X})\\ \cos (\lambda_X - \lambda_A) = (\cos d_{AX} - \sin \phi_A \sin \phi_X) / (\cos \phi_A \cos \phi_X) $$

$\cos (\lambda_X - \lambda_A)$ is not well-defined if either A or X is at the pole.

We now have the location of X in the triangle towards the pole.

Repeat the last couple of steps using $X_{azi} = A_{azi} - \alpha$ to get the location of X in the triangle away from the pole.

---

The geodesic / ellipsoidal problem

We can solve the geodesic distance triangulation problem on the WGS84 ellipsoid by using Newton's method as applied to the partial derivatives of a pair of simultaneous equations. I'll first explain that method.

Newton's method for solving simultaneous equations

Let the pair of simultaneous equations be $f(x, y) = 0$ and $g(x, y) = 0$ $$ \Delta f = \frac{\partial f}{\partial x} \Delta x + \frac{\partial f}{\partial y} \Delta y\\ \Delta g = \frac{\partial g}{\partial x} \Delta x + \frac{\partial g}{\partial y} \Delta y\\ $$

Given the partial derivatives and $\Delta f$ & $\Delta g$ we can use the standard techniques of solving linear simultaneous equations to get expressions for $\Delta x$ & $\Delta y$.

Let $$d = \frac{\partial f}{\partial x} \frac{\partial g}{\partial y} - \frac{\partial f}{\partial y} \frac{\partial g}{\partial x}$$ Then $$ d \cdot \Delta x = \Delta f \frac{\partial g}{\partial y} - \frac{\partial f}{\partial y} \Delta g\\ d \cdot \Delta y = \frac{\partial f}{\partial x} \Delta g - \Delta f \frac{\partial g}{\partial x} $$

To solve $f(x, y) = 0$ and $g(x, y) = 0$
Choose initial approximations of $x$ and $y$,
Let $\Delta f = -f(x, y)$ and $\Delta g = -g(x, y)$,
Calculate the partial derivatives for the current $x$ and $y$,
Use the above formulae to obtain $\Delta x$ and $\Delta y$,
Now let $x' = x + \Delta x$ and $y'= y + \Delta y$ be the new approximations of $x$ and $y$.
Repeat until both $\Delta x$ and $\Delta y$ are sufficiently small.

The geodesic distance triangulation problem

Computing geodesic paths on an ellipsoid is a rather more advanced topic than spherical trigonometry. :)

I refer the interested reader to the relevant Wikipedia article Geodesics on an ellipsoid and Geodesics on an ellipsoid of revolution by Charles F. F. Karney, who is the major contributor of that Wikipedia article. Karney's GeographicLib geodesic algorithms are a substantial improvement over the more well-known Vincenty's formulae, which were developed by Thaddeus Vincenty four decades ago, when computer capacity was a lot smaller than it is today. FWIW, GeographicLib is accessible from various programming languages: C / C++, Fortran, Java, JavaScript, Python, .NET, Maxima, MATLAB and Octave.

To quote from that Wikipedia article:

It is possible to reduce the various geodesic problems into one of two types. Consider two points: A at latitude $\phi_1$ and longitude $\lambda_1$ and B at latitude $\phi_2$ and longitude $\lambda_2$ (see Fig. 1). The connecting geodesic (from A to B) is AB, of length $s_{12}$, which has azimuths $\alpha_1$ and $\alpha_2$ at the two endpoints.

The two geodesic problems usually considered are:

the direct geodesic problem or first geodesic problem, given A, $\alpha_1$, and $s_{12}$, determine B and $\alpha_2$;

the inverse geodesic problem or second geodesic problem, given A and B, determine $s_{12}$, $\alpha_1$, and $\alpha_2$.

The direct problem requires solving an elliptic integral. The inverse problem is somewhat harder, since we don't know the initial azimuth $\alpha_1$, and the distance on the ellipsoid depends on the azimuth.

To use Newton's method on the ellipsoid geodesic distance triangulation problem we need a function for the geodesic distance between two points, and the partial derivatives of that function. Although the distance function itself cannot be given in closed form, the partial derivatives (which the distance function depends on) are fairly simple.

Let $a$ be the major axis of the WGS84 ellipsoid (the equatorial radius) and $e^2$ the square of the ellipsoid's eccentricity.
At point $X(\phi, \lambda)$ the meridional radius of curvature $\rho$ is given by
$\rho = a(1 - e^2)/w^3$
and the normal radius of curvature $\nu$ is given by
$\nu = a / w$
where $w = \sqrt{1 - e^2 \sin^2 \phi}$
Then $R = \nu \cos \phi$ is the radius of the circle of latitude at $\phi$.

Let $f(\phi, \lambda)$ be the geodesic distance from point A to X and
$\alpha(\phi, \lambda)$ be the azimuth at X of the path from A to X. (You can calculate these values using Vincenty's formulae, but I recommend using the Geodesic.Inverse() function from Karney's GeographicLib: it converges for all values and it's more accurate than Vincenty's formulae).

Then $$ \frac{\partial f}{\partial \phi} = \rho \cos \alpha\\ \frac{\partial f}{\partial \lambda} = R \sin \alpha $$

---

The Python 2 program GeodesicTriangulate.py implements these algorithms. (I can't paste the program here due to the answer size restriction). It can calculate great circle approximate solutions to feed to the geodesic routine, or you can supply it with approximate solutions on the command line. This program doesn't use many "Pythonisms" apart from some functions returning multiple values, so it should be fairly easy to translate to other languages, if desired.

To run this program you will need to install geographiclib. This can easily be done using pip or see the pypi page for geographiclib.

Occasionally, both great circle approximations will converge to the same point when they're fed to the geodesic routine. In such cases, I suggest try supplying various random positions on the command line - even a point at the North or South pole may work.

Update

This program is now available in a form that is compatible with both Python 2 and Python 3: GeodesicTriangulate.py

Answer 2

Let $A=(\phi_A,\theta_A)$ and $B=(\phi_B,\theta_B)$ the known points, and $X=(\phi,\theta)$ the point of unknown coordinates ( $\phi$ is the longitude and $\theta$ is the latitude).

The arc distance $\alpha$ between $A$ and $X$ is given by:

$$ \cos \alpha= \sin \theta_A\sin \theta+\cos \theta_A\cos \theta\cos(\phi_A-\phi) $$

and the same for the arc distance $BX=\beta$ $$ \cos \beta= \sin \theta_B\sin \theta+\cos \theta_B\cos \theta\cos(\phi_B-\phi) $$

(here there is a proof Great arc distance between two points on a unit sphere, but note that here the angle $\theta$ is that used in spherical coordinates, so it is the complement of the latitude) .

So we have two equations that, simultaneously solved for $(\phi,\theta)$, give the coordinate of $X$.

---

Added after comments.

I'm not expert on geography, so i think that this problem has some standard solution that I don't know. But, as an exercise I'm tempting a possible solution that is more tractable than the previous one.

The idea is to use cartesian coordinates.

1) from the spherical coordinates of $A$ and $B$ find the cartesian coordinates for the points on a sphere of radius $R$, that, without loss of generality, we can take $R=1$. Let $A=(x_A,y_A,z_A)$ and $B=(x_B,y_B,z_B)$.

2) From the arc distances $AX$ and $BX$ find the cartesian distances using the fact that the chord length between two points separated by an arc $\alpha$ in a circle of radius $R$ is $\overline{AX}=2R\sin(\alpha/2)$ .

3) Solve the system $$ \begin{cases} (x-x_A)^2+(y-y_A)^2+(z-z_A)^2=\overline {AX}^2\\ (x-x_B)^2+(y-y_B)^2+(z-z_B)^2=\overline {BX}^2\\ x^2+y^2+z^2=1 \end{cases} $$ Note that since $x_A^2+y_A^2+z_A^2=x_B^2+y_B^2+z_B^2=x^2+y^2+z^2=1$ we can reduce this system to a more tractable one ( but I've not solved it).

4) From the cartesian coordinates $(x,y,z)$ now we can find the spherical coordinates of $X$.

Answer 3

This is not an answer but it is too long for a comment.

Let us consider three equations which write $$F_i=(X-x_i)^2+(Y-y_i)^2+(Z-z_i)^2-d_i^2=0$$ $(i=1,2,3)$ where the unknown variables are $X,Y,Z$.

Write $(F_2-F_1)$ and $(F_3-F_1)$ to get respectively $$2(x_1-x_2)X+2(y_1-y_2)Y+2(z_1-z_2)Z=(x_1^2+y_1^2+z_1^2-d_1^2)-(x_2^2+y_2^2+z_2^2-d_2^2)$$ $$2(x_1-x_3)X+2(y_1-y_3)Y+2(z_1-z_3)Z=(x_1^2+y_1^2+z_1^2-d_1^2)-(x_3^2+y_3^2+z_3^2-d_3^2)$$ These two linear equations allow to express $Y$ and $Z$ as very simple linear expressions of $X$ ($Y=\alpha +\beta X$ ; $Z=\gamma+\delta X$). Plugging the results in equation $F_1$ leads to a quadratic equation in $X$ $$a X^2+2b X + c =0$$ where $$a=1+\beta^2+\delta^2$$ $$b=\alpha \beta +\gamma \delta -(x_1+\beta y_1+\delta z_1)$$ $$c=\alpha ^2+\gamma ^2+(x_1^2+y_1^2+z_1^2-d_1^2)-2 (\alpha y_1+ \gamma z_1)$$