Skip to main content

Finding The Intersection Of Two Arcs That Lie On A Sphere

Geoffrey Hunter Author

Finding The Intersection Of Two Arcs That Lie On A Sphere

This is a common problem which usually arises because you want to find if two great circles (or segments of two great circles) on the earth's surface intersect (obviously with are approximating the earth as a sphere).

Let's Start With Our Inputs To The Problem

We will choose to define each of the two arcs with two points on the surface of the sphere/earth. We could use either the geocentric (also called ECEF, coordinates are provided in terms of x, y, z) or geodetic (latitude/longitude) coordinate system to define these points. Let's use geodetic for now:

Points for arc 1 (a1a1):

Pa11=[λ11ϕ11]Pa12=[λ12ϕ12]\begin{align} \b{P_{a11}} &= \left[ {\begin{array}{c} \lambda_{11} \\ \phi_{11} \end{array} } \right] \\ \b{P_{a12}} &= \left[ {\begin{array}{c} \lambda_{12} \\ \phi_{12} \end{array} } \right] \\ \end{align}

Points for arc 2 (a2a2):

Pa21=[λ21ϕ21]Pa22=[λ22ϕ22]\begin{align} \b{P_{a21}} &= \left[ {\begin{array}{c} \lambda_{21} \\ \phi_{21} \end{array} } \right] \\ \b{P_{a22}} &= \left[ {\begin{array}{c} \lambda_{22} \\ \phi_{22} \end{array} } \right] \\ \end{align}

λ\lambda = latitude, in degrees
ϕ\phi = longitude, in degrees

We will then convert these into spherical coordinates with the following equations:

x=cos(λ)sin(ϕ)y=cos(λ)cos(ϕ)z=sin(λ)\begin{align} x = \cos(\lambda) \, \sin(\phi) \\ y = \cos(\lambda) \, \cos(\phi) \\ z = \sin(\lambda) \\ \end{align} Pa11=[x11y11z11]Pa12=[x12y12z12]Pa21=[x21y21z21]Pa22=[x22y22z22]\begin{align} \b{P_{a11}} &= \left[ {\begin{array}{c} x_{11} \\ y_{11} \\ z_{11} \end{array} } \right] \\ \b{P_{a12}} &= \left[ {\begin{array}{c} x_{12} \\ y_{12} \\ z_{12} \end{array} } \right] \\ \b{P_{a21}} &= \left[ {\begin{array}{c} x_{21} \\ y_{21} \\ z_{21} \end{array} } \right] \\ \b{P_{a22}} &= \left[ {\begin{array}{c} x_{22} \\ y_{22} \\ z_{22} \end{array} } \right] \\ \end{align}

Calculate Cross-Products

We then need to calculate a cross-product for each arc, using the two vectors which define the arc as the inputs to the cross-product:

N1=Pa11×Pa12N2=Pa21×Pa22\begin{align} \b{N_1} = \b{P_{a11}} \times \b{P_{a12}} \\ \b{N_2} = \b{P_{a21}} \times \b{P_{a22}} \\ \end{align}

Note how we switched from using "point" to "vector" terminology here. Before we were talking about two "points" that lie on the sphere and define the arc. This is the easiest way to intuitively understand how two points can fully define an arc on the surface of a sphere. But now we can describe them as vectors. When we talk about them as vectors, we are referring to the vector which starts at [0,0,0][0, 0, 0] and passes through this point.

These cross-products will be perpendicular to the plane that the arc lies in. Hence these cross-products are called plane normals.

Calculate The Cross-Product of the Cross-Product

Yup, you read that correctly. We now calculate the cross-product of the two cross-products we calculated above. Because the two arc cross-products define two plane normals, the cross-product of two plane normals gives us a vector that is in the same direction as the intersecting line.

L=N1×N2\b{L} = \b{N_1} \times \b{N_2}

Because the line must pass through [0,0,0]T[0, 0, 0]^T, this vector is also the line itself!

Find The Intersecting Points

Now that we have the line of intersection, we can now find the two points of intersection between the arc planes by normalizing L (remember, we are using a unit sphere, so normalizing the line reduces it to a vector to one of the intersecting points):

I1=LLI2=I1\b{I_1} = \frac{\b{L}}{||\b{L}||} \\ \b{I_2} = -\b{I_1}

To get the other point of intersection, just take the negative of the first one.

Check If These Intersecting Points Are Within The Original Arc Segments

So we now have the two points of intersection between the arcs. If these arcs were complete great circles, then job done, you are guaranteed that the two above points lie within your great circles. But what if you had great-circle arc segments? These intersection points may not be within both arc segments.

We can determine if these intersection points lies within arc segments by doing a angle test, using the dot product:

ab=abcos(θ)\b{a} \cdot \b{b} = ||\b{a}|| ||\b{b}|| cos (\theta)

We will re-arrange for θ\theta:

