Java, GPS Receivers, and Geocaching: Vincenty’s Formula

Save/Share Google Yahoo! Digg It Reddit del.icio.us
My Zimbio

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 Java source code from here. If you prefer C#, please see the C# version of this discussion.

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 pretense 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 for here to there?”

The Geodetic Inverse Problem

“How far?” turns out to be a fairly complicated question. For centuries, humankind has known 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 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.

Or, you can just download my Java implementation of Vincenty’s Formula.

Building the code
For those of you familiar with the C# version of the library, you’ll recall I created an Angle datatype to encapsulate the idea of angular measurements. This freed the developer from having to know when to use degrees and when to use radians. C#’s support for value datatypes made this a slick and efficient design choice.

Alas, Java does not allow for value datatypes – only reference types. Creating a reference type solely to encapsulate a simple unit conversion seemed silly to me given the garbage collection overhead it would create. Is it a premature optimization? Perhaps. At the very least, it’s another example of how Java and C# language nuances make you think differently.

So, the Angle datatype in the Java version of the geodesy library has nothing more than degree/radian conversion methods. I’ll skip the code here, but you can check it out in the full download.

The code encapsulating a location on an Ellipsoid, however, looks similar in both languages:

public class GlobalCoordinates

implements Comparable<GlobalCoordinates>, Serializable

{

private double mLatitude;

private double mLongitude;

public GlobalCoordinates(double latitude, double longitude)

{

mLatitude = latitude;

mLongitude = longitude;

}

public double getLatitude()

{

return mLatitude;

}

public void setLatitude(double latitude)

{

mLatitude = latitude;

}

public double getLongitude()

{

return mLongitude;

}

public void setLongitude(double longitude)

{

mLongitude = longitude;

}

}

Negative latitudes are in the southern hemisphere, and negative longitudes are in the western hemisphere. The Comparable<T> implementation is omitted for clarity.

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

public class GeodeticCurve implements Serializable

{

private final double mEllipsoidalDistance;

private final double mAzimuth;

private final double mReverseAzimuth;

public GeodeticCurve(double ellipsoidalDistance, double azimuth,

double reverseAzimuth)

{

mEllipsoidalDistance = ellipsoidalDistance;

mAzimuth = azimuth;

mReverseAzimuth = reverseAzimuth;

}

public double getEllipsoidalDistance()

{

return mEllipsoidalDistance;

}

public double getAzimuth()

{

return mAzimuth;

}

public double getReverseAzimuth()

{

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 180 degree turn on an ellipsoid).

The final input to Vincenty’s Formula we need 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.

public class Ellipsoid implements Serializable

{

private final double mSemiMajorAxis;

private final double mSemiMinorAxis;

private final double mFlattening;

private final double mInverseFlattening;

private Ellipsoid(double semiMajor, double semiMinor, double flattening,

double inverseFlattening)

{

mSemiMajorAxis = semiMajor;

mSemiMinorAxis = semiMinor;

mFlattening = flattening;

mInverseFlattening = inverseFlattening;

}

static public final Ellipsoid WGS84 = fromAAndInverseF(6378137.0,

298.257223563);

static public final Ellipsoid GRS80 = fromAAndInverseF(6378137.0,

298.257222101);

static public final Ellipsoid GRS67 = fromAAndInverseF(6378160.0, 298.25);

static public final Ellipsoid ANS = fromAAndInverseF(6378160.0, 298.25);

static public final 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 getSemiMajorAxis()

{

return mSemiMajorAxis;

}

public double getSemiMinorAxis()

{

return mSemiMinorAxis;

}

public double getFlattening()

{

return mFlattening;

}

public double getInverseFlattening()

{

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):

public class GeodeticCalculator

{

private final double TwoPi = 2.0 * Math.PI;

public GeodeticCurve calculateGeodeticCurve(Ellipsoid ellipsoid,

GlobalCoordinates start, GlobalCoordinates end)

{

// get constants

double a = ellipsoid.getSemiMajorAxis();

double b = ellipsoid.getSemiMinorAxis();

double f = ellipsoid.getFlattening();

// get parameters as radians

double phi1 = Angle.toRadians(start.getLatitude());

double lambda1 = Angle.toRadians(start.getLongitude());

double phi2 = Angle.toRadians(end.getLatitude());

double lambda2 = Angle.toRadians(end.getLongitude());

// 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;

boolean 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);

double alpha1;

double alpha2;

// didn’t converge? must be N/S

if (!converged)

{

if (phi1 > phi2)

{

alpha1 = 180.0;

alpha2 = 0.0;

}

else if (phi1 < phi2)

{

alpha1 = 0.0;

alpha2 = 180.0;

}

else

{

alpha1 = Double.NaN;

alpha2 = Double.NaN;

}

}

// else, it converged, so do the math

else

{

double radians;

// eq. 20

radians = Math.atan2(cosU2 * Math.sin(lambda),

(cosU1sinU2 – sinU1cosU2 * Math.cos(lambda)));

if (radians < 0.0) radians += TwoPi;

alpha1 = Angle.toDegrees(radians);

// eq. 21

radians = Math.atan2(cosU1 * Math.sin(lambda),

(-sinU1cosU2 + cosU1sinU2 * Math.cos(lambda))) + Math.PI;

if (radians < 0.0) radians += TwoPi;

alpha2 = Angle.toDegrees(radians);

}

if (alpha1 >= 360.0) alpha1 -= 360.0;

if (alpha2 >= 360.0) alpha2 -= 360.0;

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:

public class Example

{

static public void main(String[] args)

{

// instantiate the calculator

GeodeticCalculator geoCalc = new GeodeticCalculator();

// select a reference elllipsoid

Ellipsoid reference = Ellipsoid.WGS84;

// set Lincoln Memorial coordinates

GlobalCoordinates lincolnMemorial;

lincolnMemorial = new GlobalCoordinates(38.88922, -77.04978);

// set Eiffel Tower coordinates

GlobalCoordinates eiffelTower;

eiffelTower = new GlobalCoordinates(48.85889, 2.29583);

// calculate the geodetic curve

GeodeticCurve geoCurve = geoCalc.calculateGeodeticCurve(

reference, lincolnMemorial, eiffelTower

);

double ellipseKilometers = geoCurve.getEllipsoidalDistance() / 1000.0;

double ellipseMiles = ellipseKilometers * 0.621371192;

System.out.println(“Lincoln Memorial to Eiffel Tower using WGS84”);

System.out.printf(

” Ellipsoidal Distance: %1.2f kilometers (%1.2f miles)\n”,

ellipseKilometers, ellipseMiles

);

System.out.printf(” Azimuth: %1.2f degrees\n”,

geoCurve.getAzimuth()

);

System.out.printf(” Reverse Azimuth: %1.2f degrees\n”,

geoCurve.getReverseAzimuth()

);

}

}

The 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 Java library.

Save/Share Google Yahoo! Add to Technorati Favorites Digg It Reddit
del.icio.us My Zimbio

13 Responses to “Java, GPS Receivers, and Geocaching: Vincenty’s Formula”

