The following is code that I've carefully ripped off from the Interweb and
messed with. You should be able to understand /extract the bit you need. It's
a scary amount of trig. and I have absolutely no idea how it works, but seems
to be some sort of iterative process. If you feel guilty using it for free in
your project - donate some money to charity ;-)
I think that it works, but I do not have the resources to check it on a larger
scale than the Scottish tests (about 20 km distance). Perhaps others in this
list can help verify the code. I would appreciate that a lot...
=======================================================
package com.klickngo.geometry;
import org.apache.commons.lang.builder.*;
public class LatLongCoordinate {
/* These are directly accessible */
public double altitude;
public double latitudeY;
public double longitudeX;
private LatLongCoordinate() {
this(0.0, 0.0, 0.0);
}
public LatLongCoordinate(double longitudeX, double latitudeY) {
this(longitudeX, latitudeY, 0.0);
}
public LatLongCoordinate(double longitudeX, double latitudeY, double
altitude)
throws IllegalArgumentException {
if (longitudeX > 180.0 | longitudeX < -180.0) throw new
IllegalArgumentException("Longitude out of +/- 180 bounds");
if (latitudeY > 90.0 | latitudeY < -90.0) throw new
IllegalArgumentException("Latitude out of +/- 90 bounds");
if (altitude < 0.0) throw new IllegalArgumentException("Altitude below
sea level");
this.longitudeX = longitudeX;
this.latitudeY = latitudeY;
this.altitude = altitude;
}
public LatLongCoordinate(LatLongCoordinate coordinate) {
this(coordinate.longitudeX, coordinate.latitudeY, coordinate.altitude);
}
public void setCoordinate(LatLongCoordinate otherCoord) {
this.longitudeX = otherCoord.longitudeX;
this.latitudeY = otherCoord.latitudeY;
this.altitude = otherCoord.altitude;
}
/**********************************************************
*
* Did some testing in Scotland (via Memory Map).
* Two points 18 by 9 km grid squares apart, results were:
* VD = 20.125 km, Pythagarus = 20.116 km, MM = 20 km
*
* Two points 1 by 1 km grid square:-
* VD = 1.413 km, Pythagarus = 1.414 km, MM = 1.41 km
*
**********************************************************/
public double getVincentyDistance(LatLongCoordinate otherCoord) {
double retVal = 0.0D;
double a = 6378137, b = 6356752.3142, f = 1 / 298.257223563;
double l = Math.toRadians(otherCoord.longitudeX) -
Math.toRadians(this.longitudeX);
double u1 = Math.atan((1 - f) *
Math.tan(Math.toRadians(this.latitudeY)));
double u2 = Math.atan((1 - f) *
Math.tan(Math.toRadians(otherCoord.latitudeY)));
double sinU1 = Math.sin(u1), cosU1 = Math.cos(u1);
double sinU2 = Math.sin(u2), cosU2 = Math.cos(u2);
double lambda = l, lambdaP = 2 * Math.PI;
int iterLimit = 20;
double sinLambda, cosLambda;
double sinSigma = 0.0D;
double cosSigma = 0.0D;
double sigma = 0.0D;
double alpha;
double cosSqAlpha = 0.0D;
double cosToSigmaM = 0.0D;
double c;
while (Math.abs(lambda - lambdaP) > 1e-12 && --iterLimit > 0) {
sinLambda = Math.sin(lambda); cosLambda = Math.cos(lambda);
sinSigma = Math.sqrt( (cosU2 * sinLambda) * (cosU2 *
sinLambda) +
(cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) *
(cosU1 * sinU2 - sinU1 * cosU2 * cosLambda));
cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
sigma = Math.atan2(sinSigma, cosSigma);
alpha = Math.asin(cosU1 * cosU2 * sinLambda / sinSigma);
cosSqAlpha = Math.cos(alpha) * Math.cos(alpha);
cosToSigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
c = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha));
lambdaP = lambda;
lambda = l + (1 - c) * f *
Math.sin(alpha) *
(sigma + c * sinSigma *
(cosToSigmaM + c * cosSigma * (- 1 + 2 *
cosToSigmaM * cosToSigmaM)));
}
// IF no convergence
if (iterLimit == 0) {
retVal = Double.NaN;
} else {
double uSq = cosSqAlpha * (a * a - b * b) / (b * b);
double a2 = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq *
(320 - 175 * uSq)));
double b2 = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 -
47 * uSq)));
double deltaSigma = b2 * sinSigma *
(cosToSigmaM +
b2 / 4 * (cosSigma * (-1 + 2 * cosToSigmaM *
cosToSigmaM) -
b2 / 6 * cosToSigmaM *
(-3 + 4 * sinSigma * sinSigma) *
(-3 + 4 * cosToSigmaM * cosToSigmaM)));
retVal = b * a2 * (sigma - deltaSigma);
}
return retVal;
}
public String toString() {
return new ToStringBuilder(this).
append("Long.", this.longitudeX).
append("Lat.", this.latitudeY).
append("Alt.", this.altitude).
toString();
}
public boolean equals(Object obj) {
if (obj == null) { return false; }
if (obj == this) { return true; }
if (obj.getClass() != getClass()) {
return false;
}
LatLongCoordinate rhs = (LatLongCoordinate) obj;
return new EqualsBuilder()
.append(this.altitude, rhs.altitude)
.append(this.longitudeX, rhs.longitudeX)
.append(this.latitudeY, rhs.latitudeY)
.isEquals();
}
public int hashCode() {
// you pick a hard-coded, randomly chosen, non-zero, odd number
// ideally different for each class
return new HashCodeBuilder(307, 31).
append(this.altitude).
append(this.longitudeX).
append(this.latitudeY).
toHashCode();
}
/* These private accessors are only for Hibernate to use */
private double getAltitude() {
return altitude;
}
private void setAltitude(double altitude) {
this.altitude = altitude;
}
private double getLongitudeX() {
return longitudeX;
}
private void setLongitudeX(double longitudeX) {
this.longitudeX = longitudeX;
}
private double getLatitudeY() {
return latitudeY;
}
private void setLatitudeY(double latitudeY) {
this.latitudeY = latitudeY;
}
}
============================================
_______________________________________________
jts-devel mailing list
[email protected]
http://lists.refractions.net/mailman/listinfo/jts-devel
