diff --git a/src/ol/ellipsoid.js b/src/ol/ellipsoid.js new file mode 100644 index 0000000000..1669c7de5f --- /dev/null +++ b/src/ol/ellipsoid.js @@ -0,0 +1,170 @@ +goog.provide('ol.Ellipsoid'); + +goog.require('goog.math'); +goog.require('ol.Coordinate'); + + + +/** + * @constructor + * @param {number} a Major radius. + * @param {number} flattening Flattening. + */ +ol.Ellipsoid = function(a, flattening) { + + /** + * @type {number} + */ + this.a = a; + + /** + * @type {number} + */ + this.flattening = flattening; + + /** + * @type {number} + */ + this.b = this.a * (1 - this.flattening); + +}; + + +/** + * @param {ol.Coordinate} c1 Coordinate 1. + * @param {ol.Coordinate} c2 Coordinate 1. + * @param {number=} opt_minDeltaLambda Minimum delta lambda for convergence. + * @param {number=} opt_maxIterations Maximum iterations. + * @return {{distance: number, initialBearing: number, finalBearing: number}} + * Vincenty. + */ +ol.Ellipsoid.prototype.vincenty = + function(c1, c2, opt_minDeltaLambda, opt_maxIterations) { + var minDeltaLambda = goog.isDef(opt_minDeltaLambda) ? + opt_minDeltaLambda : 1e-12; + var maxIterations = goog.isDef(opt_maxIterations) ? + opt_maxIterations : 100; + var f = this.flattening; + var lat1 = goog.math.toRadians(c1.y); + var lat2 = goog.math.toRadians(c2.y); + var deltaLon = goog.math.toRadians(c2.x - c1.x); + var U1 = Math.atan((1 - f) * Math.tan(lat1)); + var cosU1 = Math.cos(U1); + var sinU1 = Math.sin(U1); + var U2 = Math.atan((1 - f) * Math.tan(lat2)); + var cosU2 = Math.cos(U2); + var sinU2 = Math.sin(U2); + var lambda = deltaLon; + var cosSquaredAlpha, sinAlpha; + var cosLambda, deltaLambda = Infinity, sinLambda; + var cos2SigmaM, cosSigma, sigma, sinSigma; + var i; + for (i = maxIterations; i > 0; --i) { + cosLambda = Math.cos(lambda); + sinLambda = Math.sin(lambda); + var x = cosU2 * sinLambda; + var y = cosU1 * sinU2 - sinU1 * cosU2 * cosLambda; + sinSigma = Math.sqrt(x * x + y * y); + if (sinSigma === 0) { + return { + distance: 0, + initialBearing: 0, + finalBearing: 0 + }; + } + cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda; + sigma = Math.atan2(sinSigma, cosSigma); + sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma; + cosSquaredAlpha = 1 - sinAlpha * sinAlpha; + cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSquaredAlpha; + if (isNaN(cos2SigmaM)) { + cos2SigmaM = 0; + } + var C = f / 16 * cosSquaredAlpha * (4 + f * (4 - 3 * cosSquaredAlpha)); + var lambdaPrime = deltaLon + (1 - C) * f * sinAlpha * (sigma + + C * sinSigma * (cos2SigmaM + + C * cosSigma * (2 * cos2SigmaM * cos2SigmaM - 1))); + deltaLambda = Math.abs(lambdaPrime - lambda); + lambda = lambdaPrime; + if (deltaLambda < minDeltaLambda) { + break; + } + } + if (i === 0) { + return { + distance: NaN, + finalBearing: NaN, + initialBearing: NaN + }; + } + var aSquared = this.a * this.a; + var bSquared = this.b * this.b; + var uSquared = cosSquaredAlpha * (aSquared - bSquared) / bSquared; + var A = 1 + uSquared / 16384 * + (4096 + uSquared * (uSquared * (320 - 175 * uSquared) - 768)); + var B = uSquared / 1024 * + (256 + uSquared * (uSquared * (74 - 47 * uSquared) - 128)); + var deltaSigma = B * sinSigma * (cos2SigmaM + B / 4 * + (cosSigma * (2 * cos2SigmaM * cos2SigmaM - 1) - + B / 6 * cos2SigmaM * (4 * sinSigma * sinSigma - 3) * + (4 * cos2SigmaM * cos2SigmaM - 3))); + cosLambda = Math.cos(lambda); + sinLambda = Math.sin(lambda); + var alpha1 = Math.atan2(cosU2 * sinLambda, + cosU1 * sinU2 - sinU1 * cosU2 * cosLambda); + var alpha2 = Math.atan2(cosU1 * sinLambda, + cosU1 * sinU2 * cosLambda - sinU1 * cosU2); + return { + distance: this.b * A * (sigma - deltaSigma), + initialBearing: goog.math.toDegrees(alpha1), + finalBearing: goog.math.toDegrees(alpha2) + }; +}; + + +/** + * Returns the distance from c1 to c2 using Vincenty. + * + * @param {ol.Coordinate} c1 Coordinate 1. + * @param {ol.Coordinate} c2 Coordinate 1. + * @param {number=} opt_minDeltaLambda Minimum delta lambda for convergence. + * @param {number=} opt_maxIterations Maximum iterations. + * @return {number} Vincenty distance. + */ +ol.Ellipsoid.prototype.vincentyDistance = + function(c1, c2, opt_minDeltaLambda, opt_maxIterations) { + var vincenty = this.vincenty(c1, c2, opt_minDeltaLambda, opt_maxIterations); + return vincenty.distance; +}; + + +/** + * Returns the final bearing from c1 to c2 using Vincenty. + * + * @param {ol.Coordinate} c1 Coordinate 1. + * @param {ol.Coordinate} c2 Coordinate 1. + * @param {number=} opt_minDeltaLambda Minimum delta lambda for convergence. + * @param {number=} opt_maxIterations Maximum iterations. + * @return {number} Initial bearing. + */ +ol.Ellipsoid.prototype.vincentyFinalBearing = + function(c1, c2, opt_minDeltaLambda, opt_maxIterations) { + var vincenty = this.vincenty(c1, c2, opt_minDeltaLambda, opt_maxIterations); + return vincenty.finalBearing; +}; + + +/** + * Returns the initial bearing from c1 to c2 using Vincenty. + * + * @param {ol.Coordinate} c1 Coordinate 1. + * @param {ol.Coordinate} c2 Coordinate 1. + * @param {number=} opt_minDeltaLambda Minimum delta lambda for convergence. + * @param {number=} opt_maxIterations Maximum iterations. + * @return {number} Initial bearing. + */ +ol.Ellipsoid.prototype.vincentyInitialBearing = + function(c1, c2, opt_minDeltaLambda, opt_maxIterations) { + var vincenty = this.vincenty(c1, c2, opt_minDeltaLambda, opt_maxIterations); + return vincenty.initialBearing; +}; diff --git a/src/ol/ellipsoid/wgs84.js b/src/ol/ellipsoid/wgs84.js new file mode 100644 index 0000000000..e5244b80d7 --- /dev/null +++ b/src/ol/ellipsoid/wgs84.js @@ -0,0 +1,10 @@ +goog.provide('ol.ellipsoid.WGS84'); + +goog.require('ol.Ellipsoid'); + + +/** + * @const + * @type {ol.Ellipsoid} + */ +ol.ellipsoid.WGS84 = new ol.Ellipsoid(6378137, 1 / 298.257223563); diff --git a/test/spec/ol/ellipsoid.test.js b/test/spec/ol/ellipsoid.test.js new file mode 100644 index 0000000000..cf65f8d23f --- /dev/null +++ b/test/spec/ol/ellipsoid.test.js @@ -0,0 +1,391 @@ +goog.provide('ol.test.Ellipsoid'); + + +describe('ol.Ellipsoid', function() { + + var expected = [ + { + c1: new ol.Coordinate(0, 0), + c2: new ol.Coordinate(0, 0), + vincentyFinalBearing: 0, + vincentyInitialBearing: 0, + vincentyDistance: 0 + }, + { + c1: new ol.Coordinate(0, 0), + c2: new ol.Coordinate(45, 45), + vincentyFinalBearing: 54.890773827979565, + vincentyInitialBearing: 35.41005890511814, + vincentyDistance: 6662472.718217184 + }, + { + c1: new ol.Coordinate(0, 0), + c2: new ol.Coordinate(45, -45), + vincentyFinalBearing: 125.10922617202044, + vincentyInitialBearing: 144.58994109488185, + vincentyDistance: 6662472.718217184 + }, + { + c1: new ol.Coordinate(0, 0), + c2: new ol.Coordinate(-45, -45), + vincentyFinalBearing: -125.10922617202044, + vincentyInitialBearing: -144.58994109488185, + vincentyDistance: 6662472.718217184 + }, + { + c1: new ol.Coordinate(0, 0), + c2: new ol.Coordinate(-45, 45), + vincentyFinalBearing: -54.890773827979565, + vincentyInitialBearing: -35.41005890511814, + vincentyDistance: 6662472.718217184 + }, + { + c1: new ol.Coordinate(0, 0), + c2: new ol.Coordinate(180, 90), + vincentyFinalBearing: 180, + vincentyInitialBearing: 4.296211503097554e-31, + vincentyDistance: 10001965.729311794 + }, + { + c1: new ol.Coordinate(0, 0), + c2: new ol.Coordinate(180, -90), + vincentyFinalBearing: 7.0164775638926606e-15, + vincentyInitialBearing: 180, + vincentyDistance: 10001965.729311794 + }, + { + c1: new ol.Coordinate(0, 0), + c2: new ol.Coordinate(-180, 90), + vincentyFinalBearing: -180, + vincentyInitialBearing: -4.296211503097554e-31, + vincentyDistance: 10001965.729311794 + }, + { + c1: new ol.Coordinate(0, 0), + c2: new ol.Coordinate(-180, 90), + vincentyFinalBearing: -180, + vincentyInitialBearing: -4.296211503097554e-31, + vincentyDistance: 10001965.729311794 + }, + { + c1: new ol.Coordinate(45, 45), + c2: new ol.Coordinate(45, 45), + vincentyFinalBearing: 0, + vincentyInitialBearing: 0, + vincentyDistance: 0 + }, + { + c1: new ol.Coordinate(45, 45), + c2: new ol.Coordinate(45, -45), + vincentyFinalBearing: 180, + vincentyInitialBearing: 180, + vincentyDistance: 9969888.755957305 + }, + { + c1: new ol.Coordinate(45, 45), + c2: new ol.Coordinate(-45, -45), + vincentyFinalBearing: -125.10922617202044, + vincentyInitialBearing: -125.10922617202044, + vincentyDistance: 13324945.436434371 + }, + { + c1: new ol.Coordinate(45, 45), + c2: new ol.Coordinate(-45, 45), + vincentyFinalBearing: -125.27390277185786, + vincentyInitialBearing: -54.726097228142166, + vincentyDistance: 6690232.932559058 + }, + { + c1: new ol.Coordinate(45, 45), + c2: new ol.Coordinate(180, 90), + vincentyFinalBearing: 135, + vincentyInitialBearing: 3.5023624896823797e-15, + vincentyDistance: 5017021.35133314 + }, + { + c1: new ol.Coordinate(45, 45), + c2: new ol.Coordinate(180, -90), + vincentyFinalBearing: 45.00000000000001, + vincentyInitialBearing: 180, + vincentyDistance: 14986910.107290443 + }, + { + c1: new ol.Coordinate(45, 45), + c2: new ol.Coordinate(-180, 90), + vincentyFinalBearing: 135.00000000000003, + vincentyInitialBearing: 3.5023624896823793e-15, + vincentyDistance: 5017021.35133314 + }, + { + c1: new ol.Coordinate(45, 45), + c2: new ol.Coordinate(-180, 90), + vincentyFinalBearing: 135.00000000000003, + vincentyInitialBearing: 3.5023624896823793e-15, + vincentyDistance: 5017021.35133314 + }, + { + c1: new ol.Coordinate(45, -45), + c2: new ol.Coordinate(45, -45), + vincentyFinalBearing: 0, + vincentyInitialBearing: 0, + vincentyDistance: 0 + }, + { + c1: new ol.Coordinate(45, -45), + c2: new ol.Coordinate(-45, -45), + vincentyFinalBearing: -54.726097228142166, + vincentyInitialBearing: -125.27390277185786, + vincentyDistance: 6690232.932559058 + }, + { + c1: new ol.Coordinate(45, -45), + c2: new ol.Coordinate(-45, 45), + vincentyFinalBearing: -54.890773827979565, + vincentyInitialBearing: -54.890773827979565, + vincentyDistance: 13324945.436434371 + }, + { + c1: new ol.Coordinate(45, -45), + c2: new ol.Coordinate(180, 90), + vincentyFinalBearing: 135, + vincentyInitialBearing: 3.5023624896823797e-15, + vincentyDistance: 14986910.107290443 + }, + { + c1: new ol.Coordinate(45, -45), + c2: new ol.Coordinate(180, -90), + vincentyFinalBearing: 45.00000000000001, + vincentyInitialBearing: 180, + vincentyDistance: 5017021.35133314 + }, + { + c1: new ol.Coordinate(45, -45), + c2: new ol.Coordinate(-180, 90), + vincentyFinalBearing: 135.00000000000003, + vincentyInitialBearing: 3.5023624896823793e-15, + vincentyDistance: 14986910.107290443 + }, + { + c1: new ol.Coordinate(45, -45), + c2: new ol.Coordinate(-180, 90), + vincentyFinalBearing: 135.00000000000003, + vincentyInitialBearing: 3.5023624896823793e-15, + vincentyDistance: 14986910.107290443 + }, + { + c1: new ol.Coordinate(-45, -45), + c2: new ol.Coordinate(-45, -45), + vincentyFinalBearing: 0, + vincentyInitialBearing: 0, + vincentyDistance: 0 + }, + { + c1: new ol.Coordinate(-45, -45), + c2: new ol.Coordinate(-45, 45), + vincentyFinalBearing: 0, + vincentyInitialBearing: 0, + vincentyDistance: 9969888.755957305 + }, + { + c1: new ol.Coordinate(-45, -45), + c2: new ol.Coordinate(180, 90), + vincentyFinalBearing: -135.00000000000003, + vincentyInitialBearing: -3.5023624896823793e-15, + vincentyDistance: 14986910.107290443 + }, + { + c1: new ol.Coordinate(-45, -45), + c2: new ol.Coordinate(180, -90), + vincentyFinalBearing: -44.999999999999986, + vincentyInitialBearing: -180, + vincentyDistance: 5017021.35133314 + }, + { + c1: new ol.Coordinate(-45, -45), + c2: new ol.Coordinate(-180, 90), + vincentyFinalBearing: -135, + vincentyInitialBearing: -3.5023624896823797e-15, + vincentyDistance: 14986910.107290443 + }, + { + c1: new ol.Coordinate(-45, -45), + c2: new ol.Coordinate(-180, 90), + vincentyFinalBearing: -135, + vincentyInitialBearing: -3.5023624896823797e-15, + vincentyDistance: 14986910.107290443 + }, + { + c1: new ol.Coordinate(-45, 45), + c2: new ol.Coordinate(-45, 45), + vincentyFinalBearing: 0, + vincentyInitialBearing: 0, + vincentyDistance: 0 + }, + { + c1: new ol.Coordinate(-45, 45), + c2: new ol.Coordinate(180, 90), + vincentyFinalBearing: -135.00000000000003, + vincentyInitialBearing: -3.5023624896823793e-15, + vincentyDistance: 5017021.35133314 + }, + { + c1: new ol.Coordinate(-45, 45), + c2: new ol.Coordinate(180, -90), + vincentyFinalBearing: -44.999999999999986, + vincentyInitialBearing: -180, + vincentyDistance: 14986910.107290443 + }, + { + c1: new ol.Coordinate(-45, 45), + c2: new ol.Coordinate(-180, 90), + vincentyFinalBearing: -135, + vincentyInitialBearing: -3.5023624896823797e-15, + vincentyDistance: 5017021.35133314 + }, + { + c1: new ol.Coordinate(-45, 45), + c2: new ol.Coordinate(-180, 90), + vincentyFinalBearing: -135, + vincentyInitialBearing: -3.5023624896823797e-15, + vincentyDistance: 5017021.35133314 + }, + { + c1: new ol.Coordinate(180, 90), + c2: new ol.Coordinate(180, 90), + vincentyFinalBearing: 0, + vincentyInitialBearing: 0, + vincentyDistance: 0 + }, + { + c1: new ol.Coordinate(180, 90), + c2: new ol.Coordinate(180, -90), + vincentyFinalBearing: 180, + vincentyInitialBearing: 180, + vincentyDistance: 20003931.458623584 + }, + { + c1: new ol.Coordinate(180, 90), + c2: new ol.Coordinate(-180, 90), + vincentyFinalBearing: 90, + vincentyInitialBearing: 90, + vincentyDistance: 9.565041537306137e-26 + }, + { + c1: new ol.Coordinate(180, 90), + c2: new ol.Coordinate(-180, 90), + vincentyFinalBearing: 90, + vincentyInitialBearing: 90, + vincentyDistance: 9.565041537306137e-26 + }, + { + c1: new ol.Coordinate(180, -90), + c2: new ol.Coordinate(180, -90), + vincentyFinalBearing: 0, + vincentyInitialBearing: 0, + vincentyDistance: 0 + }, + { + c1: new ol.Coordinate(180, -90), + c2: new ol.Coordinate(-180, 90), + vincentyFinalBearing: 7.0164775638926606e-15, + vincentyInitialBearing: 7.0164775638926606e-15, + vincentyDistance: 20003931.458623584 + }, + { + c1: new ol.Coordinate(180, -90), + c2: new ol.Coordinate(-180, 90), + vincentyFinalBearing: 7.0164775638926606e-15, + vincentyInitialBearing: 7.0164775638926606e-15, + vincentyDistance: 20003931.458623584 + }, + { + c1: new ol.Coordinate(-180, 90), + c2: new ol.Coordinate(-180, 90), + vincentyFinalBearing: 0, + vincentyInitialBearing: 0, + vincentyDistance: 0 + }, + { + c1: new ol.Coordinate(-180, 90), + c2: new ol.Coordinate(-180, 90), + vincentyFinalBearing: 0, + vincentyInitialBearing: 0, + vincentyDistance: 0 + }, + { + c1: new ol.Coordinate(-180, 90), + c2: new ol.Coordinate(-180, 90), + vincentyFinalBearing: 0, + vincentyInitialBearing: 0, + vincentyDistance: 0 + } + ]; + + describe('vincenty', function() { + + it('returns the same distances as Chris Veness\'s reference implementation', + function() { + var e, i, v; + for (i = 0; i < expected.length; ++i) { + e = expected[i]; + v = ol.ellipsoid.WGS84.vincenty(e.c1, e.c2, 1e-12, 100); + expect(v.distance).toRoughlyEqual(e.vincentyDistance, 1e-8); + expect(v.finalBearing).toRoughlyEqual(e.vincentyFinalBearing, 1e-9); + expect(v.initialBearing).toRoughlyEqual(e.vincentyInitialBearing, 1e-9); + } + }); + + }); + + describe('vincentyDistance', function() { + + it('returns the same distances as Chris Veness\'s reference implementation', + function() { + var e, i, vincentyDistance; + for (i = 0; i < expected.length; ++i) { + e = expected[i]; + vincentyDistance = + ol.ellipsoid.WGS84.vincentyDistance(e.c1, e.c2, 1e-12, 100); + expect(vincentyDistance).toRoughlyEqual(e.vincentyDistance, 1e-8); + } + }); + + }); + + describe('vincentyFinalBearing', function() { + + it('returns the same distances as Chris Veness\'s reference implementation', + function() { + var e, i, vincentyFinalBearing; + for (i = 0; i < expected.length; ++i) { + e = expected[i]; + vincentyFinalBearing = + ol.ellipsoid.WGS84.vincentyFinalBearing(e.c1, e.c2, 1e-12, 100); + expect(vincentyFinalBearing).toRoughlyEqual( + e.vincentyFinalBearing, 1e-9); + } + }); + + }); + + describe('vincentyInitialBearing', function() { + + it('returns the same distances as Chris Veness\'s reference implementation', + function() { + var e, i, vincentyInitialBearing; + for (i = 0; i < expected.length; ++i) { + e = expected[i]; + vincentyInitialBearing = + ol.ellipsoid.WGS84.vincentyInitialBearing(e.c1, e.c2, 1e-12, 100); + expect(vincentyInitialBearing).toRoughlyEqual( + e.vincentyInitialBearing, 1e-9); + } + }); + + }); + +}); + + +goog.require('ol.Coordinate'); +goog.require('ol.ellipsoid.WGS84');