15,031,219 members
Articles / Programming Languages / C#
Article
Posted 8 Aug 2007

84.8K views
79 bookmarked

# GPS Receivers, Geodesy, and Geocaching: Vincenty’s Formula

Rate me:
Vincenty's Formula is an iterative solution for calculating the distance and direction between two points along the surface of Earth.

## Introduction

Vincenty's Formula is an iterative solution for calculating the distance and direction between two points along the surface of Earth. For clarity, I've stripped out portions of the code I've put up for discussion, but you can download the entire C# source code from the link at the top of this article.

## Background

Several years ago, I stumbled on a great pastime called "geocaching[^]." It's a worldwide treasure hunting game where participants use handheld GPS receivers to find hidden "caches" - small boxes filled with prizes, trinkets, and "travel bugs[^]". The caches are hidden by other participants who post nothing more than the latitude and longitude on a website like Geocaching.com[^]. My children and I have had a blast. It's a great way for a grown man to justify playing in the woods (and buying an expensive gadget!) under the pretence of "playing with the kids." With over 420,000 caches in 222 countries on all continents (including Antarctica!) there are bound to be several near you.

Want to know more? Check out this video of my son on a geocache hunt[^].

So, being the true geek I am, I coded up a bunch of applications to record where we've been, where we'd like to go, and interesting waypoints along the way. Often, I'd ask the question "How far is it from here to there?"

## The Geodetic Inverse Problem

"How far?" turns out to be a fairly complicated question. For centuries, humankind has known that Earth is round - but everything else about the shape of Earth is less understood. In fact, an entire field called "geodesy[^]", or "geodetics" focuses on understanding the shape of Earth.

If we assume that Earth is a perfect sphere, the math is easy. Given the size of Earth and the latitude and longitude of two points, we can quickly calculate the "great circle[^]" distance and direction from one point to the next.

But, Earth is 'not' a perfect sphere. It's an irregular shape geodesists call the "geoid." Because precise measurements of the geoid and the corresponding distance calculations would be impractical, the geoid is often modeled as an ellipsoid - an object resembling a sphere but slightly flattened at the poles. It's still an inexact model, but it performs far better than the spherical model.

That's where the math gets hard. The Geodetic Inverse Problem - finding the distance and direction from one point to another along an ellipsoid - does not have a closed form solution. However, in 1975, Thaddeus Vincenty[^] developed an extremely accurate iterative solution for which he won the Department of Commerce Medal for Meritorious Service. If you want to sprain your brain, you can read the full publication of his solution[^].

## Building the Code

I'll start with a convenient type for encapsulating angles. We'll need this for latitudes and longitudes. It might seem to be an overkill since we could simply uses `double` instead of a custom `Angle` type, but humans think in "degrees" while computers do math in "radians." Rather than risk confusion and sprinkle conversion code throughout, I thought this would be easier:

C#
```namespace Gavaghan.Geodesy
{
public struct Angle : IComparable<angle />
{
private const double PiOver180 = Math.PI / 180.0;
private double mDegrees;

static public readonly Angle Zero = new Angle(0);
static public readonly Angle Angle180 = new Angle(180);

public Angle( double degrees )
{
mDegrees = degrees;
}

public double Degrees
{
get
{
return mDegrees;
}
set
{
mDegrees = value;
}
}

{
get
{
return mDegrees * PiOver180;
}
set
{
mDegrees = value / PiOver180;
}
}

//
// Equals, IComparable<T> implementation and
// comparison operators omitted for clarity.
//
}
}```

Some critical portions of the `Angle` type above are omitted for clarity - like all of the comparison operators. You can check out the complete implementation in the source code download[^].

Next, we have a class that encapsulates a location on an ellipsoid.

C#
```namespace Gavaghan.Geodesy
{
public struct GlobalCoordinates : IComparable<GlobalCoordinates>
{
private Angle mLatitude;
private Angle mLongitude;

public GlobalCoordinates( Angle latitude, Angle longitude )
{
mLatitude = latitude;
mLongitude = longitude;
}

public Angle Latitude
{
get
{
return mLatitude;
}
set
{
mLatitude = value;
}
}

public Angle Longitude
{
get
{
return mLongitude;
}
set
{
mLongitude = value;
}
}
}
}```

