Allow user to pass an optional ellipsoid definition to distance(), an…

…d if they do perform an ellipsoidal distance calculation rather than the default great circle.

packages/turf-distance/index.ts

@@ -1,13 +1,61 @@
 import { Point } from "geojson";
 import { getCoord } from "@turf/invariant";
-import { radiansToLength, degreesToRadians, Coord, Units } from "@turf/helpers";
+import {
+  radiansToLength,
+  degreesToRadians,
+  datums,
+  Coord,
+  Datum,
+  Units,
+  convertLength,
+} from "@turf/helpers";
+const LatLon = require("geodesy").LatLonEllipsoidal;
 
 //http://en.wikipedia.org/wiki/Haversine_formula
 //http://www.movable-type.co.uk/scripts/latlong.html
 
+/**
+ * Calculates the distance between two {@link Point|points}.
+ * If a specific datum is passed, uses a geodesic ellipsoid calculation.
+ * If no datum is passed, performs a great circle calculation.
+ *
+ * @name distance
+ * @param {Coord | Point} from origin point or coordinate
+ * @param {Coord | Point} to destination point or coordinate
+ * @param {Object} [options={}] Optional parameters
+ * @param {string} [options.units='kilometers'] can be degrees, radians, miles, or kilometers
+ * @param {Datum} [options.datum=datums.WGS84] datum listed in {@link Helpers}
+ * @returns {number} distance between the two points in chosen units
+ * @example
+ * var from = turf.point([-75.343, 39.984]);
+ * var to = turf.point([-75.534, 39.123]);
+ *
+ * var options = {units: 'miles'};
+ * var distance = turf.distance(from, to, options);
+ *
+ * //addToMap
+ * var addToMap = [from, to];
+ * from.properties.distance = distance;
+ * to.properties.distance = distance;
+ */
+function distance(
+  from: Coord | Point,
+  to: Coord | Point,
+  options: {
+    units?: Units;
+    datum?: Datum;
+  } = {}
+) {
+  if (options?.datum) {
+    return geodesic_ellipsoid_distance(from, to, options);
+  } else {
+    return great_circle_distance(from, to, options);
+  }
+}
+
 /**
  * Calculates the distance between two {@link Point|points} in degrees, radians, miles, or kilometers.
- * This uses the [Haversine formula](http://en.wikipedia.org/wiki/Haversine_formula) to account for global curvature.
+ * Performs a great circle calculation using the [Haversine formula](http://en.wikipedia.org/wiki/Haversine_formula) to account for global curvature.
  *
  * @name distance
  * @param {Coord | Point} from origin point or coordinate
@@ -27,7 +75,7 @@ import { radiansToLength, degreesToRadians, Coord, Units } from "@turf/helpers";
  * from.properties.distance = distance;
  * to.properties.distance = distance;
  */
