Great Circle Intersections

I was driving through Northern California a few years back, and was wondering what a particular mountain peak off in the distance was. It occurred to me that if I was able to take a bearing on it from two different points, I’d be able to triangulate its position on a map. I decided to write an iPhone app called WhatThePeak to do this.

Download WhatThePeak in the iTunes App Store.

I’m going to discuss here the math that makes the app work.

Great circle – Wikipedia, the free encyclopedia

The shortest path between two points on the surface of the earth is called a Great Circle. This is why airlines fly routes that, just looking on a flat map seem unintuitive — Chicago to Moscow going over the North Pole, for instance. If you actually look at the path on the globe, you’ll see that it is indeed the shortest way to get there, other than drilling a tunnel through the earth.

Any two distinct great circles will intersect in two places on the earth.

You can define a great circle with a latitude, longitude, and initial heading from that point.

The easy way to find the intersections of two great circles is to use two points on each of the circles to define the planes that the circles pass through. Next, calculate the normal vectors for those planes, using the cross product of the vectors defined by those points. Finally, take the cross product of those two normal vectors. This cross product describes one intersection between the two great circles. Its negative vector describes the other intersection.

As an example, I’ll use two points located in my home town, in tuples of (longitude, latitude, heading):

centennial = (-112.554006, 45.995948, 155.43)

continental = (-112.474351, 45.961395, 278.30)

The first point is looking south-southeast from Centennial blvd, and the second roughly west from continental drive.

Quaternions and spatial rotation – Wikipedia, the free encyclopedia

To turn those points into great circles, I’ll use a technique called quaternion rotation. A quaternion is a 4-element vector which is made up of an axis and a rotation about that axis. You can use it to describe any rotation in 3-space.

For each tuple of coordinate and heading, I’ll first use the coordinate to generate a cartesian (x, y, z) vector in a coordinate system defined on the earth as x-axis positive towards the prime meridian at the equator, y-axis positive towards 90 degrees east at the equator, and the z-axis pointing at the north pole. This vector is used as my axis of rotation. To finish defining the rotation quaternion, I’ll use the negative of the heading, as rotations in compass coordinates are positive in the clockwise direction, but rotations in cartesian space are positive in the counter-clockwise direction.

Next, I’ll use this quaternion to rotate a vector pointing at the north pole (0, 0, 1). The second vector will now be oriented along the great circle defined by the heading.

The cross product of these two vectors is the vector normal to the plane in which the great circle sits.

Generating two normal vectors from both sets of coordinate and heading, I use the cross product to find the normal to these normals, which points to one of the intersections of the two great circles. The other intersection is found by negating the first intersection vector.

The distance between two vectors can be determined by taking the magnitude of the vector produced by subtracting one vector from the other. I compare the distance between each of the intersections and the first coordinate, and choose whichever intersection is closer.

In the diagram below, the equator is yellow, and the prime meridian is green. The red circle is the bearing west from Continental Drive, and the blue circle is the bearing south from Centennial. The light blue and red vectors are the calculated second vectors along the great circles, and the dark vectors are the surface normals of the great circles’ planes. The magenta vector is the chosen intersection.

The top view is looking down from the north pole, and the front view is looking at the earth from 90 degrees west along the equator.

Below are two python functions that generate the great circles and calculate their intersections.

I use the cross product function from numerical python, along with functions from a quaternion library to do quaternion rotations.

gs_orient() and gs_intersect() are the interesting functions here.

import numpy as np
import math as m
import transformations as tran

rad2deg = lambda x: (x/pi)*180
deg2rad = lambda d: (d*pi)/180

def quat(axis, theta):
 return tran.quaternion_about_axis(theta, axis)

def qtransform(tpoint, q):
 tmat = tran.quaternion_matrix(q)
 point = np.array([tpoint[0], tpoint[1], tpoint[2], 0])
 return np.dot(tmat, point)[0:3]

mag = lambda v: np.sqrt(np.vdot(v, v))
norm = lambda v: v/mag(v)

north = np.array([0, 0, 1]) # north pole
#ninety = np.array([1, 0, 0]) # looking at 90 degrees west on the equator
origin = np.array([0,0,0])

def lonlat_xyz(longitude, latitude, r=1):
 lat = deg2rad(latitude)
 #print "lat_d: %s, lat_r: %s" % (latitude, lat)
 lon = deg2rad(longitude)
 #print "lon_d: %s, lon_r: %s" % (longitude, lon)
 x = r * m.cos(lat)*m.cos(lon)
 y = r*m.cos(lat)*m.sin(lon)
 z = r*m.sin(lat)
 return array([x, y, z])

def xyz_lonlat(x, y, z, r=1):
 lat = rad2deg(m.asin(z/r))
 lon = rad2deg(m.atan2(y, x))
 return array([lon, lat])

def decdeg2dms(dd):
 minutes,seconds = divmod(abs(dd)*3600,60)
 degrees,minutes = divmod(minutes,60)
 degrees = degrees if (dd >= 0) else -degrees
 return "%do%d'%f''" % (degrees,minutes,seconds)

def gs_orient(lon, lat, azimuth):
 # find a vector representing a great circle defined by a longitude, latitude, and heading

 # define an axis of rotation based on the vector representation of the (lon, lat)
 p1 = lonlat_xyz(lon, lat)
 # this axis of rotation along with the negated azimuth
 # (compass rotation is clockwise, cartesian rotation is CCW)
 # define a quaternion that can be used for rotation
 trans_quat = quat(p1, deg2rad(-azimuth))
 # use this quaternion to rotate a vector at the north pole to a
 # new position on the great circle defined by the (lon, lat) and azimuth
 p2 = qtransform(north, trans_quat)
 return p2

def gs_intersect(lla1, lla2):
 # find the intersection of two points defined by a longitude, latitude, and heading

 # turn the first point into a vector from a (lon, lat) pair
 p1 = lonlat_xyz(lla1[0], lla1[1])
 # find another vector in the great circle from the first point along the first heading
 gv1 = gs_orient(lla1[0], lla1[1], lla1[2])
 # find the normal to the plane described by these two vectors
 norm1 = np.cross(gv1, p1)

 # do the same for the second point
 p2 = lonlat_xyz(lla2[0], lla2[1])
 gv2 = gs_orient(lla2[0], lla2[1], lla2[2])
 norm2 = np.cross(gv2, p2)

 # An intersection between the two great circles is the cross product
 # of the two planes' normals
 inter = norm(np.cross(norm1, norm2))
 # The other intersection is the first vector negated
 inter2 = -1*inter

 # find a straight-line distance between the first point and both intersections
 dist1 = mag(inter-p1)
 dist2 = mag(inter2-p1)

 # return the (lon_lat pair of the closer intersection)
 if dist1 < dist2:
 int_ll = xyz_lonlat(*inter)
 else:
 int_ll = xyz_lonlat(*inter2)

 return int_ll