Negative latitudes are in the southern hemisphere, and negative longitudes are in the western hemisphere. Once again, the `IComparable<T>` implementation is omitted for clarity.

Now, we start getting to the fun part. The next type encapsulates the "geodetic curve."

C#
```namespace Gavaghan.Geodesy
{
public struct GeodeticCurve
{

public GeodeticCurve( double ellipsoidalDistance,
Angle azimuth, Angle reverseAzimuth )
{
mEllipsoidalDistance = ellipsoidalDistance;
mAzimuth = azimuth;
mReverseAzimuth = reverseAzimuth;
}

public double EllipsoidalDistance
{
get
{
return mEllipsoidalDistance;
}
}

public Angle Azimuth
{
get
{
return mAzimuth;
}
}

public Angle ReverseAzimuth
{
get
{
return mReverseAzimuth;
}
}
}
}```

A "geodetic curve" is the solution we're looking for. It describes how to get from one point on an ellipsoid to another. The ellipsoidal distance is the distance, in meters, between the two points along the surface of the ellipsoid. The azimuth is the direction of travel from the starting point to the ending point. The reverse azimuth, of course, is the direction back from the endpoint (which isn't necessarily a <nobr>180 degree turn on an ellipsoid).

The final input we need to Vincenty's Formula is an ellipsoid. We hold onto four parameters of an ellipsoid: the length of each semi-axis (in meters), the flattening ratio, and the inverse of the flattening ratio. Technically, we only need the length of one semi-axis and any one of the other three parameters. We'll record all four for convenience.

C#
```namespace Gavaghan.Geodesy
{
public struct Ellipsoid
{

private Ellipsoid(double semiMajor, double semiMinor,
double flattening, double inverseFlattening)
{
mSemiMajorAxis = semiMajor;
mSemiMinorAxis = semiMinor;
mFlattening = flattening;
mInverseFlattening = inverseFlattening;
}

static public readonly Ellipsoid WGS84 =
FromAAndInverseF( 6378137.0, 298.257223563 );
static public readonly Ellipsoid GRS80 =
FromAAndInverseF( 6378137.0, 298.257222101 );
static public readonly Ellipsoid GRS67 =
FromAAndInverseF( 6378160.0, 298.25 );
static public readonly Ellipsoid ANS   =
FromAAndInverseF( 6378160.0, 298.25 );
static public readonly Ellipsoid Clarke1880 =
FromAAndInverseF( 6378249.145, 293.465 );

static public Ellipsoid FromAAndInverseF
( double semiMajor, double inverseFlattening )
{
double f = 1.0 / inverseFlattening;
double b = (1.0 - f) * semiMajor;
return new Ellipsoid(semiMajor, b, f, inverseFlattening);
}

public double SemiMajorAxis
{
get
{
return mSemiMajorAxis;
}
}

public double SemiMinorAxis
{
get
{
return mSemiMinorAxis;
}
}

public double Flattening
{
get
{
return mFlattening;
}
}

public double InverseFlattening
{
get
{
return mInverseFlattening;
}
}
}
}```

Generally, you won't need to specify ellipsoid parameters. The `Ellipsoid` type has a number of static instances that define "reference ellipsoids." Reference ellipsoids represent some organization's consensus on the "best" ellipsoidal parameters to use to model Earth. Two of the most widely accepted reference ellipsoids are defined above: WGS84 (the 1984 World Geodetic System[^]) and GRS80 (the 1980 Geodetic Reference System[^]).

Finally, we have the class that actually implements Vincenty's Formula to solve the Geodetic Inverse Problem given a reference ellipsoid and two sets of global coordinates (equation numbers in comments relate directly to Vincenty's publication[^]):

