/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* Vector-based spherical geodetic (latitude/longitude) functions (c) Chris Veness 2011-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 Vector3d from './vector3d.js';
import Dms from './dms.js';
/**
* Tools for working with points and paths on (a spherical model of) the earth’s surface using a
* vector-based approach using ‘n-vectors’ (rather than the more common spherical trigonometry;
* a vector-based approach makes many calculations much simpler, and easier to follow, compared
* with trigonometric equivalents).
*
* Note on a spherical model earth, an n-vector is equivalent to a normalised version of an (ECEF)
* cartesian coordinate.
*
* @module latlon-nvector-spherical
* @requires vector3d
* @requires dms
*/
/* LatLonNvectorSpherical - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/**
* Latitude/longitude points on an spherical model earth, and methods for calculating distances,
* bearings, destinations, etc on great circle paths.
*/
class LatLonNvectorSpherical { // note prototype-based class not inheritance-based class
/**
* Creates a LatLon point on spherical model earth.
*
* @param {number} lat - Latitude in degrees.
* @param {number} lon - Longitude in degrees.
*
* @example
* import LatLon from 'latlon-nvector-spherical';
* var p1 = new LatLon(52.205, 0.119);
*/
constructor(lat, lon) {
this.lat = Number(lat);
this.lon = Number(lon);
}
/**
* Converts ‘this’ lat/lon point to n-vector (normal to earth's surface).
*
* @private
* @returns {Nvector} Normalised n-vector representing lat/lon point.
*
* @example
* var p = new LatLon(45, 45);
* var v = p.toNvector(); // [0.5000,0.5000,0.7071]
*/
toNvector() { // note: replicated in LatLonNvectorEllipsoidal
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 NvectorSpherical(x, y, z);
}
/**
* Nvector normal to great circle obtained by heading on given bearing from ‘this’ point.
*
* Direction of vector is such that initial bearing vector b = c × p.
*
* @private
* @param {number} bearing - Compass bearing in degrees.
* @returns {Nvector} Normalised vector representing great circle.
*
* @example
* var p1 = new LatLon(53.3206, -1.7297);
* var gc = p1.greatCircle(96.0); // [-0.794,0.129,0.594]
*/
greatCircle(bearing) {
var φ = this.lat.toRadians();
var λ = this.lon.toRadians();
var θ = Number(bearing).toRadians();
var x = Math.sin(λ) * Math.cos(θ) - Math.sin(φ) * Math.cos(λ) * Math.sin(θ);
var y = -Math.cos(λ) * Math.cos(θ) - Math.sin(φ) * Math.sin(λ) * Math.sin(θ);
var z = Math.cos(φ) * Math.sin(θ);
return new NvectorSpherical(x, y, z);
}
/**
* Returns the distance from ‘this’ point to the specified point.
*
* @param {LatLon} point - Latitude/longitude of destination point.
* @param {number} [radius=6371e3] - (Mean) radius of earth (defaults to radius in metres).
* @returns {number} Distance between this point and destination point, in same units as radius.
*
* @example
* var p1 = new LatLon(52.205, 0.119), p2 = new LatLon(48.857, 2.351);
* var d = p1.distanceTo(p2); // d.toPrecision(4): 404300
*/
distanceTo(point, radius=6371e3) {
if (!(point instanceof LatLonNvectorSpherical)) throw new TypeError('point is not LatLon object');
var p1 = this.toNvector();
var p2 = point.toNvector();
var δ = p1.angleTo(p2); // = atan2(|p1×p2|, p1⋅p2)
var d = δ * Number(radius);
return d;
}
/**
* Returns the (initial) bearing from ‘this’ point to the specified point, in compass degrees.
*
* @param {LatLon} point - Latitude/longitude of destination point.
* @returns {number} Initial bearing in degrees from North (0°..360°).
*
* @example
* var p1 = new LatLon(52.205, 0.119);
* var p2 = new LatLon(48.857, 2.351);
* var b1 = p1.bearingTo(p2); // 156.2°
*/
bearingTo(point) {
if (!(point instanceof LatLonNvectorSpherical)) throw new TypeError('point is not LatLon object');
var p1 = this.toNvector();
var p2 = point.toNvector();
var northPole = new NvectorSpherical(0, 0, 1);
var c1 = p1.cross(p2); // great circle through p1 & p2
var c2 = p1.cross(northPole); // great circle through p1 & north pole
// bearing is (signed) angle between c1 & c2
var bearing = c1.angleTo(c2, p1).toDegrees();
return (bearing + 360) % 360; // normalise to 0..360
}
/**
* Returns the midpoint between ‘this’ point and specified point.
*
* @param {LatLon} point - Latitude/longitude of destination point.
* @returns {LatLon} Midpoint between this point and destination point.
*
* @example
* var p1 = new LatLon(52.205, 0.119);
* var p2 = new LatLon(48.857, 2.351);
* var pMid = p1.midpointTo(p2); // 50.5363°N, 001.2746°E
*/
midpointTo(point) {
if (!(point instanceof LatLonNvectorSpherical)) throw new TypeError('point is not LatLon object');
var p1 = this.toNvector();
var p2 = point.toNvector();
var mid = p1.plus(p2);
return new NvectorSpherical(mid.x, mid.y, mid.z).toLatLon();
}
/**
* Returns the destination point from ‘this’ point having travelled the given distance on the
* given initial bearing (bearing will normally vary before destination is reached).
*
* @param {number} distance - Distance travelled, in same units as earth radius (default: metres).
* @param {number} bearing - Initial bearing in degrees from north.
* @param {number} [radius=6371e3] - (Mean) radius of earth (defaults to radius in metres).
* @returns {LatLon} Destination point.
*
* @example
* var p1 = new LatLon(51.4778, -0.0015);
* var p2 = p1.destinationPoint(7794, 300.7); // 51.5135°N, 000.0983°W
*/
destinationPoint(distance, bearing, radius=6371e3) {
var δ = Number(distance) / Number(radius); // angular distance in radians
// get great circle obtained by starting from 'this' point on given bearing
var c = this.greatCircle(Number(bearing));
var p1 = this.toNvector();
var x = p1.times(Math.cos(δ)); // component of p2 parallel to p1
var y = c.cross(p1).times(Math.sin(δ)); // component of p2 perpendicular to p1
var p2 = x.plus(y);
return new NvectorSpherical(p2.x, p2.y, p2.z).toLatLon();
}
/**
* Returns the point of intersection of two paths each defined by point pairs or start point and bearing.
*
* @param {LatLon} path1start - Start point of first path.
* @param {LatLon|number} path1brngEnd - End point of first path or initial bearing from first start point.
* @param {LatLon} path2start - Start point of second path.
* @param {LatLon|number} path2brngEnd - End point of second path or initial bearing from second start point.
* @returns {LatLon} Destination point (null if no unique intersection defined)
*
* @example
* var p1 = new LatLon(51.8853, 0.2545), brng1 = 108.55;
* var p2 = new LatLon(49.0034, 2.5735), brng2 = 32.44;
* var pInt = LatLon.intersection(p1, brng1, p2, brng2); // 50.9076°N, 004.5086°E
*/
static intersection(path1start, path1brngEnd, path2start, path2brngEnd) {
if (!(path1start instanceof LatLonNvectorSpherical)) throw new TypeError('path1start is not LatLon object');
if (!(path2start instanceof LatLonNvectorSpherical)) throw new TypeError('path2start is not LatLon object');
if (!(path1brngEnd instanceof LatLonNvectorSpherical) && isNaN(path1brngEnd)) throw new TypeError('path1brngEnd is not LatLon object or bearing');
if (!(path2brngEnd instanceof LatLonNvectorSpherical) && isNaN(path2brngEnd)) throw new TypeError('path2brngEnd is not LatLon object or bearing');
// if c1 & c2 are great circles through start and end points (or defined by start point + bearing),
// then candidate intersections are simply c1 × c2 & c2 × c1; most of the work is deciding correct
// intersection point to select! if bearing is given, that determines which intersection, if both
// paths are defined by start/end points, take closer intersection
var p1 = path1start.toNvector();
var p2 = path2start.toNvector();
var c1, c2, path1def, path2def;
// c1 & c2 are vectors defining great circles through start & end points; p × c gives initial bearing vector
if (path1brngEnd instanceof LatLonNvectorSpherical) { // path 1 defined by endpoint
c1 = p1.cross(path1brngEnd.toNvector());
path1def = 'endpoint';
} else { // path 1 defined by initial bearing
c1 = path1start.greatCircle(Number(path1brngEnd));
path1def = 'bearing';
}
if (path2brngEnd instanceof LatLonNvectorSpherical) { // path 2 defined by endpoint
c2 = p2.cross(path2brngEnd.toNvector());
path2def = 'endpoint';
} else { // path 2 defined by initial bearing
c2 = path2start.greatCircle(Number(path2brngEnd));
path2def = 'bearing';
}
// there are two (antipodal) candidate intersection points; we have to choose which to return
var i1 = c1.cross(c2);
var i2 = c2.cross(c1);
// am I making heavy weather of this? is there a simpler way to do it?
// selection of intersection point depends on how paths are defined (bearings or endpoints)
var intersection = null, dir1 = null, dir2 = null;
switch (path1def + '+' + path2def) {
case 'bearing+bearing':
// if c×p⋅i1 is +ve, the initial bearing is towards i1, otherwise towards antipodal i2
dir1 = Math.sign(c1.cross(p1).dot(i1)); // c1×p1⋅i1 +ve means p1 bearing points to i1
dir2 = Math.sign(c2.cross(p2).dot(i1)); // c2×p2⋅i1 +ve means p2 bearing points to i1
switch (dir1 + dir2) {
case 2: // dir1, dir2 both +ve, 1 & 2 both pointing to i1
intersection = i1;
break;
case -2: // dir1, dir2 both -ve, 1 & 2 both pointing to i2
intersection = i2;
break;
case 0: // dir1, dir2 opposite; intersection is at further-away intersection point
// take opposite intersection from mid-point of p1 & p2 [is this always true?]
intersection = p1.plus(p2).dot(i1) > 0 ? i2 : i1;
break;
}
break;
case 'bearing+endpoint': // use bearing c1 × p1
dir1 = Math.sign(c1.cross(p1).dot(i1)); // c1×p1⋅i1 +ve means p1 bearing points to i1
intersection = dir1 > 0 ? i1 : i2;
break;
case 'endpoint+bearing': // use bearing c2 × p2
dir2 = Math.sign(c2.cross(p2).dot(i1)); // c2×p2⋅i1 +ve means p2 bearing points to i1
intersection = dir2 > 0 ? i1 : i2;
break;
case 'endpoint+endpoint': // select nearest intersection to mid-point of all points
var mid = p1.plus(p2).plus(path1brngEnd.toNvector()).plus(path2brngEnd.toNvector());
intersection = mid.dot(i1) > 0 ? i1 : i2;
break;
}
return new NvectorSpherical(intersection.x, intersection.y, intersection.z).toLatLon();
}
/**
* Returns (signed) distance from ‘this’ point to great circle defined by start-point and end-point/bearing.
*
* @param {LatLon} pathStart - Start point of great circle path.
* @param {LatLon|number} pathBrngEnd - End point of great circle path or initial bearing from great circle start point.
* @param {number} [radius=6371e3] - (Mean) radius of earth (defaults to radius in metres).
* @returns {number} Distance to great circle (-ve if to left, +ve if to right of path).
*
* @example
* var pCurrent = new LatLon(53.2611, -0.7972);
*
* var p1 = new LatLon(53.3206, -1.7297), brng = 96.0;
* var d = pCurrent.crossTrackDistanceTo(p1, brng); // Number(d.toPrecision(4)): -305.7
*
* var p1 = new LatLon(53.3206, -1.7297), p2 = new LatLon(53.1887, 0.1334);
* var d = pCurrent.crossTrackDistanceTo(p1, p2); // Number(d.toPrecision(4)): -307.5
*/
crossTrackDistanceTo(pathStart, pathBrngEnd, radius=6371e3) {
if (!(pathStart instanceof LatLonNvectorSpherical)) throw new TypeError('pathStart is not LatLon object');
var p = this.toNvector();
var gc;
if (pathBrngEnd instanceof LatLonNvectorSpherical) {
// great circle defined by two points (note JavaScript is not good at method overloading)
var pathEnd = pathBrngEnd;
gc = pathStart.toNvector().cross(pathEnd.toNvector());
} else {
// great circle defined by point + bearing
var bearing = Number(pathBrngEnd);
gc = pathStart.greatCircle(bearing);
}
var α = gc.angleTo(p, p.cross(gc)); // (signed) angle between point & great-circle normal vector
α = α < 0 ? -Math.PI / 2 - α : Math.PI / 2 - α; // (signed) angle between point & great-circle
var d = α * Number(radius);
return d;
}
/**
* Returns closest point on great circle segment between point1 & point2 to ‘this’ point.
*
* If this point is ‘within’ the extent of the segment, the point is on the segment between point1 &
* point2; otherwise, it is the closer of the endpoints defining the segment.
*
* @param {LatLon} point1 - Start point of great circle segment.
* @param {LatLon} point2 - End point of great circle segment.
* @returns {number} Nearest point on segment.
*
* @example
* var p1 = new LatLon(51.0, 1.0), p2 = new LatLon(51.0, 2.0);
*
* var p0 = new LatLon(51.0, 1.9);
* var p = p0.nearestPointOnSegment(p1, p2); // 51.0004°N, 001.9000°E
* var d = p.distanceTo(p); // 42.71
*
* var p0 = new LatLon(51.0, 2.1);
* var p = p0.nearestPointOnSegment(p1, p2); // 51.0000°N, 002.0000°E
*/
nearestPointOnSegment(point1, point2) {
var v0 = this.toNvector(), v1 = point1.toNvector(), v2 = point2.toNvector();
// dot product p10⋅p12 tells us if p0 is on p2 side of p1, similarly for p20⋅p21
var p10 = v0.minus(v1), p12 = v2.minus(v1);
var p20 = v0.minus(v2), p21 = v1.minus(v2);
var extent1 = p10.dot(p12);
var extent2 = p20.dot(p21);
var withinExtent = extent1>=0 && extent2>=0;
var p = null;
if (withinExtent) {
// closer to segment than to its endpoints, find closest point on segment
var c1 = v1.cross(v2); // v1×v2 = vector representing great circle through p1, p2
var c2 = v0.cross(c1); // v0×c1 = vector representing great circle through p0 normal to c1
var v = c1.cross(c2); // c1×c2 = nearest point on c1 to v0
p = v.toLatLon();
} else {
// beyond segment extent, take closer endpoint
var d1 = this.distanceTo(point1);
var d2 = this.distanceTo(point2);
p = d1<d2 ? point1 : point2;
}
return p;
}
/**
* Tests whether ‘this’ point is enclosed by the (convex) polygon defined by a set of points.
*
* @param {LatLon[]} points - Ordered array of points defining vertices of polygon.
* @returns {bool} Whether this point is enclosed by region.
* @throws {Error} If polygon is not convex.
*
* @example
* var bounds = [ new LatLon(45,1), new LatLon(45,2), new LatLon(46,2), new LatLon(46,1) ];
* var p = new LatLon(45,1, 1.1);
* var inside = p.enclosedBy(bounds); // inside: true;
*/
enclosedBy(points) {
var v = this.toNvector(); // vector to 'this' point
// if fully closed polygon, pop last point off array
if (points[0].equals(points[points.length - 1])) points.pop();
// get great-circle vector for each segment
var c = [];
for (var n = 0; n < points.length; n++) {
var p1 = points[n].toNvector();
var p2 = points[n + 1 == points.length ? 0 : n + 1].toNvector();
c[n] = p1.cross(p2); // great circle for segment n
}
// is 'this' point on same side of each segment? (left/right depending on (anti-)clockwise)
var toLeft0 = c[0].angleTo(v) <= Math.PI / 2; // 'this' point to left of first segment?
for (var n = 1; n < points.length; n++) {
var toLeftN = c[n].angleTo(v) <= Math.PI / 2; // 'this' point to left of segment n?
if (toLeft0 != toLeftN) return false;
}
// is polygon convex? (otherwise above test is not reliable)
for (var n = 0; n < points.length; n++) {
var c1 = c[n];
var c2 = c[n + 1 == points.length ? 0 : n + 1];
var α = c1.angleTo(c2, v); // angle between great-circle vectors, signed by direction of v
if (α < 0) throw new Error('Polygon is not convex');
}
return true;
}
/**
* Returns point representing geographic mean of supplied points.
*
* @param {LatLon[]} points - Array of points to be averaged.
* @returns {LatLon} Point at the geographic mean of the supplied points.
*/
static meanOf(points) {
var m = new NvectorSpherical(0, 0, 0); // null vector
// add all vectors
for (var p = 0; p < points.length; p++) {
m = m.plus(points[p].toNvector());
}
// m is now geographic mean
return new NvectorSpherical(m.x, m.y, m.z).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), p2 = new LatLon(52.205, 0.119);
* var equal = p1.equals(p2); // equal: true
*/
equals(point) {
if (!(point instanceof LatLonNvectorSpherical)) throw new TypeError('point is not LatLon object');
if (this.lat != point.lat) return false;
if (this.lon != point.lon) return false;
return true;
}
/**
* Returns a string representation of ‘this’ point, formatted as degrees, degrees+minutes, or
* degrees+minutes+seconds.
*
* @param {string} [format=dms] - Format point as 'd', 'dm', 'dms'.
* @param {number} [dp=0|2|4] - Number of decimal places to use: default 0 for dms, 2 for dm, 4 for d.
* @returns {string} Comma-separated formatted latitude/longitude.
*/
toString(format='dms', dp=undefined) {
return Dms.toLat(this.lat, format, dp) + ', ' + Dms.toLon(this.lon, format, dp);
}
}
/* NvectorSpherical - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/**
* An n-vector is a vector normal to the Earth's surface (a non-singular position representation).
*
* For many applications, n-vectors are more convenient to work with than other position
* representations such as latitude/longitude, UTM coordinates, etc.
*
* On a spherical model earth, an n-vector is equivalent to an earth-centred earth-fixed (ECEF)
* vector.
*
* @extends vector3d
*/
class NvectorSpherical extends Vector3d { // note prototype-based class not inheritance-based class
// note commonality with latlon-nvector-ellipsoidal
/**
* 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 n = new Nvector(0.5000, 0.5000, 0.7071, 1); // p.toLatLon(): 45.000°N, 45.000°E
*/
constructor(x, y, z, h=0) {
var u = new Vector3d(x, y, z).unit(); // n-vectors are always normalised
super(u.x, u.y, u.z);
this.h = Number(h);
}
/**
* Converts ‘this’ n-vector to latitude/longitude point.
*
* @returns {LatLon} Latitude/longitude point vector points to.
*
* @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 NvectorEllipsoidal
var φ = Math.atan2(this.z, Math.sqrt(this.x * this.x + this.y * this.y));
var λ = Math.atan2(this.y, this.x);
return new LatLonNvectorSpherical(φ.toDegrees(), λ.toDegrees(), this.h);
}
/**
* 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) + ']';
}
}
// 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 { LatLonNvectorSpherical as default, NvectorSpherical as Nvector };
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */