Determining circle-line intersection on sphere

Question

There are two sets of coordinates (latitude/longitude) that form a line along the surface area of the Earth. There is a third coordinate, along with a radius, that denotes a circular region somewhere on the surface of the Earth. I need to determine whether or not the given line intersects with that given region.

Example:

Seattle, Washington's coordinates are (47.6062, -122.3321)

Orlando, Florida's coordinates are (28.5383, -81.3792)

Using Haversine's formula, we know that the distance between them is ~2554 miles.

The state of Wyoming has a radius of approximately 150 miles, and the coordinates of its geographic center are (42.5818, -107.4018), putting it in the path of the line from Seattle to Orlando.

Answer 1

First idea : (probably the simplest) Determine the plane of the great circle and the plane of the small circle and determine if their intersection line $L$ is exterior or not to the sphere (for example by testing if the shortest distance of the origin to line $L$ is larger than radius $R$ or not).

Second idea : Use stereographic projection with direct and inverse formulas you will find there.

In this document, it is also established that the image of circles, either great circles or others, are still circles in the projection plane ; your issue becomes a plane (analytical) geometry question : do two circles intersect, and if such is the case, in which points (using in a further step inverse formulas to obtain the spherical coordinates of these points).

There are exceptional cases (circles passing though the pole) that should (and can) be avoided, for example by doing a preliminary rotation that will be compensated later on.

Answer 2

a third idea is:
- take the plane of the great circle, and thus its normal vector $\bf n$ (or just determine it by the cross product of the position vectors of two points on the "line")
- take the vector of the center of the small circle $\bf c$
- determine the angle $\alpha$ between $\bf c$ and $\bf n$ ($\arccos$ of the normalized dot product)
- $R (\pi/2-\alpha)$ is the distance along the earth surface from the center of the small circle to the great circle considered.

Answer 3

My solution has a few steps, so I would first use G Cab's answer to see if it is worth going further. That is if the distance from G Cab's answer is larger than the small circle radius, then they do not intersect. However, that approach does not tell you if the great circle segment intersects the small circle, and for that a bit more work may be required. This approach has the added benefit of finding the intersection points.

The steps are:

1. Identify the planes normal to the great circle and small circle;
2. Find the line of intersection of the two planes; and
3. Find the points on the line of intersection that lie on the surface of the sphere.

Before we get there, any points are going to need to be transformed to cartesian coordinates. So for example given a latitude $\theta$ and longitude $\phi$, one would need to calculate the $[x,y,z]^T$ equivalent coordinates, e.g.: $$\mathbf{x} = R \left[\begin{array}{c} \cos\theta\cos\phi \\ \cos\theta\sin\phi \\ \sin\theta\end{array}\right],$$ where $R$ is the radius of the sphere. This coordinate system puts the center of the sphere at the origin.

Step 1: Identify the planes normal to the great circle and small circle

Great Circle

Given two points on the surface of the sphere, $\mathbf{p}_1$ and $\mathbf{p}_2$, the vector normal to the great circle is $$\mathbf{n}=\mathbf{p}_1\times\mathbf{p}_2.$$ That normal vector along with either point define the plane. (Below we use $\mathbf{p}_1$ but $\mathbf{p}_2$ could be used just as easily.)

Small Circle

The point $\mathbf{c}$ at the center of the small circle gives us the direction normal to the plane containing it. Knowing that the small circle is also defined based on some distance $r$ along the surface of the sphere, then the angle subtended from the line connecting the origin and $\mathbf{c}$ to a point on the small circle will have an arc $\alpha = r/R$ (in radians). Thus, with some simple trigonometry (SOHCAHTOA), we can show that an arbitrary point in the plane containing the small circle would be $\mathbf{p}_c = \mathbf{c}\cos\alpha = \mathbf{c}\cos(r/R)$.

2. Find the line of intersection of the two planes

From here, we can work through the mathematics to find the intersection of the points by following Krumm (2000). (Note that this part of the solution is based on findings in another post that includes a link to Matlab code if one wishes to use it.)