C#
```namespace Gavaghan.Geodesy
{
public class GeodeticCalculator
{
private const double TwoPi = 2.0 * Math.PI;

public GeodeticCurve CalculateGeodeticCurve( Ellipsoid ellipsoid,
GlobalCoordinates start, GlobalCoordinates end )
{
// get constants
double a = ellipsoid.SemiMajorAxis;
double b = ellipsoid.SemiMinorAxis;
double f = ellipsoid.Flattening;

// calculations
double a2 = a * a;
double b2 = b * b;
double a2b2b2 = (a2 - b2) / b2;

double omega = lambda2 - lambda1;

double tanphi1 = Math.Tan(phi1);
double tanU1 = (1.0 - f) * tanphi1;
double U1 = Math.Atan(tanU1);
double sinU1 = Math.Sin(U1);
double cosU1 = Math.Cos(U1);

double tanphi2 = Math.Tan(phi2);
double tanU2 = (1.0 - f) * tanphi2;
double U2 = Math.Atan(tanU2);
double sinU2 = Math.Sin(U2);
double cosU2 = Math.Cos(U2);

double sinU1sinU2 = sinU1 * sinU2;
double cosU1sinU2 = cosU1 * sinU2;
double sinU1cosU2 = sinU1 * cosU2;
double cosU1cosU2 = cosU1 * cosU2;

// eq. 13
double lambda = omega;

// intermediates we'll need to compute 's'
double A = 0.0;
double B = 0.0;
double sigma = 0.0;
double deltasigma = 0.0;
double lambda0;
bool converged = false;

for (int i = 0; i < 10; i++)
{
lambda0 = lambda;

double sinlambda = Math.Sin(lambda);
double coslambda = Math.Cos(lambda);

// eq. 14
double sin2sigma = (cosU2 * sinlambda * cosU2 * sinlambda) +
(cosU1sinU2 - sinU1cosU2 * coslambda) *
(cosU1sinU2 - sinU1cosU2 * coslambda);

double sinsigma = Math.Sqrt(sin2sigma);

// eq. 15
double cossigma = sinU1sinU2 + (cosU1cosU2 * coslambda);

// eq. 16
sigma = Math.Atan2(sinsigma, cossigma);

// eq. 17 Careful! sin2sigma might be almost 0!
double sinalpha = (sin2sigma == 0) ? 0.0 :
cosU1cosU2 * sinlambda / sinsigma;

double alpha = Math.Asin(sinalpha);
double cosalpha = Math.Cos(alpha);
double cos2alpha = cosalpha * cosalpha;

// eq. 18 Careful! cos2alpha might be almost 0!
double cos2sigmam = cos2alpha == 0.0 ? 0.0 :
cossigma - 2 * sinU1sinU2 / cos2alpha;

double u2 = cos2alpha * a2b2b2;

double cos2sigmam2 = cos2sigmam * cos2sigmam;

// eq. 3
A = 1.0 + u2 / 16384 * (4096 + u2 *
(-768 + u2 * (320 - 175 * u2)));

// eq. 4
B = u2 / 1024 * (256 + u2 * (-128 + u2 * (74 - 47 * u2)));

// eq. 6
deltasigma = B * sinsigma * (cos2sigmam + B / 4
* (cossigma * (-1 + 2 * cos2sigmam2) - B / 6
* cos2sigmam * (-3 + 4 * sin2sigma)
* (-3 + 4 * cos2sigmam2)));

// eq. 10
double C = f / 16 * cos2alpha * (4 + f * (4 - 3 * cos2alpha));

// eq. 11 (modified)
lambda = omega + (1 - C) * f * sinalpha
* (sigma + C * sinsigma * (cos2sigmam + C
* cossigma * (-1 + 2 * cos2sigmam2)));

// see how much improvement we got
double change = Math.Abs((lambda - lambda0) / lambda);

if ((i > 1) && (change < 0.0000000000001))
{
converged = true;
break;
}
}

// eq. 19
double s = b * A * (sigma - deltasigma);
Angle alpha1;
Angle alpha2;

// didn't converge? must be N/S
if (!converged)
{
if (phi1 > phi2)
{
alpha1 = Angle.Angle180;
alpha2 = Angle.Zero;
}
else if (phi1 < phi2)
{
alpha1 = Angle.Zero;
alpha2 = Angle.Angle180;
}
else
{
alpha1 = new Angle(Double.NaN);
alpha2 = new Angle(Double.NaN);
}
}
// else, it converged, so do the math
else
{
alpha1 = new Angle();
alpha2 = new Angle();

// eq. 20
(cosU1sinU2 - sinU1cosU2 * Math.Cos(lambda)));

// eq. 21
(-sinU1cosU2 + cosU1sinU2 *
Math.Cos(lambda))) + Math.PI;

}

return new GeodeticCurve(s, alpha1, alpha2);
}
}
}```

