Calculating intersection of two 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.

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.

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

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.

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.

Answer 3

Here is some R code to do this:

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

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

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

i <- intersect(s1, s2)
coordinates(i)

#        x        y
# -92.38241 38.24267
# -88.15830 36.98740

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

Answer 4

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)));          
    }

Answer 5

Thanks a lot @whuber for your brilliant work! Hereby I share the python code which I typed based on your work, in case it might help someone.

'''
FINDING THE INTERSECTION COORDINATES (LAT/LON) OF TWO CIRCLES (GIVEN THE COORDINATES OF THE CENTER AND THE RADII)

The code below is based on whuber's brilliant work here:
https://gis.stackexchange.com/questions/48937/calculating-intersection-of-two-circles 

The idea is 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.

Note that "Decimal" is used to have higher precision which is important if the distance between two points are a few
meters.
'''
from decimal import Decimal
from math import cos, sin, sqrt
import math
import numpy as np

def intersection(p1, r1_meter, p2, r2_meter):
    # p1 = Coordinates of Point 1: latitude, longitude. This serves as the center of circle 1. Ex: (36.110174,  -90.953524)
    # r1_meter = Radius of circle 1 in meters
    # p2 = Coordinates of Point 2: latitude, longitude. This serves as the center of circle 1. Ex: (36.110174,  -90.953524)
    # r2_meter = Radius of circle 2 in meters
    '''
    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_p1 = Decimal(cos(math.radians(p1[1]))*cos(math.radians(p1[0])))  # x = cos(lon)*cos(lat)
    y_p1 = Decimal(sin(math.radians(p1[1]))*cos(math.radians(p1[0])))  # y = sin(lon)*cos(lat)
    z_p1 = Decimal(sin(math.radians(p1[0])))                           # z = sin(lat)
    x1 = (x_p1, y_p1, z_p1)

    x_p2 = Decimal(cos(math.radians(p2[1]))*cos(math.radians(p2[0])))  # x = cos(lon)*cos(lat)
    y_p2 = Decimal(sin(math.radians(p2[1]))*cos(math.radians(p2[0])))  # y = sin(lon)*cos(lat)
    z_p2 = Decimal(sin(math.radians(p2[0])))                           # z = sin(lat)
    x2 = (x_p2, y_p2, z_p2)
    '''
    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).
    '''
    r1 = Decimal(math.radians((r1_meter/1852) / 60)) # r1_meter/1852 converts meter to Nautical mile.
    r2 = Decimal(math.radians((r2_meter/1852) / 60))
    '''
    3. The geodesic circle of radius r1 around x1 is the intersection of the earth's surface with an 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).
    '''
    q = Decimal(np.dot(x1, x2))

    if q**2 != 1 :
        a = (Decimal(cos(r1)) - Decimal(cos(r2))*q) / (1 - q**2)
        b = (Decimal(cos(r2)) - Decimal(cos(r1))*q) / (1 - q**2)
        '''
        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.)
        '''
        n = np.cross(x1, x2)
        '''
        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
        '''
        x0_1 = [a*f for f in x1]
        x0_2 = [b*f for f in x2]
        x0 = [sum(f) for f in zip(x0_1, x0_2)]
        '''
          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.
        '''
        if (np.dot(x0, x0) <= 1) & (np.dot(n,n) != 0): # This is to secure that (1 - np.dot(x0, x0)) / np.dot(n,n) > 0
            t = Decimal(sqrt((1 - np.dot(x0, x0)) / np.dot(n,n)))
            t1 = t
            t2 = -t

            i1 = x0 + t1*n
            i2 = x0 + t2*n
            '''
            7. Finally, we may convert these solutions back to (lat, lon) by converting geocentric (x,y,z) to geographic
            coordinates. 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").
            '''

            i1_lat = math.degrees( math.asin(i1[2]))
            i1_lon = math.degrees( math.atan2(i1[1], i1[0] ) )
            ip1 = (i1_lat, i1_lon)

            i2_lat = math.degrees( math.asin(i2[2]))
            i2_lon = math.degrees( math.atan2(i2[1], i2[0] ) )
            ip2 = (i2_lat, i2_lon)
            return [ip1, ip2]
        elif (np.dot(n,n) == 0):
            return("The centers of the circles can be neither the same point nor antipodal points.")
        else:
            return("The circles do not intersect")
    else:
        return("The centers of the circles can be neither the same point nor antipodal points.")

'''
Example: The output of below is  [(36.989311051533505, -88.15142628069133), (38.2383796094578, -92.39048549120287)]

         intersection_points = intersection((37.673442, -90.234036), 107.5*1852, (36.109997, -90.953669), 145*1852)
         print(intersection_points)
'''

Any feedback is appreciated.

Answer 6

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])),
  sin(deg2rad(P1[1]))
);

x2 <- c(
  cos(deg2rad(P2[2])) * cos(deg2rad(P2[1])),
  sin(deg2rad(P2[2])) * cos(deg2rad(P2[1])),
  sin(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=",")

Answer 7

Thanks for the great answer, @wuhber @trailblazer. And with the minor change to convert the radii to angles, it matches the matlab scxsc better (sorry, I can't add comment to the original post).

def km2deg(km, radius=6371):
    degree = km*360/(np.pi*radius*2)
    return degree

# e.g. replace the following line in the above python code from @trailblazer with
# r1 = Decimal(math.radians((r1_meter/1852) / 60))
r1 = Decimal(math.radians(km2deg(r1_meter/1000)))

Answer 8

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.