Allow optional ellipsoid definition to be passed to distance allowing for calculation types other than great circle

packages/turf-distance/index.ts

@@ -1,25 +1,33 @@
 import { Point } from "geojson";
 import { getCoord } from "@turf/invariant";
-import { radiansToLength, degreesToRadians, Coord, Units } from "@turf/helpers";
-
+import type { Coord, Datum, Units } from "@turf/helpers";
+import {
+  datums,
+  radiansToLength,
+  degreesToRadians,
+  convertLength,
+} from "@turf/helpers";
+import LatLon from "geodesy/latlon-ellipsoidal-vincenty.js";
 //http://en.wikipedia.org/wiki/Haversine_formula
 //http://www.movable-type.co.uk/scripts/latlong.html
 
 /**
- * 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.
+ * 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
- * @returns {number} distance between the two points
+ * @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 options = {units: 'miles'};
  * var distance = turf.distance(from, to, options);
  *
  * //addToMap
@@ -28,6 +36,46 @@ import { radiansToLength, degreesToRadians, Coord, Units } from "@turf/helpers";
  * to.properties.distance = distance;
  */
 function distance(
+  from: Coord | Point,
+  to: Coord | Point,
+  options: {
+    units?: Units;
+    datum?: Datum;
+  } = {}
+) {
+  if (options?.datum) {
+    return geodesicEllipsoidDistance(from, to, {
+      datum: options.datum, // Must pass a datum option
+      ...options,
+    });
+  } else {
+    return greatCircleDistance(from, to, options);
+  }
+}
+
+/**
+ * Calculates the distance between two {@link Point|points} in degrees, radians, miles, or kilometers.
+ * Performs a great circle calculation using the [Haversine formula](http://en.wikipedia.org/wiki/Haversine_formula) to account for global curvature.
+ *
+ * @name greatCircleDistance
+ * @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
+ * @returns {number} distance between the two points
+ * @example
+ * var from = turf.point([-75.343, 39.984]);
+ * var to = turf.point([-75.534, 39.123]);
+ * var options = {units: 'miles'};
+ *
+ * var distance = turf.greatCircleDistance(from, to, options);
+ *
+ * //addToMap
+ * var addToMap = [from, to];
+ * from.properties.distance = distance;
+ * to.properties.distance = distance;
+ */
+function greatCircleDistance(
   from: Coord | Point,
   to: Coord | Point,
   options: {
@@ -51,4 +99,51 @@ 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 geodesicEllipsoidDistance
+ * @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.geodesicEllipsoidDistance(from, to, options);
+ *
+ * //addToMap
+ * var addToMap = [from, to];
+ * from.properties.distance = distance;
+ * to.properties.distance = distance;
+ */
+function geodesicEllipsoidDistance(
+  from: Coord | Point,
+  to: Coord | Point,
+  options: {
+    units?: Units;
+    datum: Datum;
+  } = { datum: datums.WGS84 }
+) {
+  const fromCoord = getCoord(from);
+  const toCoord = getCoord(to);
+
+  const datum = options.datum;
+
+  const fromLatLon = new LatLon(fromCoord[1], fromCoord[0]);
+  // datum on from point sets the tone.
+  fromLatLon.datum = datum;
+  const toLatLon = new LatLon(toCoord[1], toCoord[0]);
+
+  const meters = fromLatLon.distanceTo(toLatLon);
+  // geodesy lib result is in meters
+  return convertLength(meters, "meters", options.units);
+}
+
+export { greatCircleDistance, geodesicEllipsoidDistance };
 export default distance;

packages/turf-distance/package.json

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

packages/turf-distance/test.js

@@ -3,8 +3,9 @@ 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");
+const { datums, point } = require("@turf/helpers");
 const distance = require("./index").default;
+const { geodesicEllipsoidDistance } = require("./index");
 
 const directories = {
   in: path.join(__dirname, "test", "in") + path.sep,
@@ -56,3 +57,59 @@ 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 = 0.01; // 1 cm expressed as m
+  t.true(
+    Math.abs(overallDistance - expected) < tolerance,
+    `${overallDistance} within ${tolerance} of ${expected}`
+  );
+
+  t.end();
+});
+
+test("distance -- calling ellipsoid distance with empty options", (t) => {
+  // Make sure calling geodesicEllipsoidDistance directly without a datum
+  // option defaults to WGS84.
+  const p1 = point([-33.6, 81.1]);
+  const p2 = point([64.5, -80.8]);
+
+  let overallDistance = geodesicEllipsoidDistance(p1, p2);
+
+  const expected = 18682.436875; // m from QGIS as km
+  const tolerance = 0.00001; // 1 cm expressed as km
+  t.true(
+    Math.abs(overallDistance - expected) < tolerance,
+    `${overallDistance} within ${tolerance} of ${expected}`
+  );
+
+  t.end();
+});

packages/turf-helpers/index.ts

@@ -14,6 +14,8 @@ import {
   Position,
   GeoJsonProperties,
 } from "geojson";
+import type { Datum } from "geodesy";
+import LatLonEllipsoidalVincenty from "geodesy/latlon-ellipsoidal-vincenty.js";
 
 import { Id } from "./lib/geojson";
 export * from "./lib/geojson";
@@ -115,6 +117,19 @@ export const areaFactors: Record<AreaUnits, number> = {
   yards: 1.195990046,
 };
 
+/**
+ * Common datum and ellipsoid definitions. Re-export verbatim from geodesy.
+ *
+ * @memberof helpers
+ * @type {Object}
+ */
+const datums = LatLonEllipsoidalVincenty.datums;
+export { datums };
+
+// Re-export type from geodesy so clients don't need to refer directly to
+// geodesy types.
+export type { Datum };
+
 /**
  * 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": "^2.2.3",
     "@types/tape": "*",
     "benchmark": "*",
     "npm-run-all": "*",
@@ -62,6 +63,7 @@
     "typescript": "*"
   },
   "dependencies": {
+    "geodesy": "^2.3.0",
     "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 type { 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