There you have it! You have the tools you need to know the answer to the question "How far is it from here to there?"

Here's an example of using the code to calculate how far it is from the Lincoln Memorial[^] in Washington, D.C. to the Eiffel Tower[^] in Paris:

C#
```using System;
using System.Text;
using Gavaghan.Geodesy;

namespace Gavaghan.Geodesy.Example
{
public class Example
{
static void Main()
{
// instantiate the calculator
GeodeticCalculator geoCalc = new GeodeticCalculator();

// select a reference elllipsoid
Ellipsoid reference = Ellipsoid.WGS84;

// set Lincoln Memorial coordinates
GlobalCoordinates lincolnMemorial;
lincolnMemorial = new GlobalCoordinates( new Angle(38.88922),
new Angle(-77.04978) );

// set Eiffel Tower coordinates
GlobalCoordinates eiffelTower;
eiffelTower = new GlobalCoordinates( new Angle(48.85889),
new Angle(2.29583) );

// calculate the geodetic curve
GeodeticCurve geoCurve = geoCalc.CalculateGeodeticCurve(
reference, lincolnMemorial, eiffelTower );

double ellipseKilometers = geoCurve.EllipsoidalDistance / 1000.0;

Console.WriteLine("Lincoln Memorial to Eiffel Tower using WGS84");
Console.WriteLine("  Ellipsoidal Distance: {0:0.00} kilometers",
ellipseKilometers );
Console.WriteLine("  Azimuth: {0:0.00} degrees",
geoCurve.Azimuth.Degrees );
Console.WriteLine("  Reverse Azimuth: {0:0.00} degrees",
geoCurve.ReverseAzimuth.Degrees );
}
}
}```

The downloadable library also supports 3-D geodetic calculations (measurements that account for the elevation above or below the reference ellipsoid). For complete source code, documentation, and examples, download the entire C# library[^].

## Next Steps

The Vincenty implementation is just one of many useful tools for building a wide variety of applications to support the use of GPS receivers. I think mobile Web apps would be a great choice for GPS users traipsing through the outdoors. I plan to blog about[^] additional tools and ideas I'll be putting together based on my own experiences and reader feedback. Until then, buy your own GPS receiver and go play in the woods!

## Share

 United States
Mike Gavaghan opines on C# and .Net in his blog Talk Nerdy To Me[^]. He is a Microsoft Certified Professional Developer working as a C#/.Net software consultant based in Dallas, Texas.

Since 1992, he has supported clients in diverse businesses including financial services, travel, airlines, and telecom. He has consulted at both mammoth enterprises and small startups (and sees merits and problems in both).

