基于经纬度计算方向，距离等

The code encapsulating a location on an Ellipsoid,

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.

reference[url]http://www.gavaghan.org/blog/2007/11/16/java-gps-receivers-and-geocaching-vincentys-formula/[/url]
[url]http://www.gavaghan.org/blog/free-source-code/geodesy-library-vincentys-formula-java/[/url]