Add ellipsoid geometry functions

2013-02-18 12:26:10 +01:00
parent ea6bf83c61
commit 36a3d02816
3 changed files with 571 additions and 0 deletions
--- a/src/ol/ellipsoid.js
+++ 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;
+};
--- a/src/ol/ellipsoid/wgs84.js
+++ 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);
--- a/test/spec/ol/ellipsoid.test.js
+++ 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');