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.

For example, 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.

  • 1
    If all your radii are going to be this small (less than several kilometers), then the earth is essentially flat at this scale and you might as well choose an accurate, simple projection and perform the usual Euclidean calculations. Make sure you compute the intersection to more than three decimal places--the imprecision in the third decimal place is as large as either of your radii! – whuber Jan 24 '13 at 19:16
  • 1
    I should have added units, those radius are in NM so it's still a small distance relative to the earth's surface but larger than a few km. How does that scale affect distortion? I'm trying to find a solution accurate to less than <1nm, so it doesn't have to be super precise. Thanks! – Will Jan 24 '13 at 20:54
  • This is all good to know, because it shows you can use a spherical model of the earth--the more complicated ellipsoidal models are unnecessary. – whuber Jan 24 '13 at 22:58
  • @whuber Does this imply the problem could be restated as: find the intersection of 3 spheres where one of the spheres is the earth, and the other two are centered on the points with their respective radii? – Kirk Kuykendall Feb 8 '13 at 4:38
  • @Kirk Yes, that's the way to do it, assuming a spherical model of the earth's surface. After some preliminary calculations that reduces this to a special case of the Trilateration problem in 3D. (The calculations are needed to convert distance along spherical arcs to distances along spherical chords, which become the radii of the two smaller spheres.) – whuber Feb 8 '13 at 5:00

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 (which is pi/180 * 1/60 = 0.0002908888 radians). 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 = ax1 + 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.

3D figure

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 degrees 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.

| improve this answer | |
  • Bill, this is awesome. One clarification you could add, based on someone who was trying to implement it. In step 2 you don't explicitly give the conversion from degrees to radians. – Jersey Andy Jun 20 '14 at 21:09
  • @Jersey Thank you for your suggested edit. I changed it a little to avoid the redundancy and to keep the formulas as clear as possible. Having read the thread you are referring to, I also inserted a link to explain the dot product. – whuber Jun 20 '14 at 21:46

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.

| improve this answer | |

Here is some R code to do this:

p1 <- cbind(-90.234036, 37.673442) 
p2 <- cbind(-90.953669, 36.109997 )

steps <- seq(0, 360, 0.1)
c1 <- destPoint(p1, steps, 107.5 * 1852)
c2 <- destPoint(p2, steps, 145 * 1852)

s1 <- spLines(c1)
s2 <- spLines(c2)

i <- intersect(s1, s2)

#        x        y
# -92.38241 38.24267
# -88.15830 36.98740

s <- bind(s1, s2)
crs(s) <- "+proj=longlat +datum=WGS84"
points(i, col='red', pch=20, cex=2)
| improve this answer | |

Following from @whuber's answer, here's some Java code which is useful for two reasons:

  • it highlights a gotcha regarding ArcTan (for Java, and maybe other languages?)
  • it handles the possible edge cases, including one not mentioned in @whuber's answer.

It's not optimised or complete (I've left out obvious classes like Point), but should do the trick.

public static List<Point> intersection(EarthSurfaceCircle c1, EarthSurfaceCircle c2) {

    List<Point> intersections = new ArrayList<Point>();

    // project to (x,y,z) with unit radius
    UnitVector x1 = UnitVector.toPlanar(c1.lat, c1.lon);
    UnitVector x2 = UnitVector.toPlanar(c2.lat, c2.lon);

    // convert radii to radians:
    double r1 = c1.radius / RadiusEarth;
    double r2 = c2.radius / RadiusEarth;

    // compute the unique point x0
    double q = UnitVector.dot(x1, x2);
    double q2 = q * q;
    if (q2 == 1) {
        // no solution: circle centers are either the same or antipodal
        return intersections;
    double a = (Math.cos(r1) - q * Math.cos(r2)) / (1 - q2);
    double b = (Math.cos(r2) - q * Math.cos(r1)) / (1 - q2);
    UnitVector x0 = UnitVector.add(UnitVector.scale(x1, a), UnitVector.scale(x2, b));

    // we only have a solution if x0 is within the sphere - if not,
    // the circles are not touching.
    double x02 = UnitVector.dot(x0, x0);
    if (x02 > 1) {
        // no solution: circles not touching
        return intersections;

    // get the normal vector:
    UnitVector n = UnitVector.cross(x1, x2);
    double n2 = UnitVector.dot(n, n);
    if (n2 == 0) {
        // no solution: circle centers are either the same or antipodal
        return intersections;

    // find intersections:
    double t = Math.sqrt((1 - UnitVector.dot(x0, x0)) / n2);
    intersections.add(UnitVector.toPolar(UnitVector.add(x0, UnitVector.scale(n, t))));
    if (t > 0) {
        // there's only multiple solutions if t > 0
        intersections.add(UnitVector.toPolar(UnitVector.add(x0, UnitVector.scale(n, -t))));
    return intersections;

Also, importantly, note the use of atan2 - it's the reverse of what you'd expect from @whuber's answer (I don't know why, but it works):

    public static Point toPolar(UnitVector a) {
        return new Point(
                Math.toDegrees(Math.atan2(a.z, Math.sqrt(a.x * a.x + a.y * a.y))),
                Math.toDegrees(Math.atan2(a.y, a.x)));          
| improve this answer | |

Working 'R' code for @wuhber answer.

P1 <- c(37.673442, -90.234036)
P2 <- c(36.109997, -90.953669) 

#1 NM nautical-mile is 1852 meters
R1 <- 107.5
R2 <- 145

x1 <- c(
  cos(deg2rad(P1[2])) * cos(deg2rad(P1[1])),  
  sin(deg2rad(P1[2])) * cos(deg2rad(P1[1])),

x2 <- c(
  cos(deg2rad(P2[2])) * cos(deg2rad(P2[1])),
  sin(deg2rad(P2[2])) * cos(deg2rad(P2[1])),

r1 = R1 * (pi/180) * (1/60)
r2 = R2 * (pi/180) * (1/60)

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

n <- cross(x1,x2)

x0 = a*x1 + b*x2

t = sqrt((1 - dot(x0, x0))/dot(n,n))

point1 = x0 + (t * n)
point2 = x0 - (t * n)

lat1 = rad2deg(atan2(point1[2] ,point1[1]))
lon1= rad2deg(asin(point1[3]))
paste(lat1, lon1, sep=",")

lat2 = rad2deg(atan2(point2[2] ,point2[1]))
lon2 = rad2deg(asin(point2[3]))
paste(lat2, lon2, sep=",")
| improve this answer | |

If one of circle is the Nortstar, then there is a easiest way with unit sphere.

You can measure your latitude with Nortstar. Then you have a relative position on this sphere. v1(0,sin(la),cos(la)) You know the position(angle) of the another star (star2), from almanach. v2(sin(lo2)*cos(la2),sin(la2), cos(lo2)*cos(la2)) Its vectors. From equation of sphere.

lo2 is the relative longitude. Its unknow.

The angle between you and star2, you can measure too, (m) And you know, innerproduct of two unit vector is cos(angle) of between. cos(m)=dot(v1,v2) u can calculate now the relative longitude (lo2). lo2=acos((cos(m)-sin(la)*sin(la2))/(cos(la)*cos(la2)))

After all you add the real longitude of star2 to lo2. (or sub, depend on its on west side from you, or east.) lo2 now is your longitude.

Sorry for my english, i never learn this language.

2 things: Northstar mean pole star.

Another. Because the angle measured to horizont relatively, always need 90-angle correction . Its valid to m angle too.

ps: real angle mean: star position - time correction.

| improve this answer | |
  • It is not apparent how this answers the question. – whuber Dec 24 '13 at 15:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.