Add functions for spherical and Vincenty distance.

OpenGeoscience · Mar 21, 2018 · 1603fe3 · 1603fe3
1 parent 3b4070c
commit 1603fe3
Show file tree

Hide file tree

Showing 3 changed files with 216 additions and 1 deletion.
diff --git a/src/transform.js b/src/transform.js
@@ -486,4 +486,131 @@ transform.affineInverse = function (def, coords) {
   return coords;
 };
 
+/**
+ * Compute the distance on the surface on a sphere.  The sphere is the major
+ * radius of a specified ellipsoid.  Altitude is ignored.
+ *
+ * @param {geo.geoPosition} pt1 The first point.
+ * @param {geo.geoPosition} pt2 The second point.
+ * @param {string|geo.transform} [gcs] `undefined` to use the same gcs as the
+ *    ellipsoid, otherwise the gcs of the points.
+ * @param {string|geo.transform} [baseGcs='EPSG:4326'] the gcs of the
+ *    ellipsoid.
+ * @param {object} [ellipsoid=proj4.WGS84] An object with at least `a` and one
+ *    of `b`, `f`, or `rf` (1 / `f`) -- this works with  proj4 ellipsoid
+ *    definitions.
+ * @param {number} [maxIterations=100] Maximum number of iterations to use
+ *    to test convergence.
+ * @returns {number} The distance in meters (or whatever units the ellipsoid
+ *    was specified in.
+ */
+transform.sphericalDistance = function (pt1, pt2, gcs, baseGcs, ellipsoid) {
+  baseGcs = baseGcs || 'EPSG:4326';
+  ellipsoid = ellipsoid || proj4.WGS84;
+  gcs = gcs || baseGcs;
+  if (gcs !== baseGcs) {
+    var pts = transform.transformCoordinates(gcs, baseGcs, [pt1, pt2]);
+    pt1 = pts[0];
+    pt2 = pts[1];
+  }
+  // baseGcs must be in degrees or this will be wrong
+  var phi1 = pt1.y * Math.PI / 180,
+      phi2 = pt2.y * Math.PI / 180,
+      lambda = (pt2.x - pt1.x) * Math.PI / 180,
+      sinphi1 = Math.sin(phi1), cosphi1 = Math.cos(phi1),
+      sinphi2 = Math.sin(phi2), cosphi2 = Math.cos(phi2);
+  var sigma = Math.atan2(
+    Math.pow(
+      Math.pow(cosphi2 * Math.sin(lambda), 2) +
+      Math.pow(cosphi1 * sinphi2 - sinphi1 * cosphi2 * Math.cos(lambda), 2), 0.5),
+    sinphi1 * sinphi2 + cosphi1 * cosphi2 * Math.cos(lambda)
+  );
+  return ellipsoid.a * sigma;
+};
+
+/**
+ * Compute the Vincenty distance on the surface on an ellipsoid.  Altitude is
+ * ignored.
+ *
+ * @param {geo.geoPosition} pt1 The first point.
+ * @param {geo.geoPosition} pt2 The second point.
+ * @param {string|geo.transform} [gcs] `undefined` to use the same gcs as the
+ *    ellipsoid, otherwise the gcs of the points.
+ * @param {string|geo.transform} [baseGcs='EPSG:4326'] the gcs of the
+ *    ellipsoid.
+ * @param {object} [ellipsoid=proj4.WGS84] An object with at least `a` and one
+ *    of `b`, `f`, or `rf` (1 / `f`) -- this works with  proj4 ellipsoid
+ *    definitions.
+ * @param {number} [maxIterations=100] Maximum number of iterations to use
+ *    to test convergence.
+ * @returns {object} An object with `distance` in meters (or whatever units the
+ *    ellipsoid was specified in), `alpha1` and `alpha2`, the azimuths at the
+ *    two points in radians.  The result may be `undefined` if the formula
+ *    fails to converge, which can happen near antipodal points.
+ */
+transform.vincentyDistance = function (pt1, pt2, gcs, baseGcs, ellipsoid, maxIterations) {
+  baseGcs = baseGcs || 'EPSG:4326';
+  ellipsoid = ellipsoid || proj4.WGS84;
+  maxIterations = maxIterations || 100;
+  gcs = gcs || baseGcs;
+  if (gcs !== baseGcs) {
+    var pts = transform.transformCoordinates(gcs, baseGcs, [pt1, pt2]);
+    pt1 = pts[0];
+    pt2 = pts[1];
+  }
+  var a = ellipsoid.a,
+      b = ellipsoid.b || ellipsoid.a * (1.0 - (ellipsoid.f || 1.0 / ellipsoid.rf)),
+      f = ellipsoid.f || (ellipsoid.rf ? 1.0 / ellipsoid.rf : 1.0 - b / a),
+      // baseGcs must be in degrees or this will be wrong
+      phi1 = pt1.y * Math.PI / 180,
+      phi2 = pt2.y * Math.PI / 180,
+      L = (((pt2.x - pt1.x) % 360 + 360) % 360) * Math.PI / 180,
+      U1 = Math.atan((1 - f) * Math.tan(phi1)),  // reduced latitude
+      U2 = Math.atan((1 - f) * Math.tan(phi2)),
+      sinU1 = Math.sin(U1), cosU1 = Math.cos(U1),
+      sinU2 = Math.sin(U2), cosU2 = Math.cos(U2),
+      lambda = L, lastLambda = L + Math.PI * 2,
+      sinSigma, cosSigma, sigma, sinAlpha, cos2alpha, cos2sigmasubm, C,
+      u2, A, B, deltaSigma, iter;
+  if (phi1 === phi2 && !L) {
+    return {
+      distance: 0,
+      alpha1: 0,
+      alpha2: 0
+    };
+  }
+  for (iter = maxIterations; iter > 0 && Math.abs(lambda - lastLambda) > 1e-12; iter -= 1) {
+    sinSigma = Math.pow(
+        Math.pow(cosU2 * Math.sin(lambda), 2) +
+        Math.pow(cosU1 * sinU2 - sinU1 * cosU2 * Math.cos(lambda), 2), 0.5);
+    cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * Math.cos(lambda);
+    sigma = Math.atan2(sinSigma, cosSigma);
+    sinAlpha = cosU1 * cosU2 * Math.sin(lambda) / sinSigma;
+    cos2alpha = 1 - Math.pow(sinAlpha, 2);
+    // cos2alpha is zero only when phi1 and phi2 are nearly zero.  In this
+    // case, sinU1 and sinU2 are nearly zero and the the second term can be
+    // dropped
+    cos2sigmasubm = cosSigma - (cos2alpha ? 2 * sinU1 * sinU2 / cos2alpha : 0);
+    C = f / 16 * cos2alpha * (4 + f * (4 - 3 * cos2alpha));
+    lastLambda = lambda;
+    lambda = L + (1 - C) * f * sinAlpha * (sigma + C * sinSigma * (
+        cos2sigmasubm + C * cosSigma * (-1 + 2 * Math.pow(cos2sigmasubm, 2))));
+  }
+  if (!iter) { // failure to converge
+    return;
+  }
+  u2 = cos2alpha * (a * a - b * b) / (b * b);
+  A = 1 + u2 / 16384 * (4096 + u2 * (-768 + u2 * (320 - 175 * u2)));
+  B = u2 / 1024 * (256 + u2 * (-128 + u2 * (74 - 47 * u2)));
+  deltaSigma = B * sinSigma * (cos2sigmasubm + B / 4 * (
+    cosSigma * (-1 + 2 * Math.pow(cos2sigmasubm, 2)) -
+    B / 6 * cos2sigmasubm * (-3 + 4 * sinSigma * sinSigma) *
+    (-3 + 4 * Math.pow(cos2sigmasubm, 2))));
+  return {
+    distance: b * A * (sigma - deltaSigma),
+    alpha1: Math.atan2(cosU2 * Math.sin(lambda), cosU1 * sinU2 - sinU1 * cosU2 * Math.cos(lambda)),
+    alpha2: Math.atan2(cosU1 * Math.sin(lambda), -sinU1 * cosU2 + cosU1 * sinU2 * Math.cos(lambda))
+  };
+};
+
 module.exports = transform;
