refactor LngLat distance calculations (mapbox#9202)

src/geo/lng_lat.js

@@ -3,6 +3,13 @@
 import {wrap} from '../util/util';
 import LngLatBounds from './lng_lat_bounds';
 
+/*
+* Approximate radius of the earth in meters.
+* Uses the WGS-84 approximation. The radius at the equator is ~6378137 and at the poles is ~6356752. https://en.wikipedia.org/wiki/World_Geodetic_System#WGS84
+* 6371008.8 is one published "average radius" see https://en.wikipedia.org/wiki/Earth_radius#Mean_radius, or ftp://athena.fsv.cvut.cz/ZFG/grs80-Moritz.pdf p.4
+*/
+export const earthRadius = 6371008.8;
+
 /**
  * A `LngLat` object represents a given longitude and latitude coordinate, measured in degrees.
  *
@@ -73,6 +80,27 @@ class LngLat {
         return `LngLat(${this.lng}, ${this.lat})`;
     }
 
+    /**
+     * Returns the approximate distance between a pair of coordinates in meters
+     * Uses the Haversine Formula (from R.W. Sinnott, "Virtues of the Haversine", Sky and Telescope, vol. 68, no. 2, 1984, p. 159)
+     *
+     * @param {LngLat} lngLat coordinates to compute the distance to
+     * @returns {number} Distance in meters between the two coordinates.
+     * @example
+     * var new_york = new mapboxgl.LngLat(-74.0060, 40.7128);
+     * var los_angeles = new mapboxgl.LngLat(-118.2437, 34.0522);
+     * new_york.distanceTo(los_angeles); // = 3935751.690893987, "true distance" using a non-spherical approximation is ~3966km
+     */
+    distanceTo(lngLat: LngLat) {
+        const rad = Math.PI / 180;
+        const lat1 = this.lat * rad;
+        const lat2 = lngLat.lat * rad;
+        const a = Math.sin(lat1) * Math.sin(lat2) + Math.cos(lat1) * Math.cos(lat2) * Math.cos((lngLat.lng - this.lng) * rad);
+
+        const maxMeters = earthRadius * Math.acos(Math.min(a, 1));
+        return maxMeters;
+    }
+
     /**
      * Returns a `LngLatBounds` from the coordinates extended by a given `radius`. The returned `LngLatBounds` completely contains the `radius`.
      *

src/geo/mercator_coordinate.js

@@ -1,18 +1,18 @@
 // @flow
 
-import LngLat from '../geo/lng_lat';
+import LngLat, {earthRadius} from '../geo/lng_lat';
 import type {LngLatLike} from '../geo/lng_lat';
 
 /*
- * The circumference of the world in meters at the equator.
+ * The average circumference of the world in meters.
  */
-const circumferenceAtEquator = 2 * Math.PI * 6378137;
+const earthCircumfrence = 2 * Math.PI * earthRadius; // meters
 
 /*
- * The circumference of the world in meters at the given latitude.
+ * The circumference at a line of latitude in meters.
  */
 function circumferenceAtLatitude(latitude: number) {
-    return circumferenceAtEquator * Math.cos(latitude * Math.PI / 180);
+    return earthCircumfrence * Math.cos(latitude * Math.PI / 180);
 }
 
 export function mercatorXfromLng(lng: number) {
@@ -142,7 +142,7 @@ class MercatorCoordinate {
      */
     meterInMercatorCoordinateUnits() {
         // 1 meter / circumference at equator in meters * Mercator projection scale factor at this latitude
-        return 1 / circumferenceAtEquator * mercatorScale(latFromMercatorY(this.y));
+        return 1 / earthCircumfrence * mercatorScale(latFromMercatorY(this.y));
     }
 
 }

src/ui/control/scale_control.js

@@ -92,7 +92,9 @@ function updateScale(map, container, options) {
     const maxWidth = options && options.maxWidth || 100;
 
     const y = map._container.clientHeight / 2;
-    const maxMeters = getDistance(map.unproject([0, y]), map.unproject([maxWidth, y]));
+    const left = map.unproject([0, y]);
+    const right = map.unproject([maxWidth, y]);
+    const maxMeters = left.distanceTo(right);
     // The real distance corresponding to 100px scale length is rounded off to
     // near pretty number and the scale length for the same is found out.
     // Default unit of the scale is based on User's locale.
@@ -121,21 +123,6 @@ function setScale(container, maxWidth, maxDistance, unit) {
     container.innerHTML = distance + unit;
 }
 
-function getDistance(latlng1, latlng2) {
-    // Uses spherical law of cosines approximation.
-    const R = 6371000;
-
-    const rad = Math.PI / 180,
-        lat1 = latlng1.lat * rad,
-        lat2 = latlng2.lat * rad,
-        a = Math.sin(lat1) * Math.sin(lat2) +
-          Math.cos(lat1) * Math.cos(lat2) * Math.cos((latlng2.lng - latlng1.lng) * rad);
-
-    const maxMeters = R * Math.acos(Math.min(a, 1));
-    return maxMeters;
-
-}
-
 function getDecimalRoundNum(d) {
     const multiplier = Math.pow(10, Math.ceil(-Math.log(d) / Math.LN10));
     return Math.round(d * multiplier) / multiplier;

test/integration/render-tests/zoomed-fill/negative-zoom/style.json

@@ -2,6 +2,7 @@
   "version": 8,
   "metadata": {
     "test": {
+      "allowed": 0.0002,
       "width": 256,
       "height": 256
     }

test/unit/geo/lng_lat.test.js

@@ -56,6 +56,33 @@ test('LngLat', (t) => {
         t.end();
     });
 
+    t.test('#distanceTo', (t) => {
+        const newYork = new LngLat(-74.0060, 40.7128);
+        const losAngeles = new LngLat(-118.2437, 34.0522);
+        const d = newYork.distanceTo(losAngeles); // 3935751.690893987, "true distance" is 3966km
+        t.ok(d > 3935750, "New York should be more than 3935750m from Los Angeles");
+        t.ok(d < 3935752, "New York should be less than 3935752m from Los Angeles");
+        t.end();
+    });
+
+    t.test('#distanceTo to pole', (t) => {
+        const newYork = new LngLat(-74.0060, 40.7128);
+        const northPole = new LngLat(-135, 90);
+        const d = newYork.distanceTo(northPole); // 5480494.158486183 , "true distance" is 5499km
+        t.ok(d > 5480493, "New York should be more than 5480493m from the North Pole");
+        t.ok(d < 5480495, "New York should be less than 5480495m from the North Pole");
+        t.end();
+    });
+
+    t.test('#distanceTo to Null Island', (t) => {
+        const newYork = new LngLat(-74.0060, 40.7128);
+        const nullIsland = new LngLat(0, 0);
+        const d = newYork.distanceTo(nullIsland); // 8667080.125666846 , "true distance" is 8661km
+        t.ok(d > 8667079, "New York should be more than 8667079m from Null Island");
+        t.ok(d < 8667081, "New York should be less than 8667081m from Null Island");
+        t.end();
+    });
+
     t.test('#toBounds', (t) => {
         t.deepEqual(new LngLat(0, 0).toBounds(10).toArray(), [[-0.00008983152770714982, -0.00008983152770714982], [0.00008983152770714982, 0.00008983152770714982]]);
         t.deepEqual(new LngLat(-73.9749, 40.7736).toBounds(10).toArray(), [[-73.97501862141328, 40.77351016847229], [-73.97478137858673, 40.77368983152771]]);

test/unit/geo/mercator_coordinate.test.js

@@ -35,7 +35,7 @@ test('LngLat', (t) => {
 
     t.test('#meterInMercatorCoordinateUnits', (t) => {
         const nullIsland = new LngLat(0, 0);
-        t.equal(MercatorCoordinate.fromLngLat(nullIsland).meterInMercatorCoordinateUnits(), 2.495320233665337e-8, 'length of 1 meter in MercatorCoordinate units at the equator');
+        t.equal(MercatorCoordinate.fromLngLat(nullIsland).meterInMercatorCoordinateUnits(), 2.4981121214570498e-8, 'length of 1 meter in MercatorCoordinate units at the equator');
         t.end();
     });