You may also view his profile on LinkedIn[^].

 First Prev Next
 Vincenty Formula Sumal.V1-May-12 3:54 Sumal.V 1-May-12 3:54
 My vote of 5 see_seA15-Jun-11 5:11 see_seA 15-Jun-11 5:11
 Well done! see_seA15-Jun-11 5:10 see_seA 15-Jun-11 5:10
 3D direct calculation dash832-Apr-09 6:02 dash83 2-Apr-09 6:02
 PB to understand alpha1 and alpha2 of the article of Vincentry estrade30-Jul-08 10:33 estrade 30-Jul-08 10:33
 First rate job! Dave1975199528-Sep-07 5:41 Dave19751995 28-Sep-07 5:41
 Very nice work + an addition iq-man14-Sep-07 7:20 iq-man 14-Sep-07 7:20
 Thank you, your classes have helped me a lot. What I needed was the ability to calculate a new position from a starting position, a bearing and a distance. Your classes allowed me to do this quite easy. Here is the function i have added to your GeodeticCalculator class. It too uses the algoritm from Vincenty. You can certainly add it if you like, to make your litle Geodetic framework more complete. /// /// Calculates the new position when supplied with a starting position, the distance travelled in meters /// and the angle. /// /// The Ellipsoid standard to use in calculation /// The start position /// The bearing /// The distance travelled /// The calculated position public GlobalCoordinates CalculateNewPosition(Ellipsoid refEllipsoid, GlobalCoordinates start, Angle brng, double dist) { double a = refEllipsoid.SemiMajorAxis; double b = refEllipsoid.SemiMinorAxis; double f = refEllipsoid.Flattening; double s = dist; double alpha1 = brng.Radians; double sinAlpha1 = Math.Sin(alpha1), cosAlpha1 = Math.Cos(alpha1); double tanU1 = (1 - f) * Math.Tan(start.Latitude.Radians); double cosU1 = 1 / Math.Sqrt((1 + tanU1 * tanU1)), sinU1 = tanU1 * cosU1; double sigma1 = Math.Atan2(tanU1, cosAlpha1); double sinAlpha = cosU1 * sinAlpha1; double cosSqAlpha = 1 - sinAlpha * sinAlpha; double uSq = cosSqAlpha * (a * a - b * b) / (b * b); double A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq))); double B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq))); double sigma = s / (b * A), sigmaP = 2 * Math.PI; double cos2SigmaM = 0; double sinSigma = 0; double deltaSigma = 0; double cosSigma = 0; while (Math.Abs(sigma - sigmaP) > 1e-12) { cos2SigmaM = Math.Cos(2 * sigma1 + sigma); sinSigma = Math.Sin(sigma); cosSigma = Math.Cos(sigma); deltaSigma = B * sinSigma * (cos2SigmaM + B / 4 * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) - B / 6 * cos2SigmaM * (-3 + 4 * sinSigma * sinSigma) * (-3 + 4 * cos2SigmaM * cos2SigmaM))); sigmaP = sigma; sigma = s / (b * A) + deltaSigma; } double tmp = sinU1 * sinSigma - cosU1 * cosSigma * cosAlpha1; double lat2 = Math.Atan2(sinU1 * cosSigma + cosU1 * sinSigma * cosAlpha1, (1 - f) * Math.Sqrt(sinAlpha * sinAlpha + tmp * tmp)); double lambda = Math.Atan2(sinSigma * sinAlpha1, cosU1 * cosSigma - sinU1 * sinSigma * cosAlpha1); double C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha)); double L = lambda - (1 - C) * f * sinAlpha * (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM))); double revAz = Math.Atan2(sinAlpha, -tmp); // final bearing double resultLat = lat2 * (180 / Math.PI); double resultLon = start.Longitude.Degrees + L * (180 / Math.PI); return new GlobalCoordinates(new Angle(resultLat), new Angle(resultLon)); }
 error estimation Friedhelm Schuetz13-Aug-07 23:01 Friedhelm Schuetz 13-Aug-07 23:01
 Re: error estimation Mike Gavaghan14-Aug-07 1:48 Mike Gavaghan 14-Aug-07 1:48
 Looks like GFC# Sam Blackburn9-Aug-07 2:45 Sam Blackburn 9-Aug-07 2:45
 Nice one there... Paul Selormey8-Aug-07 18:12 Paul Selormey 8-Aug-07 18:12
 Last Visit: 31-Dec-99 18:00     Last Update: 18-Sep-21 18:53 Refresh 1