θ=arccos(abab)\theta = \arccos \left( \frac{\b{a} \cdot \b{b}}{||\b{a}|| ||\b{b}|| } \right)

The angle from the start of arc 1 to intersecting point 1, and the angle from the end of arc 1 to intersecting point 1, and compare it against the angle between the start and end of arc 1.

The intersecting point is within arc 1 if:

θa11,i1+θa12,i1=θa11,a12\theta_{a11,i1} + \theta_{a12,i1} = \theta_{a11,a12}

Take note that if calculating this result using any data type that can lose precision (e.g. floats, doubles), you will have to check it is close to equal rather than exactly equal. This can be done by adding some epsilon. Usually 1101^{-10} is sufficient.

This test has to be repeated between potential intersecting point 1 and arc 2. If both arcs intersect this potential intersecting point, then with have just confirmed that both arcs do intersect at this point!

The same process has to be applied to potential intersecting point 2 (remember, normalizing the line of intersection between the two planes gives us TWO potential points of intersection).


Let's define two arcs a1a1 and a2a2 using two points each in geodetic (lat/lon) form (latitude and longitude are in degrees):

Pa11=[1020]Pa12=[6090]Pa21=[5010]Pa22=[580]\begin{align} P_{a11} = \left[ {\begin{array}{c} 10 \\ 20 \end{array} } \right] P_{a12} = \left[ {\begin{array}{c} 60 \\ 90 \end{array} } \right] \\ P_{a21} = \left[ {\begin{array}{c} 50 \\ 10 \end{array} } \right] P_{a22} = \left[ {\begin{array}{c} 5 \\ 80 \end{array} } \right] \end{align}

Then convert them to spherical coordinates:

a11=[x11y11z11]=[cos(θ)cos(ϕ)cos(θ)sin(ϕ)sin(θ)]=[cos(10)cos(20)cos(10)sin(20)sin(10)]=\pAOneStart\begin{align} \newcommand{\pAOneStart}{\left[ {\begin{array}{c} 0.9254 \\ 0.3368 \\ 0.1736 \end{array} } \right]} \b{a_{11}} &= \left[ {\begin{array}{c} x_{11} \\ y_{11} \\ z_{11} \end{array} } \right] \\ &= \left[ {\begin{array}{c} \cos(\theta) \cos(\phi) \\ \cos(\theta) \sin(\phi) \\ \sin(\theta) \end{array} } \right] \\ &= \left[ {\begin{array}{c} \cos(10) \cos(20) \\ \cos(10) \sin(20) \\ \sin(10) \end{array} } \right] \\ &= \pAOneStart \\ \end{align}

For simplicity, we are working on a unit sphere. If we were not (e.g. working on the Earth's surface, when converting to spherical coordinates we'd have to multiply each x, y and z component by the radius). This of course pretending the Earth is a sphere, which isn't perfect -- but better than thinking it's flat. The other option is to keep working on the unit sphere and scale the intersection point at the end.

And do the same for the other three points:

a12=[3.0616e175.0000e018.6603e01]a21=[0.6330.11160.766]a22=[0.1730.98110.0872]\begin{align} \newcommand{\pAOneEnd}{\left[ {\begin{array}{c} 3.0616e{-}17 \\ 5.0000e{-}01 \\ 8.6603e{-}01 \end{array} } \right]} \b{a_{12}} = \pAOneEnd \b{a_{21}} = \left[ {\begin{array}{c} 0.633 \\ 0.1116 \\ 0.766 \end{array} } \right] \b{a_{22}} = \left[ {\begin{array}{c} 0.173 \\ 0.9811 \\ 0.0872 \end{array} } \right] \end{align}

Now calculate the normal vectors for arc 1 and 2 by taking the cross-product:

N1=a11×a12=\pAOneStart×\pAOneEnd=[0.20490.80140.4627]N2=a21×a22=[0.6330.11160.766]×[0.1730.98110.0872]=[0.74180.07730.6017]\begin{align} \b{N_1} &= \b{a_{11}} \times \b{a_{12}} \\ &= \pAOneStart \times \pAOneEnd \\ &= \left[ {\begin{array}{c} 0.2049 \\ -0.8014 \\ 0.4627 \end{array} } \right] \\ \b{N_2} &= \b{a_{21}} \times \b{a_{22}} \\ &= \left[ {\begin{array}{c} 0.633 \\ 0.1116 \\ 0.766 \end{array} } \right] \times \left[ {\begin{array}{c} 0.173 \\ 0.9811 \\ 0.0872 \end{array} } \right] \\ &= \left[ {\begin{array}{c} -0.7418 \\ 0.0773 \\ 0.6017 \end{array} } \right] \end{align}

Now calculate the "normal of the normals":