The direction of the line of intersection of the two planes is $\mathbf{s} = \mathbf{c} \times \mathbf{n}$. This plus a point on the line gives us all we need to define the line of intersection between the planes.

Finding a point on the line is a bit more involved, as one must use Lagrange multipliers. Here we find a point $\mathbf{p}$ that is on the line and closest to an arbitrary point. For simplicity I used the origin, but it is straightforward to choose a different point as described below.

As the point $\mathbf{p}$ is on the line of intersection and thus is also in both planes, we have two constraints: \begin{array}{ll} \displaystyle\left(\mathbf{p}-\mathbf{p}_1\right)\cdot\mathbf{n} = 0; & \mathrm{and} \\ \\ \displaystyle\left(\mathbf{p}-\mathbf{p}_c\right)\cdot\mathbf{c} = 0. \\ \end{array}

Following Krumm (2000), looking to find the point on the line that is closest to the origin, we are looking to minimize $||\mathbf{p}||$ subject to these constraints, which we can combine into the function $w$: $$w = ||\mathbf{p}||^2 + \lambda\left(\mathbf{p}-\mathbf{p}_1\right)\cdot\mathbf{n} + \mu \left(\mathbf{p}-\mathbf{p}_c\right)\cdot\mathbf{c}.$$ We must then take the derivates of $w$ and set them to zero:

$$\begin{array}{l} \displaystyle \frac{\partial w}{\partial p_x} = 2p_x + \lambda n_x + \mu c_x = 0 \\ \\ \displaystyle \frac{\partial w}{\partial p_y} = 2p_y + \lambda n_y + \mu c_y = 0 \\ \\ \displaystyle \frac{\partial w}{\partial p_z} = 2p_z + \lambda n_z + \mu c_x = 0 \\ \\ \displaystyle \frac{\partial w}{\partial \lambda} = \mathbf{p} \cdot \mathbf{n} - \mathbf{p}_1 \cdot \mathbf{n} = 0 \\ \\ \displaystyle \frac{\partial w}{\partial \mu} = \mathbf{p} \cdot \mathbf{c} - \mathbf{p}_c \cdot \mathbf{c} = 0 \end{array}$$ where we have decomposed $\mathbf{p}=\left[p_x, p_y, p_z\right]^T$, $\mathbf{n} = \left[n_x, n_y, n_z\right]$, and $\mathbf{c} = \left[c_x, c_y, c_z\right]^T$. This gives us five equations and five unknowns which can be represented as, $$ \left[\begin{array}{ccccc} 2 & 0 & 0 & n_x & c_x \\ 0 & 2 & 0 & n_y & c_y \\ 0 & 0 & 2 & n_z & c_z \\ n_x & n_y & n_z & 0 & 0 \\ c_x & c_y & c_z & 0 & 0 \end{array}\right] \left[\begin{array}{c} p_x \\ p_y \\ p_z \\ \lambda \\ \mu \end{array}\right] = \left[\begin{array}{c} 0 \\ 0 \\ 0 \\ \mathbf{p}_1 \cdot \hat{\mathbf{n}} \\ \mathbf{p}_c \cdot \hat{\mathbf{c}} \end{array}\right] $$

(Note that if a location other than the origin had been picked, say instead some other point $\mathbf{q}$, then we would instead be looking to minimize $||\mathbf{p}-\mathbf{q}||$ and instead of zeros the first three terms on the right hand side of the equation would be $2q_x$, $2q_y$, and $2q_z$, respectively.)

3. Find the points on the line of intersection that lie on the surface of the sphere

Given the point $\mathbf{p}$ and the direction $\mathbf{s}$, we have an equation for the line $$ \mathbf{g}(t)=\mathbf{p}+\mathbf{s}t, $$ where $t$ is a parametric variable. From here, to see if the great circle crosses into the small circle, one needs just to find the values for $t$ such that $||g||^2=R^2$. Note that this will be quadratic in $t$: $$||\mathbf{s}||^2t^2 + 2(p_xs_x + p_ys_y + p_zs_z)t + ||\mathbf{p}||^2 - R^2 = 0.$$

If two (real) values exist, then the great circle crosses through the small circle. If one exists, then the great circle is tangent to a spot on the small circle. If none exists, then they do not intersect at all.

Finally, say you have two points at $t_1$ and $t_2$ that give you solutions. To identify if they are in the path of the great circle segment between $\mathbf{p}_1$ and $\mathbf{p}_2$, we can define a new parameterization along the great circle by defining $$\mathbf{p}_{1\perp} = \hat{\mathbf{n}}\times\mathbf{p}_1$$ This vector will be in the great circle plane and be perpendicular to $\mathbf{p}_1$, and have the property $||\mathbf{p}_1||=||\mathbf{p}_{1\perp}||=R$.

There exists some angle $0<\beta<2\pi$ such that $$\mathbf{p}_2 = \mathbf{p}_1\cos\beta_{p_2} + \mathbf{p}_{1\perp}\sin\beta_{p_2},$$ and there are similar angles $\beta_{t_1}$ and $\beta_{t_2}$. If $\beta_{t_1}<\beta_{p_2}$, then it is along the great circle segment. Of course, the same is true for $\beta_{t_2}$.

Answer 4

You can find the haversine cross-track distance, which is the closest distance from Point 3 (the center of the circle) to the line. If the cross-track distance is less than or equal to your circle's radius, then you know the circle intersects the line.

In general, the cross-track distance formula is:

$D_{xt} = asin(sin(D_{13})⋅sin(B_{13}−B_{12}))⋅R$

where:

$D_{13}= acos(sin(lat_1)⋅sin(lat_3) + cos(lat_1)⋅cos(lat_3)⋅cos(lon_3-lon_1))$ $B_{13}= atan2(sin(lon_3-lon_1)⋅cos(lat_3), cos(lat_1)⋅sin(lat_3)− sin(lat_1)⋅cos(lat_3)⋅cos(lon_3-lon_1))$ $B_{12}= atan2(sin(lon_2-lon_1)⋅cos(lat_2), cos(lat_1)⋅sin(lat_2)− sin(lat_1)⋅cos(lat_2)⋅cos(lon_2-lon_1))$ $R = EarthRadius$

($D_{13}$ is (angular) distance from 1st line point to 3rd point. $B_{13}$ is (initial) bearing from 1st line point to 3rd point. $B_{12}$ is (initial) bearing from 1st line point to the 2nd line point. $R$ is the earth’s radius, 3958.8 miles)

In Excel, the whole formula looks like this:

=ASIN(SIN(ACOS(SIN(lat_1)*SIN(lat_3) + COS(lat_1)*COS(lat_3)*COS(lon_3-lon_1)))*SIN(ATAN2(COS(lat_1)*SIN(lat_3)-SIN(lat_1)*COS(lat_3)*COS(lon_3-lon_1),SIN(lon_3-lon_1)*COS(lat_3))-ATAN2(COS(lat_1)*SIN(lat_2)-SIN(lat_1)*COS(lat_2)*COS(lon_2-lon_1),SIN(lon_2-lon_1)*COS(lat_2))))*3958.8

You'll need to convert the coordinates to radians before plugging them in for lat_1, lat_3, etc. For example "=RADIANS(47.6062)" yields 0.83088493, so use 0.83088493 for lat_1.

The equation gives us $52.8527$ miles. (I corroborated this in Google Maps.) And since 52.9 miles ≤ 150 miles (your approximate radius of Wyoming), this shows Wyoming is indeed a flyover state−just joking!−from Seattle to Orlando.

Note: The cross-track distance formula will find the nearest distance between the 3rd point and ANY part of the great-circle line, not just between the two points on the line, so be careful there.

More Haversine formulas at https://www.movable-type.co.uk/scripts/latlong.html

closest distance from (42.5818, -107.4018) to path