diff --git a/src/util/index.js b/src/util/index.js
@@ -1,4 +1,5 @@
 var $ = require('jquery');
+var proj4 = require('proj4');
 
 var throttle = require('./throttle');
 var mockVGL = require('./mockVGL');
@@ -1320,7 +1321,7 @@ var util = module.exports = {
    * Radius of the earth in meters, from the equatorial radius of SRID 4326.
    * @memberof geo.util
    */
-  radiusEarth: 6378137,
+  radiusEarth: proj4.WGS84.a,
 
   /**
    * A regular expression string that will parse a number (integer or floating

diff --git a/tests/cases/transform.js b/tests/cases/transform.js
@@ -324,4 +324,91 @@ describe('geo.transform', function () {
       expect(res[1]).toEqual({x: 0, y: -4 / 3, z: 6 / 4});
     });
   });
+
+  describe('vincentyDistance', function () {
+    it('test distance measurement', function () {
+      var result;
+      result = geo.transform.vincentyDistance(
+        {x: -71.0693514, y: 42.3541165},  // Boston
+        {x: -73.9680804, y: 40.7791472}   // New York
+      );
+      expect(result.distance).toBeCloseTo(298396.057);
+      expect(result.alpha1).toBeCloseTo(-2.180);
+      expect(result.alpha2).toBeCloseTo(-2.213);
+      expect(geo.transform.vincentyDistance(
+        {x: -73.9680804, y: 40.7791472},
+        {x: -71.0693514, y: 42.3541165}
+      ).distance).toBeCloseTo(298396.057);
+      result = geo.transform.vincentyDistance(
+        {x: -74, y: 42},
+        {x: -71, y: 42}
+      );
+      expect(result.alpha1).toBeCloseTo(1.553);
+      expect(result.alpha2).toBeCloseTo(1.588);
+      // test equal points
+      expect(geo.transform.vincentyDistance(
+        {x: -71.0693514, y: 42.3541165},
+        {x: -71.0693514, y: 42.3541165}
+      ).distance).toBe(0);
+      // test convergence failure
+      expect(geo.transform.vincentyDistance(
+        {x: 0, y: 0},
+        {x: -179.5, y: 0.5}
+      )).toBe(undefined);
+      expect(geo.transform.vincentyDistance(
+        {x: 0, y: 0},
+        {x: -179.5, y: 0.5},
+        undefined, undefined, undefined, 200
+      ).distance).toBeCloseTo(19936288.579);
+      // test near-equator distances
+      expect(geo.transform.vincentyDistance(
+        {x: 0, y: 1e-7},
+        {x: 90, y: 1e-7}
+      ).distance).toBeCloseTo(geo.util.radiusEarth * Math.PI / 2);
+      // test using a different ellipsoid
+      expect(geo.transform.vincentyDistance(
+        {x: -71.0693514, y: 42.3541165},
+        {x: -73.9680804, y: 40.7791472},
+        'EPSG:4326',
+        '+proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs',
+        {a: 6378206.4, b: 6356583.8}
+      ).distance).toBeCloseTo(298394.412);
+    });
+  });
+
+  describe('sphericalDistance', function () {
+    it('test distance measurement', function () {
+      expect(geo.transform.sphericalDistance(
+        {x: -71.0693514, y: 42.3541165},  // Boston
+        {x: -73.9680804, y: 40.7791472}   // New York
+      )).toBeCloseTo(298342.833);
+      expect(geo.transform.sphericalDistance(
+        {x: -73.9680804, y: 40.7791472},
+        {x: -71.0693514, y: 42.3541165}
+      )).toBeCloseTo(298342.833);
+      // test equal points
+      expect(geo.transform.sphericalDistance(
+        {x: -71.0693514, y: 42.3541165},
+        {x: -71.0693514, y: 42.3541165}
+      )).toBe(0);
+      // test near antipodal points
+      expect(geo.transform.sphericalDistance(
+        {x: 0, y: 0},
+        {x: -179.5, y: 0.5}
+      )).toBeCloseTo(19958794.076);
+      // test near-equator distances
+      expect(geo.transform.sphericalDistance(
+        {x: 0, y: 1e-7},
+        {x: 90, y: 1e-7}
+      )).toBeCloseTo(geo.util.radiusEarth * Math.PI / 2);
+      // test using a different ellipsoid
+      expect(geo.transform.sphericalDistance(
+        {x: -71.0693514, y: 42.3541165},
+        {x: -73.9680804, y: 40.7791472},
+        'EPSG:4326',
+        '+proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs',
+        {a: 6378206.4, b: 6356583.8}
+      )).toBeCloseTo(298340.559);
+    });
+  });
 });