  1. Javier Says:

    Fantastic article about the application of Vincenty’s formula. It can be very useful.

  2. paul Says:

    Hi,
    Yours is the only c# implementation of Vincenty’s direct formula that I’ve come across on the net and I’ve a little problem that you might be able to help me with. Before I found this implementation I recoded the NGS downloadable fortan into c# myself and while both my and your implementations both work to give the same results, I’m finding them wrong and am wondering if you might know what the problem could be. I’m a PhD student in Ireland, Europe and if I enter the following simple values into either my or your implementation I get the same results to 10 decimal places, these are wrong though. However, if I enter these values into the NGS online calculator I get the correct results – http://www.ngs.noaa.gov/cgi-bin/Inv_Fwd/forward2.prl . Our implementations are producing an incorrect Longitude and I can’t seem to isolate the problem.
    Values:
    Latitude: 53.387125N
    Longitude: 6.58569W
    Azimuth: 218.512579
    Distance: 59.352571

    Any help you might be able to offer would be greatly appreciated.
    Thanks in advance.

    Paul

  3. Paul Says:

    My complete apologies, there is nothing wrong with either your or my code implementation. I’m just a brainless idiot who has spent far too long looking for a big problem when a minus sign should have been obvious. Sorry.

    Paul

  4. Peter Says:

    Can I use this code to add a certain distance to a coordinate? So calculating a coordinate that is e.g. 100miles north of a given coordinate?

  5. Jacky Says:

    This API is awesome. Very useful for me. Can I include this API in my system? Does any license I need follow?

  6. Michael Says:

    I tried the GeodeticCalculator for use with Nasa’s open source WorldWindJava and after a quick set of tests that converted WorldWind Java’s LatLon, Position and Angle to the correct input required by the GeodeticCalculator the results were very accurate.

  7. rick Says:

    I used the Geodetic Calculator to determine the azimuth given 2 points at the same latitude using the WGS84 ellipsoid:
    GlobalPosition(45, 100, 0.0);
    GlobalPosition(45, 140, 0.0);

    However, the azimuth returned was not 90 degrees as I expected, but about 75 degrees. As I change the longitude, I get different angles for the azimuth. Maybe I don’t understand this correctly, but the answer seems counter-intuitive. Could someone explain?

    thanks

    rick

  8. Mike Gavaghan Says:

    Hi Rick,

    The reason is because the shortest path between two points at the same latitude doesn’t follow the latitudinal line (unless you’re on the equator). If you’re not on the equator, the shortest path will bow out toward the nearest pole.

    Try GlobalPosition(45, 100, 0.0) and GlobalPosition(45, 101, 0.0). Because these points are very close, the azimuth is 89.65 (very nearly 90).

    Now, try GlobalPosition(45, 100, 0.0) and GlobalPosition(45, -75, 0.0). These two points are on nearly opposite sides. The azimuth is 3.53. In order words, the shortest path takes you right by the north pole.

    Make sense?

    –Mike

  9. rick Says:

    That makes perfect sense – thanks.

    rick

  10. Noel Grandin Says:

    Wow, thanks. This is the only code I’ve found so far (and I’ve tried about 5 today) that matches up nicely with what Google Earth thinks is the bearing between 2 points.
    Saved me a ton of headache-inducing mathematics 🙂

  11. Ukho Says:

    Wow thanks! This is awesome. I’m overwhelmed by this source code! I am a grade 11 student from South Africa. I’m doing an I.T PAT using JAVA and I chose to do a geocaching game with geometry used in real life situations. So reading this explanation made me have a better understanding, even though I couldn’t comprehend the whole source code. Thanks to the compiler! I hope I’ll manage with my small mini coding game. 🙂

  12. anonymous Says:

    Hi,

    can this algorithm be used for highly fluctuating GPS data, i mean does it provides accuracy too, as we know that GPS data fluctuates very heavily.

    Thanks

  13. anonymous Says:

    Actually this code gives us the displacement between 2 coordinates, can you help me with the calculation of driving distance between 2 points.

Leave a Reply