For an upcoming post I wanted to be able to plot the shortest routes between various positions on the Earth using the Bing Maps Silverlight control. Although since I started working on the problem Bing have provided a similar feature with their Distance Calculator App, the functions are not available for reuse via the public API. Interested developers may just want to skip the maths and just download the code.
Geodesic source code for Silverlight 4.0
Geodesics
The shortest path between two points on an arbitrary surface is called a Geodesic, and on a sphere, it is a Great Circle. Modelling the surface of the earth as a perfect sphere, the shortest distance between any two locations on the surface is then described by a section of a Great Circle, ie an arc that lies on the plane that is described by the vectors between its start and end points and the Earth’s centre ( see figure 1 ).
With this information, one way ( and the way I have adopted ) to plot such a curve is as follows:
- Generate the points of the curve in two dimensions using the parametric equation of a circle.
- Transform the plane of the 2d curve into 3D space such that it intersects the end points on the sphere, and the sphere’s centre.
- Project the transformed points back into 2D space using the Mercator projection equations.
[silverlight: Geodesic.xap,520,384,false]
A Parametric Representation of a Great Circle
The Mercator projection gives the 2D rectilinear coordinates (x,y) as a function of the latitude and longitude of a point on a sphere. However, it is easier to draw the point using a typical drawing API, if we have a representation that gives each point of the curve in terms of a single parameter. To derive such a function, we observe that the parametric equation for a circle is given by
\[
\begin{pmatrix} x \\ y \end{pmatrix} = \begin{pmatrix} r \cos t \\ r \sin t \end{pmatrix}
\]
Since a great circle is a rigid transformation of a circle in 3D space, it can also be represented as a function of a single parameter
\[
\begin{pmatrix} x \\ y \\ z \end{pmatrix} = \mathbf{R}\cdot\begin{pmatrix} r \cos t \\ r \sin t \end{pmatrix}
\]
where \(\mathbf{R}\) is a 2×3 matrix that transforms the plane circle to a location on a sphere. Using spherical coordinates, the Mercator projection of the curve specified above is then
\[
\mathbf{C}(t) =\begin{pmatrix} \lambda \\ \tanh^{-1} \left(\sin \phi \right)\end{pmatrix} = \begin{pmatrix} \tan^{-1}\left(y/x\right)\\ \tanh^{-1}z \end{pmatrix}
\]
Where \(\lambda\) and \(\phi\) are the longitude and latitude of the point to be projected, respectively. Writing the equation out in full gives,
\[
\mathbf{C}(t) = \begin{pmatrix}
\tan^{-1}\left(\frac{R_{2,1}\cos t + R_{2,2}\sin t}{R_{1,1}\cos t + R_{1,2}\sin t}\right)\\
\tanh^{-1}(R_{3,1}\cos t +R_{3,2}\sin t)
\end{pmatrix}
\]
Now that we have a suitable parametric equation, we can draw the geodesic with a series of connected line segments by varying the parameter, t.
Mathematical note – the parameterization given by this expression is highly non-uniform, meaning that there are many more points generated in some parts of the curve than in others. The mathematics of generating uniform ( or natural ) parameterizations belongs to the field of differential curve geometry and is beyond the scope of this article ( and my brain ).
Implementation notes
MapShapeBase class to create new shape overlays and secondly, how to handle drawing the curves when they ‘wrap’ beyond the map’s viewable area ( this is easier to illustrate than to describe – see the figure opposite ).
Inheriting from MapShapeBase
I must confess that I cheated slightly in implementing the MapGeodesicPath class, in that I used Reflector to peer into the implementation of the base class. The existing derivations of this class simply defer to MapShapeBase for most of the work, which they can do since for a MapPolygon and MapPolyline there is a one-to-one relationship between Locations, latitude and longitude points on the map and Points, the actual 2D cartesian coordinates used to draw the shape. For the Geodesic, this is not the case, because we only want to specify the start and ending points of the curve, not every point in between. A solution is to delegate the point generation code to a secondary class, one that can be independently tested.
Splitting the curves at the map boundary
When the curves wrap around the map projection, they need to be split at the boundary. This is done by finding the parameter t for the longitude value of the boundary. Where the longitude value is +/- 180, this is straightforward as the equation above reduces to,
\[
t = -\tan^{-1}\frac{R_{1,2}}{R_{2,2}}
\]
For other longitude values, we simply offset the longitude values by the required amount, and calculate the value of t for the new matrix.