-function distance(
+function great_circle_distance(
   from: Coord | Point,
   to: Coord | Point,
   options: {
@@ -51,4 +99,53 @@ function distance(
   );
 }
 
+/**
+ * Calculates the distance between two {@link Point|points}.
+ * Performs a geodesic ellipsoid calculation using [Vincenty's formulae](https://en.wikipedia.org/wiki/Vincenty%27s_formulae) to account for speroidal curvature.
+ *
+ * @name geodesic_ellipsoid_distance
+ * @param {Coord | Point} from origin point or coordinate
+ * @param {Coord | Point} to destination point or coordinate
+ * @param {Object} [options={}] Optional parameters
+ * @param {string} [options.units='kilometers'] can be degrees, radians, miles, or kilometers
+ * @param {Datum} [options.datum=datums.WGS84] datum listed in {@link Helpers}
+ * @returns {number} distance between the two points in chosen units
+ * @example
+ * var from = turf.point([-75.343, 39.984]);
+ * var to = turf.point([-75.534, 39.123]);
+ * var options = {units: 'miles', datum: datums.WGS84};
+ *
+ * var distance = turf.distance(from, to, options);
+ *
+ * //addToMap
+ * var addToMap = [from, to];
+ * from.properties.distance = distance;
+ * to.properties.distance = distance;
+ */
+function geodesic_ellipsoid_distance(
+  from: Coord | Point,
+  to: Coord | Point,
+  options: {
+    units?: Units;
+    datum?: Datum;
+  } = {}
+) {
+  const fromCoord = getCoord(from);
+  const toCoord = getCoord(to);
+  const fromLatLon = new LatLon(fromCoord[1], fromCoord[0]);
+  const toLatLon = new LatLon(toCoord[1], toCoord[0]);
+
+  if (options?.datum) {
+    // datum on from point sets the tone.
+    fromLatLon.datum = options.datum;
+  } else {
+    fromLatLon.datum = datums.WGS84;
+  }
+
+  const meters = fromLatLon.distanceTo(toLatLon);
+  // geodesy lib result is in meters
+  return convertLength(meters, "meters", options.units);
+}
+
+export { great_circle_distance, geodesic_ellipsoid_distance };
 export default distance;

packages/turf-distance/package.json

@@ -47,6 +47,7 @@
     "test:tape": "tsx test.js"
   },
   "devDependencies": {
+    "@types/geodesy": "1.1.3",
     "@types/tape": "*",
     "benchmark": "*",
     "load-json-file": "*",
@@ -60,6 +61,7 @@
   "dependencies": {
     "@turf/helpers": "^7.0.0-alpha.0",
     "@turf/invariant": "^7.0.0-alpha.0",
+    "geodesy": "1.1.3",
     "tslib": "^2.3.0"
   }
 }

packages/turf-distance/test.js

@@ -3,7 +3,7 @@ const path = require("path");
 const test = require("tape");
 const load = require("load-json-file");
 const write = require("write-json-file");
-const { point } = require("@turf/helpers");
+import { datums, point } from "@turf/helpers";
 const distance = require("./index").default;
 
 const directories = {
@@ -56,3 +56,41 @@ test("distance -- throws", (t) => {
   );
   t.end();
 });
+
+test("distance -- Issue #1726 line between poles", (t) => {
+  const p1 = point([-33.6, 81.1]);
+  const p2 = point([64.5, -80.8]);
+
+  let overallDistance = distance(p1, p2, {
+    units: "meters",
+    datum: datums.WGS84,
+  });
+
+  const expected = 18682436.875; // m from QGIS
+  const tolerance = 0.01; // 1 cm expressed as m
+  t.true(
+    Math.abs(overallDistance - expected) < tolerance,
+    `${overallDistance} within ${tolerance} of ${expected}`
+  );
+
+  t.end();
+});
+
+test("distance -- Issue #1726 line near equator", (t) => {
+  const p1 = point([34, 15.9]);
+  const p2 = point([21, 0.2]);
+
+  let overallDistance = distance(p1, p2, {
+    units: "meters",
+    datum: datums.WGS84,
+  });
+
+  const expected = 2248334.18; // m from QGIS
+  const tolerance = 1; // 1 cm expressed as m
+  t.true(
+    Math.abs(overallDistance - expected) < tolerance,
+    `${overallDistance} within ${tolerance} of ${expected}`
+  );
+
+  t.end();
+});

packages/turf-helpers/index.ts

@@ -14,6 +14,7 @@ import {
   Position,
   GeoJsonProperties,
 } from "geojson";
+import { Datums, LatLonEllipsoidal } from "geodesy";
 
 import { Id } from "./lib/geojson";
 export * from "./lib/geojson";
@@ -115,6 +116,18 @@ export const areaFactors: Record<AreaUnits, number> = {
   yards: 1.195990046,
 };
 
+/**
+ * Common datum and ellipsoid definitions. Re-export verbatim from geodesy.
+ *
+ * @memberof helpers
+ * @type {Object}
+ */
+export const datums: Datums = LatLonEllipsoidal.datum;
+
+// Re-export type from geodesy so clients don't need to refer directly to
+// geodesy types.
+export type { Datum } from "geodesy";
+
 /**
  * Wraps a GeoJSON {@link Geometry} in a GeoJSON {@link Feature}.
  *

packages/turf-helpers/package.json

@@ -54,6 +54,7 @@
     "test:types": "tsc --esModuleInterop --noEmit --strict types.ts"
   },
   "devDependencies": {
+    "@types/geodesy": "1.1.3",
     "@types/tape": "*",
     "benchmark": "*",
     "npm-run-all": "*",
@@ -63,6 +64,7 @@
     "typescript": "*"
   },
   "dependencies": {
+    "geodesy": "1.1.3",
     "tslib": "^2.3.0"
   }
 }

packages/turf-length/index.ts

@@ -1,6 +1,6 @@
 import { Feature, FeatureCollection, GeometryCollection } from "geojson";
 import distance from "@turf/distance";
-import { Units } from "@turf/helpers";
+import { Datum, Units } from "@turf/helpers";
 import { segmentReduce } from "@turf/meta";
 
 /**
@@ -23,6 +23,7 @@ export default function length(
   geojson: Feature<any> | FeatureCollection<any> | GeometryCollection,
   options: {
     units?: Units;
+    datum?: Datum;
   } = {}
 ): number {
   // Calculate distance from 2-vertex line segments