L=N1×N2=[0.20490.80140.4627]×[0.74180.07730.6017]=[0.5180.46650.5787]\begin{align} \b{L} &= \b{N_1} \times \b{N_2} \\ &= \left[ {\begin{array}{c} 0.2049 \\ -0.8014 \\ 0.4627 \end{array} } \right] \times \left[ {\begin{array}{c} -0.7418 \\ 0.0773 \\ 0.6017 \end{array} } \right] \\ &= \left[ {\begin{array}{c} -0.518 \\ -0.4665 \\ -0.5787 \end{array} } \right] \end{align}

Now find the two points of intersection of the arc planes:

I1=LL=[8.535e+147.686e+149.533e+14]11.493e15=\iOneI2=I1=\iTwo\begin{align} \newcommand{\iOne}{\left[ {\begin{array}{c} -0.5718 \\ -0.5149 \\ -0.6387 \end{array} } \right]} \b{I_1} &= \frac{ \b{L} }{ || \b{L}|| } \\ &= \left[ {\begin{array}{c} -8.535e+14 \\ -7.686e+14 \\ -9.533e+14 \end{array} } \right] \cdot \frac{1}{ 1.493e15 } \\ &= \iOne \\ \b{I_2} &= -\b{I_1} \\ \newcommand{\iTwo}{\left[ {\begin{array}{c} 0.5718 \\ 0.5149 \\ 0.6387 \end{array} } \right]} &= \iTwo \\ \end{align}

Check if these intersecting points are within the original arc segments. Let's first check if I1\b{I_1} intersects with the arc a1\b{a_1} defined by the points Pa1_start\b{P_{a1\_start}}\quad and Pa1_end\b{P_{a1\_end}}\quad.

θa1_start,i1+θa1_end,i1=θa1_start,a1_end\begin{align} \theta_{a1\_start,i1} + \theta_{a1\_end,i1} = \theta_{a1\_start,a1\_end} \end{align}

We need to calculate θa1_start,i1 \theta_{a1\_start,i1}, θa1_end,i1 \theta_{a1\_end,i1}, and θa1_start,a1_end\theta_{a1\_start,a1\_end}. Lets find θa1_start,i1\theta_{a1\_start,i1} first:

θa1_start,i1=arccos(Pa1_startPi1Pa1_startPi1)=arccos(\pAOneStart\iOne\pAOneStart\iOne)=144.4\begin{align} \theta_{a1\_start,i1} &= \arccos \left( \frac{\b{P_{a1\_start}} \cdot \b{P_{i1}}}{||\b{P_{a1\_start}}|| ||\b{P_{i1}}|| } \right) \\ &= \arccos \left( \frac{ \pAOneStart \cdot \iOne }{|| \pAOneStart || || \iOne || } \right) \\ &= 144.4 \end{align}

We can use the same equation to find θa1_end,i1\theta_{a1\_end,i1}\quad, and θa1_start,a1_end\theta_{a1\_start,a1\_end}\quad:

θa1_end,i1=144.2θa1_start,a2_start=71.4\begin{align} \theta_{a1\_end,i1} &= 144.2 \\ \theta_{a1\_start,a2\_start} &= 71.4 \\ \end{align}

Now we can check the equality:

θa1_start,i1+θa1_end,i1=θa1_start,a1_end144.4+144.2=71.4288.6=71.4\begin{align} \theta_{a1\_start,i1} + \theta_{a1\_end,i1} &= \theta_{a1\_start,a1\_end} \\ 144.4 + 144.2 &= 71.4 \\ 288.6 &= 71.4 \\ \end{align}

Obviously, this equality does not hold true! Therefore, potential intersection point I1\b{I_1} does not lie on the arc a1\b{a_1} and we can rule it out as an intersection point (we do not need to test whether the potential intersection points lies on arc a2\b{a_2}, as if either of the equalities is false, we can immediately rule it out).

Now we test if potential intersection point I2\b{I_2} lies on both a1a_{1} and a2a_{2}.

θa1_start,i2+θa1_end,i2=θa1_start,a1_end35.6+35.8=71.471.4=71.4\begin{align} \theta_{a1\_start,i2} + \theta_{a1\_end,i2} &= \theta_{a1\_start,a1\_end} \\ 35.6 + 35.8 &= 71.4 \\ 71.4 &= 71.4 \\ \end{align}

The equality holds true, I2\b{I_2} lies on the first arc! Lets see if it lies on the second arc:

θa2_start,i2+θa2_end,i2=θa2_start,a2_end24.7+48.7=73.473.4=73.4\begin{align} \theta_{a2\_start,i2} + \theta_{a2\_end,i2} &= \theta_{a2\_start,a2\_end} \\ 24.7 + 48.7 &= 73.4 \\ 73.4 &= 73.4 \\ \end{align}

This equality also holds true, I2\b{I_2} also lies on the second arc! Therefore we know that the point I2=\iTwo\b{I_2} = \iTwo is a valid intersection point between the two arcs. Problem solved!


External Resources is a great page showing how to calculate the intersection between a arc and a Rhumb line (similar to above).