/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* Vincenty Direct and Inverse Solution of Geodesics on the Ellipsoid (c) Chris Veness 2002-2021 */
/* MIT Licence */
/* www.movable-type.co.uk/scripts/latlong-vincenty.html */
/* www.movable-type.co.uk/scripts/geodesy-library.html#latlon-ellipsoidal-vincenty */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
import LatLonEllipsoidal, { Dms } from './latlon-ellipsoidal.js';
const π = Math.PI;
const ε = Number.EPSILON;
/**
* Distances & bearings between points, and destination points given start points & initial bearings,
* calculated on an ellipsoidal earth model using ‘direct and inverse solutions of geodesics on the
* ellipsoid’ devised by Thaddeus Vincenty.
*
* From: T Vincenty, "Direct and Inverse Solutions of Geodesics on the Ellipsoid with application of
* nested equations", Survey Review, vol XXIII no 176, 1975. www.ngs.noaa.gov/PUBS_LIB/inverse.pdf.
*
* @module latlon-ellipsoidal-vincenty
*/
/* LatLonEllipsoidal_Vincenty - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/**
* Extends LatLonEllipsoidal with methods for calculating distances and bearings between points, and
* destination points given distances and initial bearings, accurate to within 0.5mm distance,
* 0.000015″ bearing.
*
* By default, these calculations are made on a WGS-84 ellipsoid. For geodesic calculations on other
* ellipsoids, monkey-patch the LatLon point by setting the datum of ‘this’ point to make it appear
* as a LatLonEllipsoidal_Datum or LatLonEllipsoidal_ReferenceFrame point: e.g.
*
* import LatLon, { Dms } from '../latlon-ellipsoidal-vincenty.js';
* import { datums } from '../latlon-ellipsoidal-datum.js';
* const le = new LatLon(50.065716, -5.713824); // in OSGB-36
* const jog = new LatLon(58.644399, -3.068521); // in OSGB-36
* le.datum = datums.OSGB36; // source point determines ellipsoid to use
* const d = le.distanceTo(jog); // = 969982.014; 27.848m more than on WGS-84 ellipsoid
*
* @extends LatLonEllipsoidal
*/
class LatLonEllipsoidal_Vincenty extends LatLonEllipsoidal {
/**
* Returns the distance between ‘this’ point and destination point along a geodesic on the
* surface of the ellipsoid, using Vincenty inverse solution.
*
* @param {LatLon} point - Latitude/longitude of destination point.
* @returns {number} Distance in metres between points or NaN if failed to converge.
*
* @example
* const p1 = new LatLon(50.06632, -5.71475);
* const p2 = new LatLon(58.64402, -3.07009);
* const d = p1.distanceTo(p2); // 969,954.166 m
*/
distanceTo(point) {
try {
const dist = this.inverse(point).distance;
return Number(dist.toFixed(3)); // round to 1mm precision
} catch (e) {
if (e instanceof EvalError) return NaN; // λ > π or failed to converge
throw e;
}
}
/**
* Returns the initial bearing to travel along a geodesic from ‘this’ point to the given point,
* using Vincenty inverse solution.
*
* @param {LatLon} point - Latitude/longitude of destination point.
* @returns {number} Initial bearing in degrees from north (0°..360°) or NaN if failed to converge.
*
* @example
* const p1 = new LatLon(50.06632, -5.71475);
* const p2 = new LatLon(58.64402, -3.07009);
* const b1 = p1.initialBearingTo(p2); // 9.1419°
*/
initialBearingTo(point) {
try {
const brng = this.inverse(point).initialBearing;
return Number(brng.toFixed(7)); // round to 0.001″ precision
} catch (e) {
if (e instanceof EvalError) return NaN; // λ > π or failed to converge
throw e;
}
}
/**
* Returns the final bearing having travelled along a geodesic from ‘this’ point to the given
* point, using Vincenty inverse solution.
*
* @param {LatLon} point - Latitude/longitude of destination point.
* @returns {number} Final bearing in degrees from north (0°..360°) or NaN if failed to converge.
*
* @example
* const p1 = new LatLon(50.06632, -5.71475);
* const p2 = new LatLon(58.64402, -3.07009);
* const b2 = p1.finalBearingTo(p2); // 11.2972°
*/
finalBearingTo(point) {
try {
const brng = this.inverse(point).finalBearing;
return Number(brng.toFixed(7)); // round to 0.001″ precision
} catch (e) {
if (e instanceof EvalError) return NaN; // λ > π or failed to converge
throw e;
}
}
/**
* Returns the destination point having travelled the given distance along a geodesic given by
* initial bearing from ‘this’ point, using Vincenty direct solution.
*
* @param {number} distance - Distance travelled along the geodesic in metres.
* @param {number} initialBearing - Initial bearing in degrees from north.
* @returns {LatLon} Destination point.
*
* @example
* const p1 = new LatLon(-37.95103, 144.42487);
* const p2 = p1.destinationPoint(54972.271, 306.86816); // 37.6528°S, 143.9265°E
*/
destinationPoint(distance, initialBearing) {
return this.direct(Number(distance), Number(initialBearing)).point;
}
/**
* Returns the final bearing having travelled along a geodesic given by initial bearing for a
* given distance from ‘this’ point, using Vincenty direct solution.
* TODO: arg order? (this is consistent with destinationPoint, but perhaps less intuitive)
*
* @param {number} distance - Distance travelled along the geodesic in metres.
* @param {LatLon} initialBearing - Initial bearing in degrees from north.
* @returns {number} Final bearing in degrees from north (0°..360°).
*
* @example
* const p1 = new LatLon(-37.95103, 144.42487);
* const b2 = p1.finalBearingOn(54972.271, 306.86816); // 307.1736°
*/
finalBearingOn(distance, initialBearing) {
const brng = this.direct(Number(distance), Number(initialBearing)).finalBearing;
return Number(brng.toFixed(7)); // round to 0.001″ precision
}
/**
* Returns the point at given fraction between ‘this’ point and given point.
*
* @param {LatLon} point - Latitude/longitude of destination point.
* @param {number} fraction - Fraction between the two points (0 = this point, 1 = specified point).
* @returns {LatLon} Intermediate point between this point and destination point.
*
* @example
* const p1 = new LatLon(50.06632, -5.71475);
* const p2 = new LatLon(58.64402, -3.07009);
* const pInt = p1.intermediatePointTo(p2, 0.5); // 54.3639°N, 004.5304°W
*/
intermediatePointTo(point, fraction) {
if (fraction == 0) return this;
const inverse = this.inverse(point);
const dist = inverse.distance;
const brng = inverse.initialBearing;
return isNaN(brng) ? this : this.destinationPoint(dist*fraction, brng);
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/**
* Vincenty direct calculation.
*
* Ellipsoid parameters are taken from datum of 'this' point. Height is ignored.
*
* @private
* @param {number} distance - Distance along bearing in metres.
* @param {number} initialBearing - Initial bearing in degrees from north.
* @returns (Object} Object including point (destination point), finalBearing.
* @throws {RangeError} Point must be on surface of ellipsoid.
* @throws {EvalError} Formula failed to converge.
*/
direct(distance, initialBearing) {
if (isNaN(distance)) throw new TypeError(`invalid distance ${distance}`);
if (distance == 0) return { point: this, finalBearing: NaN, iterations: 0 };
if (isNaN(initialBearing)) throw new TypeError(`invalid bearing ${initialBearing}`);
if (this.height != 0) throw new RangeError('point must be on the surface of the ellipsoid');
const φ1 = this.lat.toRadians(), λ1 = this.lon.toRadians();
const α1 = Number(initialBearing).toRadians();
const s = Number(distance);
// allow alternative ellipsoid to be specified
const ellipsoid = this.datum ? this.datum.ellipsoid : LatLonEllipsoidal.ellipsoids.WGS84;
const { a, b, f } = ellipsoid;
const sinα1 = Math.sin(α1);
const cosα1 = Math.cos(α1);
const tanU1 = (1-f) * Math.tan(φ1), cosU1 = 1 / Math.sqrt((1 + tanU1*tanU1)), sinU1 = tanU1 * cosU1;
const σ1 = Math.atan2(tanU1, cosα1); // σ1 = angular distance on the sphere from the equator to P1
const sinα = cosU1 * sinα1; // α = azimuth of the geodesic at the equator
const cosSqα = 1 - sinα*sinα;
const uSq = cosSqα * (a*a - b*b) / (b*b);
const A = (((uSq*-175 + 320)*uSq - 768)*uSq + 4096)*uSq/16384 + 1;
const B = (((uSq*-47 + 74)*uSq - 128)*uSq + 256)*uSq/1024;
let σ = s / (b*A), sinσ = null, cosσ = null; // σ = angular distance P₁ P₂ on the sphere
let cos2σₘ = null; // σₘ = angular distance on the sphere from the equator to the midpoint of the line
let σʹ = null, iterations = 0;
do {
cos2σₘ = Math.cos(2*σ1 + σ);
sinσ = Math.sin(σ);
cosσ = Math.cos(σ);
const Δσ = B*sinσ*(cos2σₘ+B/4*(cosσ*(-1+2*cos2σₘ*cos2σₘ)-
B/6*cos2σₘ*(-3+4*sinσ*sinσ)*(-3+4*cos2σₘ*cos2σₘ)));
σʹ = σ;
σ = s / (b*A) + Δσ;
} while (Math.abs(σ-σʹ) > 1e-12 && ++iterations<100);
if (iterations >= 100) throw new EvalError('Vincenty formula failed to converge'); // not possible?
const x = sinU1*sinσ - cosU1*cosσ*cosα1;
const φ2 = Math.atan2(sinU1*cosσ + cosU1*sinσ*cosα1, (1-f)*Math.sqrt(sinα*sinα + x*x));
const λ = Math.atan2(sinσ*sinα1, cosU1*cosσ - sinU1*sinσ*cosα1);
const C = f/16*cosSqα*(4+f*(4-3*cosSqα));
const L = λ - (1-C) * f * sinα * (σ + C*sinσ*(cos2σₘ+C*cosσ*(-1+2*cos2σₘ*cos2σₘ)));
const λ2 = λ1 + L;
const α2 = Math.atan2(sinα, -x);
const destinationPoint = new LatLonEllipsoidal_Vincenty(φ2.toDegrees(), λ2.toDegrees(), 0, this.datum);
return {
point: destinationPoint,
finalBearing: Dms.wrap360(α2.toDegrees()),
iterations: iterations,
};
}
/**
* Vincenty inverse calculation.
*
* Ellipsoid parameters are taken from datum of 'this' point. Height is ignored.
*
* @private
* @param {LatLon} point - Latitude/longitude of destination point.
* @returns {Object} Object including distance, initialBearing, finalBearing.
* @throws {TypeError} Invalid point.
* @throws {RangeError} Points must be on surface of ellipsoid.
* @throws {EvalError} Formula failed to converge.
*/
inverse(point) {
if (!(point instanceof LatLonEllipsoidal)) throw new TypeError(`invalid point ‘${point}’`);
if (this.height!=0 || point.height!=0) throw new RangeError('point must be on the surface of the ellipsoid');
const p1 = this, p2 = point;
const φ1 = p1.lat.toRadians(), λ1 = p1.lon.toRadians();
const φ2 = p2.lat.toRadians(), λ2 = p2.lon.toRadians();
// allow alternative ellipsoid to be specified
const ellipsoid = this.datum ? this.datum.ellipsoid : LatLonEllipsoidal.ellipsoids.WGS84;
const { a, b, f } = ellipsoid;
const L = λ2 - λ1; // L = difference in longitude, U = reduced latitude, defined by tan U = (1-f)·tanφ.
const tanU1 = (1-f) * Math.tan(φ1), cosU1 = 1 / Math.sqrt((1 + tanU1*tanU1)), sinU1 = tanU1 * cosU1;
const tanU2 = (1-f) * Math.tan(φ2), cosU2 = 1 / Math.sqrt((1 + tanU2*tanU2)), sinU2 = tanU2 * cosU2;
const antipodal = Math.abs(L) > π/2 || Math.abs(φ2-φ1) > π/2;
let λ = L, sinλ = null, cosλ = null; // λ = difference in longitude on an auxiliary sphere
let σ = antipodal ? π : 0, sinσ = 0, cosσ = antipodal ? -1 : 1, sinSqσ = null; // σ = angular distance P₁ P₂ on the sphere
let cos2σₘ = 1; // σₘ = angular distance on the sphere from the equator to the midpoint of the line
let cosSqα = 1; // α = azimuth of the geodesic at the equator
let λʹ = null, iterations = 0;
do {
sinλ = Math.sin(λ);
cosλ = Math.cos(λ);
sinSqσ = (cosU2*sinλ)**2 + (cosU1*sinU2-sinU1*cosU2*cosλ)**2;
if (Math.abs(sinSqσ) < ε) break; // co-incident/antipodal points (falls back on λ/σ = L)
sinσ = Math.sqrt(sinSqσ);
cosσ = sinU1*sinU2 + cosU1*cosU2*cosλ;
σ = Math.atan2(sinσ, cosσ);
const sinα = cosU1 * cosU2 * sinλ / sinσ;
cosSqα = 1 - sinα*sinα;
cos2σₘ = (cosSqα != 0) ? (cosσ - 2*sinU1*sinU2/cosSqα) : 0; // on equatorial line cos²α = 0 (§6)
const C = f/16*cosSqα*(4+f*(4-3*cosSqα));
λʹ = λ;
λ = L + (1-C) * f * sinα * (σ + C*sinσ*(cos2σₘ+C*cosσ*(-1+2*cos2σₘ*cos2σₘ)));
const iterationCheck = antipodal ? Math.abs(λ)-π : Math.abs(λ);
if (iterationCheck > π) throw new EvalError('λ > π');
} while (Math.abs(λ-λʹ) > 1e-12 && ++iterations<1000);
if (iterations >= 1000) throw new EvalError('Vincenty formula failed to converge');
const uSq = cosSqα * (a*a - b*b) / (b*b);
const A = (((uSq*-175 + 320)*uSq - 768)*uSq + 4096)*uSq/16384 + 1;
const B = (((uSq*-47 + 74)*uSq - 128)*uSq + 256)*uSq/1024;
const Δσ = B*sinσ*(cos2σₘ+B/4*(cosσ*(-1+2*cos2σₘ*cos2σₘ)-B/6*cos2σₘ*(-3+4*sinσ*sinσ)*(-3+4*cos2σₘ*cos2σₘ)));
const s = b*A*(σ-Δσ); // s = length of the geodesic
// note special handling of exactly antipodal points where sin²σ = 0 (due to discontinuity
// atan2(0, 0) = 0 but atan2(ε, 0) = π/2 / 90°) - in which case bearing is always meridional,
// due north (or due south!)
// α = azimuths of the geodesic; α2 the direction P₁ P₂ produced
const α1 = Math.abs(sinSqσ) < ε ? 0 : Math.atan2(cosU2*sinλ, cosU1*sinU2-sinU1*cosU2*cosλ);
const α2 = Math.abs(sinSqσ) < ε ? π : Math.atan2(cosU1*sinλ, -sinU1*cosU2+cosU1*sinU2*cosλ);
return {
distance: s,
initialBearing: Math.abs(s) < ε ? NaN : Dms.wrap360(α1.toDegrees()),
finalBearing: Math.abs(s) < ε ? NaN : Dms.wrap360(α2.toDegrees()),
iterations: iterations,
};
}
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
export { LatLonEllipsoidal_Vincenty as default, Dms };