Sunday, October 19, 2014

Understanding Great Circle Arcs Intersection Algorithm

Understanding Great Circle Arcs Intersection Algorithm

The following material is the result of my attempt to understand the nice example from Jason Davies.
I was puzzled about the origin of the algorithm used to find the intersection of two great circle arcs.
Google helped and I discovered Roger Stafford’s post in Matlab newsgroup and the relevant Python’s implementation in the Spherical Geometry Toolkit.

The algorithm

You have two great circle arcs on a sphere, a from point a0 to a1, and b from b0 to b1, whose coordinates are expressed as longitude θ (positive going East of Greenwich) and latitude ϕ (positive going North).
Transform theses coordinates over to Cartesian coordinates using the equations:

xyz=cos(θ)cos(ϕ)=sin(θ)cos(ϕ)=sin(ϕ)

where

πθππ2ϕπ2.

These Cartesian coordinates correspond to a hypothetical spherical “earth” of unit radius, but that does not interfere in the following computations.

Let a0, a1, b0, and b1 be vectors of the Cartesian coordinate endpoints for the two arcs a (a0a1) and b (b0b1) obtained in this way. Carry out the following computations:

p=a0×a1 is the vector normal to the plane going through the arc a and the center of the Earth.

q=b0×b1 is the vector normal to the plane going through the arc b and the center of the Earth.

t=normalized(p×q) is along the line of intersection of the planes above. (The normalization was not mentioned in Roger’s post but it is implemented in the Spherical Geometry Toolkit and by Jason’s example.)

Then, let’s define the following quantities:

s1s2s3s4=(a0×p)t=(a1×p)t=(b0×q)t=(b1×q)t

(These quantities are crucial: they represent the projection of t along the arcs a and b.)
The arcs a and b will intersect s1, s2, s3, and s4 are all of the same sign. In that case they intersect along +t if they are all positive or along t if all are negative.
(Jason tests against ϵ=106, I implemented the test against the sign.)

If they do intersect, you can transform the corresponding vector, t or t, back into longitude and latitude (without worrying about its length.) Letting x, y, z be t’s Cartesian coordinates this reverse transformation can be accomplished this way:

θϕ=arctan2(y,x)=arctan2(z,x2+y2)

References

Written and published with StackEdit.

No comments:

Post a Comment