/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* Vector-based ellipsoidal geodetic (latitude/longitude) functions (c) Chris Veness 2015-2016 */
/* MIT Licence */
/* www.movable-type.co.uk/scripts/latlong-vectors.html */
/* www.movable-type.co.uk/scripts/geodesy/docs/module-latlon-nvector-spherical.html */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
import LatLon, { Cartesian } from './latlon-ellipsoidal.js';
import Vector3d from './vector3d.js';
/**
* Tools for working with points on (ellipsoidal models of) the earth’s surface using a vector-based
* approach using ‘n-vectors’ (rather than the more common spherical trigonometry).
*
* Based on Kenneth Gade’s ‘Non-singular Horizontal Position Representation’.
*
* Note that these formulations take x => 0°N,0°E, y => 0°N,90°E, z => 90°N; Gade uses x => 90°N,
* y => 0°N,90°E, z => 0°N,0°E.
*
* @module latlon-nvector-ellipsoidal
* @requires latlon-ellipsoidal
* @requires vector3d
* @requires dms
*/
/* LatLon_NvectorEllipsoidal - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/**
* Latitude/longitude points on an ellipsoidal model earth augmented with methods for calculating
* delta vectors between points, and converting to n-vectors.
*/
class LatLon_NvectorEllipsoidal extends LatLon {
/**
* Creates a LatLon point on an ellipsoidal model earth.
*
* @param {number} lat - Latitude in degrees.
* @param {number} lon - Longitude in degrees.
* @param {number} [height=0] - Height above ellipsoid in metres.
* @param {LatLon.datums} [datum=WGS84] - Datum this point is defined within.
*
* @example
* import LatLon from 'latlon-nvector-ellipsoidal';
* var p1 = new LatLon(52.205, 0.119);
*/
constructor(lat, lon, height, datum) {
super(lat, lon, height, datum);
}
/**
* Calculates delta from ‘this’ point to supplied point.
*
* The delta is given as a north-east-down NED vector. Note that this is a linear delta,
* unrelated to a geodesic on the ellipsoid.
*
* Points need not be defined on the same datum.
*
* @param {LatLon} point - Point delta is to be determined to.
* @returns {Ned} Delta from ‘this’ point to supplied point in local tangent plane of this point.
*
* @example
* var a = new LatLon(49.66618, 3.45063);
* var b = new LatLon(48.88667, 2.37472);
* var delta = a.deltaTo(b); // [N:-86126,E:-78900,D:1069]
* var dist = delta.length; // 116807.681 m
* var brng = delta.bearing; // 222.493°
* var elev = delta.elevation; // -0.5245°
*/
deltaTo(point) {
if (!(point instanceof LatLon)) throw new TypeError('point is not LatLon object');
// get delta in cartesian frame
var c1 = this.toCartesian();
var c2 = point.toCartesian();
var δc = c2.minus(c1);
// get local (n-vector) coordinate frame
var n1 = this.toNvector();
var a = new Vector3d(0, 0, 1); // axis vector pointing to 90°N
var d = n1.negate(); // down (pointing opposite to n-vector)
var e = a.cross(n1).unit(); // east (pointing perpendicular to the plane)
var n = e.cross(d); // north (by right hand rule)
// rotation matrix is built from n-vector coordinate frame axes (using row vectors)
var r = [
[ n.x, n.y, n.z ],
[ e.x, e.y, e.z ],
[ d.x, d.y, d.z ],
];
// apply rotation to δc to get delta in n-vector reference frame
var δn = new Cartesian_Nvector(
r[0][0]*δc.x + r[0][1]*δc.y + r[0][2]*δc.z,
r[1][0]*δc.x + r[1][1]*δc.y + r[1][2]*δc.z,
r[2][0]*δc.x + r[2][1]*δc.y + r[2][2]*δc.z
);
return new Ned(δn.x, δn.y, δn.z);
}
/**
* Calculates destination point using supplied delta from ‘this’ point.
*
* @param {Ned} delta - Delta from ‘this’ point to supplied point in local tangent plane of this point.
* @returns {LatLon} Destination point.
*
* @example
* var a = new LatLon(49.66618, 3.45063);
* var delta = Ned.fromDistanceBearingElevation(116807.681, 222.493, -0.5245); // [N:-86126,E:-78900,D:1069]
* var b = a.destinationPoint(delta); // 48.88667°N, 002.37472°E
*/
destinationPoint(delta) {
if (!(delta instanceof Ned)) throw new TypeError('delta is not Ned object');
// convert North-East-Down delta to standard x/y/z vector in coordinate frame of n-vector
var δn = new Vector3d(delta.north, delta.east, delta.down);
// get local (n-vector) coordinate frame
var n1 = this.toNvector();
var a = new Vector3d(0, 0, 1); // axis vector pointing to 90°N
var d = n1.negate(); // down (pointing opposite to n-vector)
var e = a.cross(n1).unit(); // east (pointing perpendicular to the plane)
var n = e.cross(d); // north (by right hand rule)
// rotation matrix is built from n-vector coordinate frame axes (using column vectors)
var r = [
[ n.x, e.x, d.x ],
[ n.y, e.y, d.y ],
[ n.z, e.z, d.z ],
];
// apply rotation to δn to get delta in cartesian (ECEF) coordinate reference frame
var δc = new Cartesian_Nvector(
r[0][0]*δn.x + r[0][1]*δn.y + r[0][2]*δn.z,
r[1][0]*δn.x + r[1][1]*δn.y + r[1][2]*δn.z,
r[2][0]*δn.x + r[2][1]*δn.y + r[2][2]*δn.z
);
// apply (cartesian) delta to c1 to obtain destination point as cartesian coordinate
var c1 = this.toCartesian();
var c2 = c1.plus(δc);
// return destination cartesian coordinate as latitude/longitude
return c2.toLatLon();
}
/**
* Checks if another point is equal to ‘this’ point.
*
* @private
* @param {LatLon} point - Point to be compared against this point.
* @returns {bool} True if points are identical.
*
* @example
* var p1 = new LatLon(52.205, 0.119);
* var p2 = new LatLon(52.205, 0.119);
* var equal = p1.equals(p2); // true
*/
equals(point) {
if (!(point instanceof LatLon)) throw new TypeError('point is not LatLon object');
for (var p in this) {
if (this[p] != point[p]) return false; // lat, lon, height, datum
}
return true;
}
/**
* Converts ‘this’ lat/lon point to n-vector (normal to earth's surface).
*
* @returns {Nvector} N-vector representing lat/lon point.
*
* @example
* var p = new LatLon(45, 45);
* var n = p.toNvector(); // [0.5000, 0.5000, 0.7071]
*/
toNvector() { // note: replicated in LatLonNvectorSpherical
var φ = this.lat.toRadians();
var λ = this.lon.toRadians();
// right-handed vector: x -> 0°E,0°N; y -> 90°E,0°N, z -> 90°N
var x = Math.cos(φ) * Math.cos(λ);
var y = Math.cos(φ) * Math.sin(λ);
var z = Math.sin(φ);
return new NvectorEllipsoidal(x, y, z, this.height, this.datum);
}
/**
* Converts ‘this’ point from (geodetic) latitude/longitude coordinates to (geocentric) cartesian
* (x/y/z) coordinates.
*
* @returns {Cartesian} Cartesian point equivalent to lat/lon point, with x, y, z in metres from
* earth centre.
*
* TODO this simply replicates LatLon.toCartesian() from latlon-ellipsoidal.js but returns
* Cartesian_Nvector in place of Cartesian, so that toNvector() is available; there must be a
* better way of doing this!
*/
toCartesian() {
var φ = this.lat.toRadians(), λ = this.lon.toRadians();
var h = this.height; // height above ellipsoid
var a = this.datum.ellipsoid.a, b = this.datum.ellipsoid.b;
var sinφ = Math.sin(φ), cosφ = Math.cos(φ);
var sinλ = Math.sin(λ), cosλ = Math.cos(λ);
var eSq = (a*a - b*b) / (a*a); // 1st eccentricity squared
var ν = a / Math.sqrt(1 - eSq*sinφ*sinφ);
var x = (ν+h) * cosφ * cosλ;
var y = (ν+h) * cosφ * sinλ;
var z = ((1-eSq)*ν + h) * sinφ;
var p = new Cartesian_Nvector(x, y, z);
return p;
}
}
/* NvectorEllipsoidal - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/**
* An n-vector is a position representation using a (unit) vector normal to the Earth ellipsoid.
* Unlike latitude/longitude points, n-vectors have no singularities or discontinuities.
*
* For many applications, n-vectors are more convenient to work with than other position
* representations such as latitude/longitude, earth-centred earth-fixed (ECEF) vectors, UTM
* coordinates, etc.
*
* @extends vector3d
*/
class NvectorEllipsoidal extends Vector3d { // note prototype-based class not inheritance-based class
// note commonality with latlon-nvector-spherical
/**
* Creates a 3d n-vector normal to the Earth's surface.
*
* @param {number} x - x component.
* @param {number} y - y component.
* @param {number} z - z component.
* @param {number} [h=0] - height above ellipsoid surface in metres.
*
* @example
* var p1 = new Nvector(0.5000, 0.5000, 0.7071, 1); // p1.toLatLon() => 45.0000°N, 045.0000°E
*/
constructor(x, y, z, h=0, datum=LatLon.datums.WGS84) {
var u = new Vector3d(x, y, z).unit(); // n-vectors are always normalised
super(u.x, u.y, u.z);
this.h = Number(h);
this.datum = datum;
}
/**
* Converts ‘this’ n-vector to latitude/longitude point.
*
* @returns {LatLon} Latitude/longitude point equivalent to this n-vector.
*
* @example
* var v = new Nvector(0.5000, 0.5000, 0.7071);
* var p = v.toLatLon(); // 45.0°N, 45.0°E
*/
toLatLon() { // note: replicated in NvectorSpherical
var φ = Math.atan2(this.z, Math.sqrt(this.x * this.x + this.y * this.y));
var λ = Math.atan2(this.y, this.x);
return new LatLon_NvectorEllipsoidal(φ.toDegrees(), λ.toDegrees(), this.h, this.datum);
}
/**
* Converts ‘this’ n-vector to cartesian coordinate.
*
* qv Gade 2010 ‘A Non-singular Horizontal Position Representation’ eqn 22
*
* @returns {Cartesian} Cartesian coordinate equivalent to this n-vector.
*
* @example
* var v = new Nvector(0.5000, 0.5000, 0.7071);
* var c = v.toCartesian(); // [3194434,3194434,4487327]
* var p = c.toLatLon(); // 45.0°N, 45.0°E
*/
toCartesian() {
var b = this.datum.ellipsoid.b;
var f = this.datum.ellipsoid.f;
var x = this.x;
var y = this.y;
var z = this.z;
var h = this.h;
var m = (1-f) * (1-f); // (1−f)² = b²/a²
var n = b / Math.sqrt(x*x/m + y*y/m + z*z);
var xʹ = n * x / m + x*h;
var yʹ = n * y / m + y*h;
var zʹ = n * z + z*h;
return new Cartesian_Nvector(xʹ, yʹ, zʹ);
}
/**
* Returns a string representation of ‘this’ n-vector.
*
* @param {number} [dp=3] - Number of decimal places to display.
* @returns {string} Comma-separated x, y, z, h values.
*/
toString(dp=3) {
return '[' + this.x.toFixed(dp) + ',' + this.y.toFixed(dp) + ',' + this.z.toFixed(dp) + ',' + this.h.toFixed(dp) + ']';
}
}
/* Cartesian_Nvector - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/**
* Cartesian_Nvector extends Cartesian with method to convert cartesian coordinates to n-vectors, &
* methods for vector operations.
*
* @extends Vector3d,Cartesian
*/
class Cartesian_Nvector extends mixin(Vector3d, Cartesian) {
// note mixin makes Cartesian_Nvector instanceof Vector3d but not instanceof Cartesian; this is
// required to satisfy type checks in Vector3d
/**
* Converts ‘this’ cartesian coordinate to an n-vector.
*
* qv Gade 2010 ‘A Non-singular Horizontal Position Representation’ eqn 23
*
* @param {LatLon.datums} [datum=WGS84] - Datum to use for conversion.
* @returns {Nvector} N-vector equivalent to this cartesian coordinate.
*
* @example
* TODO
*/
toNvector(datum=LatLon.datums.WGS84) {
var a = datum.ellipsoid.a;
var f = datum.ellipsoid.f;
var x = this.x;
var y = this.y;
var z = this.z;
var e2 = 2*f - f*f; // e² = 1st eccentricity squared ≡ (a²-b²)/a²
var e4 = e2*e2; // e⁴
var p = (x*x + y*y) / (a*a);
var q = z*z * (1-e2) / (a*a);
var r = (p + q - e4) / 6;
var s = (e4*p*q) / (4*r*r*r);
var t = Math.cbrt(1 + s + Math.sqrt(2*s+s*s));
var u = r * (1 + t + 1/t);
var v = Math.sqrt(u*u + e4*q);
var w = e2 * (u + v - q) / (2*v);
var k = Math.sqrt(u + v + w*w) - w;
var d = k * Math.sqrt(x*x + y*y) / (k + e2);
var tmp = 1 / Math.sqrt(d*d + z*z);
var xʹ = tmp * k/(k+e2) * x;
var yʹ = tmp * k/(k+e2) * y;
var zʹ = tmp * z;
var h = (k + e2 - 1)/k * Math.sqrt(d*d + z*z);
var n = new NvectorEllipsoidal(xʹ, yʹ, zʹ, h, datum);
return n;
}
}
/* Ned - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/**
* North east down (NED), also known as local tangent plane (LTP), is a vector in the local
* coordinate frame of a body.
*/
class Ned {
/**
* Creates North-East-Down vector.
*
* @param {number} north - North component in metres.
* @param {number} east - East component in metres.
* @param {number} down - Down component (normal to the surface of the ellipsoid) in metres.
*/
constructor(north, east, down) {
this.north = north;
this.east = east;
this.down = down;
}
/**
* Length of NED vector.
*
* @returns {number} Length of NED vector in metres.
*/
get length() {
var l = Math.sqrt(this.north * this.north + this.east * this.east + this.down * this.down);
return l
}
/**
* Bearing of NED vector.
*
* @returns {number} Bearing of NED vector in degrees from north.
*/
get bearing() {
var θ = Math.atan2(this.east, this.north);
return (θ.toDegrees()+360) % 360; // normalise to 0..360°
}
/**
* Elevation of NED vector.
*
* @returns {number} Elevation of NED vector in degrees from horizontal (ie tangent to ellipsoid surface).
*/
get elevation() {
var α = Math.asin(this.down/this.length);
return -α.toDegrees();
}
/**
* Creates North-East-Down vector from distance, bearing, elevation (in local coordinate system).
*
* @param {number} dist - Length of NED vector.
* @param {number} brng - Bearing (in degrees from north) of NED vector .
* @param {number} elev - Elevation (in degrees from local coordinate frame horizontal) of NED vector.
* @returns {Ned} North-East-Down vector equivalent to distance, bearing, elevation.
*/
static fromDistanceBearingElevation(dist, brng, elev) {
var θ = Number(brng).toRadians();
var α = Number(elev).toRadians();
dist = Number(dist);
var n = Math.cos(θ) * dist*Math.cos(α);
var e = Math.sin(θ) * dist*Math.cos(α);
var d = Math.sin(-α) * dist;
return new Ned(n, e, d);
}
/**
* Returns a string representation of ‘this’ NED vector.
*
* @param {number} [dp=0] - Number of decimal places to display.
* @returns {string} Comma-separated (labelled) n, e, d values.
*/
toString(dp=0) {
return '[N:' + this.north.toFixed(dp) + ',E:' + this.east.toFixed(dp) + ',D:' + this.down.toFixed(dp) + ']';
}
}
// stackoverflow.com/questions/30732241
function mixin(target, source) {
target = target.prototype;
source = source.prototype;
Object.getOwnPropertyNames(source).forEach(function(name) {
if (name !== 'constructor') {
Object.defineProperty(target, name, Object.getOwnPropertyDescriptor(source, name));
}
});
return target.constructor;
}
// also qv hacks.mozilla.org/2015/08/es6-in-depth-subclassing
// Extend Number object with methods to convert between degrees & radians
Number.prototype.toRadians = function() { return this * Math.PI / 180; };
Number.prototype.toDegrees = function() { return this * 180 / Math.PI; };
export { LatLon_NvectorEllipsoidal as default, NvectorEllipsoidal as Nvector, Cartesian_Nvector as Cartesian, Ned };
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */