Bilateration - Intersection of 2 Circles

Question

I'm trying to figure out how to mathematically derive the common points of two intersecting circles on earth's surface given a center Lat/Lon and a radius for each point.

Ex.

Given:

• Lat/Lon (37.673442, -90.234036) Radius 107.5 NM
• Lat/Lon (36.109997, -90.953669) Radius 145 NM

I should find two intersection points with one of them being (36.948, -088.158).

It would be trivially easy to solve this on a flat plane but I don't have any experience solving equations on an imperfect sphere such as the earth' surface. Any help would be appreciated!

Answer 1

It's not much harder on the sphere than on the plane, once you recognize that

1. The points in question are the mutual intersections of three spheres: a sphere centered beneath location x1 (on the earth's surface) of a given radius, a sphere centered beneath location x2 (on the earth's surface) of a given radius, and the earth itself, which is a sphere centered at O = (0,0,0) of a given radius.

2. The intersection of each of the first two spheres with the earth's surface is a circle, which defines two planes. The mutual intersections of all three spheres therefore lies on the intersection of those two planes: a line.

Consequently, the problem is reduced to intersecting a line with a sphere, which is easy.

---

Here are the details. The inputs are points P1 = (lat1, lon1) and P2 = (lat2, lon2) on the earth's surface, considered as a sphere, and two corresponding radii r1 and r2.

1. Convert (lat, lon) to (x,y,z) geocentric coordinates. As usual, because we may choose units of measurement in which the earth has a unit radius,

   x = cos(lon) cos(lat)
   y = sin(lon) cos(lat)
   z = sin(lat).
   

   In the example, P1 = (-90.234036 Degree, 37.673442 Degree) has geocentric coordinates x1 = (-0.00323306, -0.7915, 0.61116) and P2 = (-90.953669 Degree, 36.109997 Degree) has geocentric coordinates x2 = (-0.0134464, -0.807775, 0.589337).

2. Convert the radii r1 and r2 (which are measured along the sphere) to angles along the sphere. By definition, one nautical mile (NM) is 1/60 degree of arc. Therefore, as angles,

   r1 = 107.5 / 60 Degree = 0.0312705 radian
   r2 = 145 / 60 Degree = 0.0421788 radian
   

3. The geodesic circle of radius r1 around x1 is the intersection of the earth's surface with a Euclidean sphere of radius sin(r1) centered at cos(r1)*x1.

4. The plane determined by the intersection of the sphere of radius sin(r1) around cos(r1)*x1 and the earth's surface is perpendicular to x1 and passes through the point cos(r1)*x1, whence its equation is x.x1 = cos(r1) (the "." represents the usual dot product); likewise for the other plane. There will be a unique point x0 on the intersection of those two planes that is a linear combination of x1 and x2. Writing x0 = a*x1 + b*x2 the two planar equations are

   cos(r1) = x.x1 = (a*x1 + b*x2).x1 = a + b*(x2.x1)
   cos(r2) = x.x2 = (a*x1 + b*x2).x2 = a*(x1.x2) + b
   

   Using the fact that x2.x1 = x1.x2, which I shall write as q, the solution (if it exists) is given by

   a = (cos(r1) - cos(r2)*q) / (1 - q^2),
   b = (cos(r2) - cos(r1)*q) / (1 - q^2).
   

   In the running example, I compute a = 0.973503 and b = 0.0260194.

   Evidently we need q^2 != 1. This means that x1 and x2 can be neither the same point nor antipodal points.

5. Now all other points on the line of intersection of the two planes differ from x0 by some multiple of a vector n which is mutually perpendicular to both planes. The cross product

   n = x1~Cross~x2
   

   does the job provided n is nonzero: once again, this means that x1 and x2 are neither coincident nor diametrically opposite. (We need to take care to compute the cross product with high precision, because it involves subtractions with a lot of cancellation when x1 and x2 are close to each other.) In the example, n = (0.0272194, -0.00631254, -0.00803124).

6. Therefore, we seek up to two points of the form x0 + t*n which lie on the earth's surface: that is, their length equals 1. Equivalently, their squared length is 1:

   1 = squared length = (x0 + t*n).(x0 + t*n) = x0.x0 + 2t*x0.n + t^2*n.n = x0.x0 + t^2*n.n
   

   The term with x0.n disappears because x0 (being a linear combination of x1 and x2) is perpendicular to n. The two solutions easily are

   t = sqrt((1 - x0.x0)/n.n)
   

   and its negative. Once again high precision is called for, because when x1 and x2 are close, x0.x0 is very close to 1, leading to some loss of floating point precision. In the example, t = 1.07509 or t = -1.07509. The two points of intersection therefore equal

   x0 + t*n = (0.0257661, -0.798332, 0.601666)
   x0 - t*n = (-0.0327606, -0.784759, 0.618935)
   

7. Finally, we may convert these solutions back to (lat, lon) by converting geocentric (x,y,z) to geographic coordinates:

   lon = ArcTan(x,y)
   lat = ArcTan(Sqrt[x^2+y^2], z)
   

   For the longitude, use the generalized arctangent returning values in the range -180 to 180 degrees (in computing applications, this function takes both x and y as arguments rather than just the ratio y/x; it is sometimes called "ATan2").

   I obtain the two solutions (-88.151426, 36.989311) and (-92.390485, 38.238380), shown in the figure as yellow dots.

The axes display the geocentric (x,y,z) coordinates. The gray patch is the portion of the earth's surface from -95 to -87 degrees longitude, 33 to 40 degree latitude (marked off with a one degree graticule). The earth's surface has been made partly transparent to show all three spheres. The correctness of the computed solutions is evident by how the yellow points sit at the intersections of the spheres.

Answer 2

The ellipsoidal case:

This problem is a generalization of the one of finding maritime boundaries defined as "median lines" and there's an extensive literature on this topic. My solution to this problem is to leverage the equidistant azimuthal projection:

1. Guess at the intersection point
2. Project the two base points using this guessed intersection point as the the center of an equidistant azimuthal projection,
3. Solve the intersection problem in the 2d projected space.
4. It the new intersection point is too far from the old one, go back to step 2.

This algorithm converges quadratically and yields an accurate solution on an ellipsoid. (Accuracy is required in the case of maritime boundaries, since it determines the fishing, oil, and mineral rights.)

Formulas are given in Section 14 of Geodesics on an ellipsoid of revolution. The ellipsoidal equidistant azimuthal projection is provided by GeographicLib. A MATLAB version is available at Geodesic projections for an ellipsoid.