diff --git a/src/index.html b/src/index.html index 3670874..5d47e03 100644 --- a/src/index.html +++ b/src/index.html @@ -5,18 +5,7 @@ Планирование маршрутов - - - - - - - - - - - diff --git a/src/js/cldr.min.js b/src/js/cldr.min.js deleted file mode 100644 index c0ffad5..0000000 --- a/src/js/cldr.min.js +++ /dev/null @@ -1,5 +0,0 @@ -/*! - * CLDR JavaScript Library v0.5.4 2020-10-22T15:56Z MIT license © Rafael Xavier - * http://git.io/h4lmVg - */ -!function(e,t){"function"==typeof define&&define.amd?define(t):"object"==typeof module&&"object"==typeof module.exports?module.exports=t():e.Cldr=t()}(this,(function(){var e,t=Array.isArray||function(e){return"[object Array]"===Object.prototype.toString.call(e)},n=function(e,n){if(t(e)&&(e=e.join("/")),"string"!=typeof e)throw new Error('invalid path "'+e+'"');return(e=(e=e.replace(/^\//,"").replace(/^cldr\//,"")).replace(/{[a-zA-Z]+}/g,(function(e){return e=e.replace(/^{([^}]*)}$/,"$1"),n[e]}))).split("/")},r=function(e,t){var n,r;if(e.some)return e.some(t);for(n=0,r=e.length;n4, 2=>8, 3=>16) - const cardinal = cardinals[Math.round(bearing*n/360)%n * 16/n]; - - return cardinal; - } - - - /** - * Constrain degrees to range -90..+90 (for latitude); e.g. -91 => -89, 91 => 89. - * - * @private - * @param {number} degrees - * @returns degrees within range -90..+90. - */ - static wrap90(degrees) { - if (-90<=degrees && degrees<=90) return degrees; // avoid rounding due to arithmetic ops if within range - - // latitude wrapping requires a triangle wave function; a general triangle wave is - // f(x) = 4a/p ⋅ | (x-p/4)%p - p/2 | - a - // where a = amplitude, p = period, % = modulo; however, JavaScript '%' is a remainder operator - // not a modulo operator - for modulo, replace 'x%n' with '((x%n)+n)%n' - const x = degrees, a = 90, p = 360; - return 4*a/p * Math.abs((((x-p/4)%p)+p)%p - p/2) - a; - } - - /** - * Constrain degrees to range -180..+180 (for longitude); e.g. -181 => 179, 181 => -179. - * - * @private - * @param {number} degrees - * @returns degrees within range -180..+180. - */ - static wrap180(degrees) { - if (-180<=degrees && degrees<=180) return degrees; // avoid rounding due to arithmetic ops if within range - - // longitude wrapping requires a sawtooth wave function; a general sawtooth wave is - // f(x) = (2ax/p - p/2) % p - a - // where a = amplitude, p = period, % = modulo; however, JavaScript '%' is a remainder operator - // not a modulo operator - for modulo, replace 'x%n' with '((x%n)+n)%n' - const x = degrees, a = 180, p = 360; - return (((2*a*x/p - p/2)%p)+p)%p - a; - } - - /** - * Constrain degrees to range 0..360 (for bearings); e.g. -1 => 359, 361 => 1. - * - * @private - * @param {number} degrees - * @returns degrees within range 0..360. - */ - static wrap360(degrees) { - if (0<=degrees && degrees<360) return degrees; // avoid rounding due to arithmetic ops if within range - - // bearing wrapping requires a sawtooth wave function with a vertical offset equal to the - // amplitude and a corresponding phase shift; this changes the general sawtooth wave function from - // f(x) = (2ax/p - p/2) % p - a - // to - // f(x) = (2ax/p) % p - // where a = amplitude, p = period, % = modulo; however, JavaScript '%' is a remainder operator - // not a modulo operator - for modulo, replace 'x%n' with '((x%n)+n)%n' - const x = degrees, a = 180, p = 360; - return (((2*a*x/p)%p)+p)%p; - } - -} - - -// 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 default Dms; diff --git a/src/js/geo/latlon-ellipsoidal-datum.js b/src/js/geo/latlon-ellipsoidal-datum.js deleted file mode 100644 index 3cc4c13..0000000 --- a/src/js/geo/latlon-ellipsoidal-datum.js +++ /dev/null @@ -1,402 +0,0 @@ -/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ -/* Geodesy tools for conversions between (historical) datums (c) Chris Veness 2005-2022 */ -/* MIT Licence */ -/* www.movable-type.co.uk/scripts/latlong-convert-coords.html */ -/* www.movable-type.co.uk/scripts/geodesy-library.html#latlon-ellipsoidal-datum */ -/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - -import LatLonEllipsoidal, { Cartesian, Dms } from './latlon-ellipsoidal.js'; - - -/** - * Historical geodetic datums: a latitude/longitude point defines a geographic location on or - * above/below the earth’s surface, measured in degrees from the equator & the International - * Reference Meridian and metres above the ellipsoid, and based on a given datum. The datum is - * based on a reference ellipsoid and tied to geodetic survey reference points. - * - * Modern geodesy is generally based on the WGS84 datum (as used for instance by GPS systems), but - * previously various reference ellipsoids and datum references were used. - * - * This module extends the core latlon-ellipsoidal module to include ellipsoid parameters and datum - * transformation parameters, and methods for converting between different (generally historical) - * datums. - * - * It can be used for UK Ordnance Survey mapping (OS National Grid References are still based on the - * otherwise historical OSGB36 datum), as well as for historical purposes. - * - * q.v. Ordnance Survey ‘A guide to coordinate systems in Great Britain’ Section 6, - * www.ordnancesurvey.co.uk/docs/support/guide-coordinate-systems-great-britain.pdf, and also - * www.ordnancesurvey.co.uk/blog/2014/12/2. - * - * @module latlon-ellipsoidal-datum - */ - - -/* - * Ellipsoid parameters; exposed through static getter below. - */ -const ellipsoids = { - WGS84: { a: 6378137, b: 6356752.314245, f: 1/298.257223563 }, - Airy1830: { a: 6377563.396, b: 6356256.909, f: 1/299.3249646 }, - AiryModified: { a: 6377340.189, b: 6356034.448, f: 1/299.3249646 }, - Bessel1841: { a: 6377397.155, b: 6356078.962822, f: 1/299.15281285 }, - Clarke1866: { a: 6378206.4, b: 6356583.8, f: 1/294.978698214 }, - Clarke1880IGN: { a: 6378249.2, b: 6356515.0, f: 1/293.466021294 }, - GRS80: { a: 6378137, b: 6356752.314140, f: 1/298.257222101 }, - Intl1924: { a: 6378388, b: 6356911.946128, f: 1/297 }, // aka Hayford - WGS72: { a: 6378135, b: 6356750.52, f: 1/298.26 }, -}; - - -/* - * Datums; exposed through static getter below. - */ -const datums = { - // transforms: t in metres, s in ppm, r in arcseconds tx ty tz s rx ry rz - ED50: { ellipsoid: ellipsoids.Intl1924, transform: [ 89.5, 93.8, 123.1, -1.2, 0.0, 0.0, 0.156 ] }, // epsg.io/1311 - ETRS89: { ellipsoid: ellipsoids.GRS80, transform: [ 0, 0, 0, 0, 0, 0, 0 ] }, // epsg.io/1149; @ 1-metre level - Irl1975: { ellipsoid: ellipsoids.AiryModified, transform: [ -482.530, 130.596, -564.557, -8.150, 1.042, 0.214, 0.631 ] }, // epsg.io/1954 - NAD27: { ellipsoid: ellipsoids.Clarke1866, transform: [ 8, -160, -176, 0, 0, 0, 0 ] }, - NAD83: { ellipsoid: ellipsoids.GRS80, transform: [ 0.9956, -1.9103, -0.5215, -0.00062, 0.025915, 0.009426, 0.011599 ] }, - NTF: { ellipsoid: ellipsoids.Clarke1880IGN, transform: [ 168, 60, -320, 0, 0, 0, 0 ] }, - OSGB36: { ellipsoid: ellipsoids.Airy1830, transform: [ -446.448, 125.157, -542.060, 20.4894, -0.1502, -0.2470, -0.8421 ] }, // epsg.io/1314 - Potsdam: { ellipsoid: ellipsoids.Bessel1841, transform: [ -582, -105, -414, -8.3, 1.04, 0.35, -3.08 ] }, - TokyoJapan: { ellipsoid: ellipsoids.Bessel1841, transform: [ 148, -507, -685, 0, 0, 0, 0 ] }, - WGS72: { ellipsoid: ellipsoids.WGS72, transform: [ 0, 0, -4.5, -0.22, 0, 0, 0.554 ] }, - WGS84: { ellipsoid: ellipsoids.WGS84, transform: [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ] }, -}; -/* sources: - * - ED50: www.gov.uk/guidance/oil-and-gas-petroleum-operations-notices#pon-4 - * - Irl1975: www.osi.ie/wp-content/uploads/2015/05/transformations_booklet.pdf - * - NAD27: en.wikipedia.org/wiki/Helmert_transformation - * - NAD83: www.uvm.edu/giv/resources/WGS84_NAD83.pdf [strictly, WGS84(G1150) -> NAD83(CORS96) @ epoch 1997.0] - * (note NAD83(1986) ≡ WGS84(Original); confluence.qps.nl/pages/viewpage.action?pageId=29855173) - * - NTF: Nouvelle Triangulation Francaise geodesie.ign.fr/contenu/fichiers/Changement_systeme_geodesique.pdf - * - OSGB36: www.ordnancesurvey.co.uk/docs/support/guide-coordinate-systems-great-britain.pdf - * - Potsdam: kartoweb.itc.nl/geometrics/Coordinate%20transformations/coordtrans.html - * - TokyoJapan: www.geocachingtoolbox.com?page=datumEllipsoidDetails - * - WGS72: www.icao.int/safety/pbn/documentation/eurocontrol/eurocontrol wgs 84 implementation manual.pdf - * - * more transform parameters are available from earth-info.nga.mil/GandG/coordsys/datums/NATO_DT.pdf, - * www.fieldenmaps.info/cconv/web/cconv_params.js - */ -/* note: - * - ETRS89 reference frames are coincident with WGS-84 at epoch 1989.0 (ie null transform) at the one metre level. - */ - - -// freeze static properties -Object.keys(ellipsoids).forEach(e => Object.freeze(ellipsoids[e])); -Object.keys(datums).forEach(d => { Object.freeze(datums[d]); Object.freeze(datums[d].transform); }); - - -/* LatLon - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - - -/** - * Latitude/longitude points on an ellipsoidal model earth, with ellipsoid parameters and methods - * for converting between datums and to geocentric (ECEF) cartesian coordinates. - * - * @extends LatLonEllipsoidal - */ -class LatLonEllipsoidal_Datum extends LatLonEllipsoidal { - - /** - * Creates a geodetic latitude/longitude point on an ellipsoidal model earth using given datum. - * - * @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 - Datum this point is defined within. - * - * @example - * import LatLon from '/js/geodesy/latlon-ellipsoidal-datum.js'; - * const p = new LatLon(53.3444, -6.2577, 17, LatLon.datums.Irl1975); - */ - constructor(lat, lon, height=0, datum=datums.WGS84) { - if (!datum || datum.ellipsoid==undefined) throw new TypeError(`unrecognised datum ‘${datum}’`); - - super(lat, lon, height); - - this._datum = datum; - } - - - /** - * Datum this point is defined within. - */ - get datum() { - return this._datum; - } - - - /** - * Ellipsoids with their parameters; semi-major axis (a), semi-minor axis (b), and flattening (f). - * - * Flattening f = (a−b)/a; at least one of these parameters is derived from defining constants. - * - * @example - * const a = LatLon.ellipsoids.Airy1830.a; // 6377563.396 - */ - static get ellipsoids() { - return ellipsoids; - } - - - /** - * Datums; with associated ellipsoid, and Helmert transform parameters to convert from WGS-84 - * into given datum. - * - * Note that precision of various datums will vary, and WGS-84 (original) is not defined to be - * accurate to better than ±1 metre. No transformation should be assumed to be accurate to - * better than a metre, for many datums somewhat less. - * - * This is a small sample of commoner datums from a large set of historical datums. I will add - * new datums on request. - * - * @example - * const a = LatLon.datums.OSGB36.ellipsoid.a; // 6377563.396 - * const tx = LatLon.datums.OSGB36.transform; // [ tx, ty, tz, s, rx, ry, rz ] - * const availableDatums = Object.keys(LatLon.datums).join(', '); // ED50, Irl1975, NAD27, ... - */ - static get datums() { - return datums; - } - - - // note instance datum getter/setters are in LatLonEllipsoidal - - - /** - * Parses a latitude/longitude point from a variety of formats. - * - * Latitude & longitude (in degrees) can be supplied as two separate parameters, as a single - * comma-separated lat/lon string, or as a single object with { lat, lon } or GeoJSON properties. - * - * The latitude/longitude values may be numeric or strings; they may be signed decimal or - * deg-min-sec (hexagesimal) suffixed by compass direction (NSEW); a variety of separators are - * accepted. Examples -3.62, '3 37 12W', '3°37′12″W'. - * - * Thousands/decimal separators must be comma/dot; use Dms.fromLocale to convert locale-specific - * thousands/decimal separators. - * - * @param {number|string|Object} lat|latlon - Geodetic Latitude (in degrees) or comma-separated lat/lon or lat/lon object. - * @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. - * @returns {LatLon} Latitude/longitude point on ellipsoidal model earth using given datum. - * @throws {TypeError} Unrecognised datum. - * - * @example - * const p = LatLon.parse('51.47736, 0.0000', 0, LatLon.datums.OSGB36); - */ - static parse(...args) { - let datum = datums.WGS84; - - // if the last argument is a datum, use that, otherwise use default WGS-84 - if (args.length==4 || (args.length==3 && typeof args[2] == 'object')) datum = args.pop(); - - if (!datum || datum.ellipsoid==undefined) throw new TypeError(`unrecognised datum ‘${datum}’`); - - const point = super.parse(...args); - - point._datum = datum; - - return point; - } - - - /** - * Converts ‘this’ lat/lon coordinate to new coordinate system. - * - * @param {LatLon.datums} toDatum - Datum this coordinate is to be converted to. - * @returns {LatLon} This point converted to new datum. - * @throws {TypeError} Unrecognised datum. - * - * @example - * const pWGS84 = new LatLon(51.47788, -0.00147, 0, LatLon.datums.WGS84); - * const pOSGB = pWGS84.convertDatum(LatLon.datums.OSGB36); // 51.4773°N, 000.0001°E - */ - convertDatum(toDatum) { - if (!toDatum || toDatum.ellipsoid==undefined) throw new TypeError(`unrecognised datum ‘${toDatum}’`); - - const oldCartesian = this.toCartesian(); // convert geodetic to cartesian - const newCartesian = oldCartesian.convertDatum(toDatum); // convert datum - const newLatLon = newCartesian.toLatLon(); // convert cartesian back to geodetic - - return newLatLon; - } - - - /** - * Converts ‘this’ point from (geodetic) latitude/longitude coordinates to (geocentric) cartesian - * (x/y/z) coordinates, based on the same datum. - * - * Shadow of LatLonEllipsoidal.toCartesian(), returning Cartesian augmented with - * LatLonEllipsoidal_Datum methods/properties. - * - * @returns {Cartesian} Cartesian point equivalent to lat/lon point, with x, y, z in metres from - * earth centre, augmented with reference frame conversion methods and properties. - */ - toCartesian() { - const cartesian = super.toCartesian(); - const cartesianDatum = new Cartesian_Datum(cartesian.x, cartesian.y, cartesian.z, this.datum); - return cartesianDatum; - } - -} - - -/* Cartesian - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - - -/** - * Augments Cartesian with datum the cooordinate is based on, and methods to convert between datums - * (using Helmert 7-parameter transforms) and to convert cartesian to geodetic latitude/longitude - * point. - * - * @extends Cartesian - */ -class Cartesian_Datum extends Cartesian { - - /** - * Creates cartesian coordinate representing ECEF (earth-centric earth-fixed) point, on a given - * datum. The datum will identify the primary meridian (for the x-coordinate), and is also - * useful in transforming to/from geodetic (lat/lon) coordinates. - * - * @param {number} x - X coordinate in metres (=> 0°N,0°E). - * @param {number} y - Y coordinate in metres (=> 0°N,90°E). - * @param {number} z - Z coordinate in metres (=> 90°N). - * @param {LatLon.datums} [datum] - Datum this coordinate is defined within. - * @throws {TypeError} Unrecognised datum. - * - * @example - * import { Cartesian } from '/js/geodesy/latlon-ellipsoidal-datum.js'; - * const coord = new Cartesian(3980581.210, -111.159, 4966824.522); - */ - constructor(x, y, z, datum=undefined) { - if (datum && datum.ellipsoid==undefined) throw new TypeError(`unrecognised datum ‘${datum}’`); - - super(x, y, z); - - if (datum) this._datum = datum; - } - - - /** - * Datum this point is defined within. - */ - get datum() { - return this._datum; - } - set datum(datum) { - if (!datum || datum.ellipsoid==undefined) throw new TypeError(`unrecognised datum ‘${datum}’`); - this._datum = datum; - } - - - /** - * Converts ‘this’ (geocentric) cartesian (x/y/z) coordinate to (geodetic) latitude/longitude - * point (based on the same datum, or WGS84 if unset). - * - * Shadow of Cartesian.toLatLon(), returning LatLon augmented with LatLonEllipsoidal_Datum - * methods convertDatum, toCartesian, etc. - * - * @returns {LatLon} Latitude/longitude point defined by cartesian coordinates. - * @throws {TypeError} Unrecognised datum - * - * @example - * const c = new Cartesian(4027893.924, 307041.993, 4919474.294); - * const p = c.toLatLon(); // 50.7978°N, 004.3592°E - */ - toLatLon(deprecatedDatum=undefined) { - if (deprecatedDatum) { - console.info('datum parameter to Cartesian_Datum.toLatLon is deprecated: set datum before calling toLatLon()'); - this.datum = deprecatedDatum; - } - const datum = this.datum || datums.WGS84; - if (!datum || datum.ellipsoid==undefined) throw new TypeError(`unrecognised datum ‘${datum}’`); - - const latLon = super.toLatLon(datum.ellipsoid); // TODO: what if datum is not geocentric? - const point = new LatLonEllipsoidal_Datum(latLon.lat, latLon.lon, latLon.height, this.datum); - return point; - } - - - /** - * Converts ‘this’ cartesian coordinate to new datum using Helmert 7-parameter transformation. - * - * @param {LatLon.datums} toDatum - Datum this coordinate is to be converted to. - * @returns {Cartesian} This point converted to new datum. - * @throws {Error} Undefined datum. - * - * @example - * const c = new Cartesian(3980574.247, -102.127, 4966830.065, LatLon.datums.OSGB36); - * c.convertDatum(LatLon.datums.Irl1975); // [??,??,??] - */ - convertDatum(toDatum) { - // TODO: what if datum is not geocentric? - if (!toDatum || toDatum.ellipsoid == undefined) throw new TypeError(`unrecognised datum ‘${toDatum}’`); - if (!this.datum) throw new TypeError('cartesian coordinate has no datum'); - - let oldCartesian = null; - let transform = null; - - if (this.datum == undefined || this.datum == datums.WGS84) { - // converting from WGS 84 - oldCartesian = this; - transform = toDatum.transform; - } - if (toDatum == datums.WGS84) { - // converting to WGS 84; use inverse transform - oldCartesian = this; - transform = this.datum.transform.map(p => -p); - } - if (transform == null) { - // neither this.datum nor toDatum are WGS84: convert this to WGS84 first - oldCartesian = this.convertDatum(datums.WGS84); - transform = toDatum.transform; - } - - const newCartesian = oldCartesian.applyTransform(transform); - newCartesian.datum = toDatum; - - return newCartesian; - } - - - /** - * Applies Helmert 7-parameter transformation to ‘this’ coordinate using transform parameters t. - * - * This is used in converting datums (geodetic->cartesian, apply transform, cartesian->geodetic). - * - * @private - * @param {number[]} t - Transformation to apply to this coordinate. - * @returns {Cartesian} Transformed point. - */ - applyTransform(t) { - // this point - const { x: x1, y: y1, z: z1 } = this; - - // transform parameters - const tx = t[0]; // x-shift in metres - const ty = t[1]; // y-shift in metres - const tz = t[2]; // z-shift in metres - const s = t[3]/1e6 + 1; // scale: normalise parts-per-million to (s+1) - const rx = (t[4]/3600).toRadians(); // x-rotation: normalise arcseconds to radians - const ry = (t[5]/3600).toRadians(); // y-rotation: normalise arcseconds to radians - const rz = (t[6]/3600).toRadians(); // z-rotation: normalise arcseconds to radians - - // apply transform - const x2 = tx + x1*s - y1*rz + z1*ry; - const y2 = ty + x1*rz + y1*s - z1*rx; - const z2 = tz - x1*ry + y1*rx + z1*s; - - return new Cartesian_Datum(x2, y2, z2); - } -} - - -/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - -export { LatLonEllipsoidal_Datum as default, Cartesian_Datum as Cartesian, datums, Dms }; diff --git a/src/js/geo/latlon-ellipsoidal-referenceframe-txparams.js b/src/js/geo/latlon-ellipsoidal-referenceframe-txparams.js deleted file mode 100644 index 0da6656..0000000 --- a/src/js/geo/latlon-ellipsoidal-referenceframe-txparams.js +++ /dev/null @@ -1,148 +0,0 @@ -/* Helmert transform parameters tx(mm) ty(mm) tz(mm) s(ppb) rx(mas) ry(mas) rz(mas) */ -export default { - /* eslint-disable key-spacing, indent */ - 'ITRF2014→ITRF2008': { epoch: '2010.0', - params: [ 1.6, 1.9, 2.4, -0.02, 0.00, 0.00, 0.00 ], - rates: [ 0.0, 0.0, -0.1, 0.03, 0.00, 0.00, 0.00 ] }, - 'ITRF2014→ITRF2005': { epoch: '2010.0', - params: [ 2.6, 1.0, -2.3, 0.92, 0.00, 0.00, 0.00 ], - rates: [ 0.3, 0.0, -0.1, 0.03, 0.00, 0.00, 0.00 ] }, - 'ITRF2014→ITRF2000': { epoch: '2010.0', - params: [ 0.7, 1.2, -26.1, 2.12, 0.00, 0.00, 0.00 ], - rates: [ 0.1, 0.1, -1.9, 0.11, 0.00, 0.00, 0.00 ] }, - 'ITRF2014→ITRF97': { epoch: '2010.0', - params: [ 7.4, -0.5, -62.8, 3.80, 0.00, 0.00, 0.26 ], - rates: [ 0.1, -0.5, -3.3, 0.12, 0.00, 0.00, 0.02 ] }, - 'ITRF2014→ITRF96': { epoch: '2010.0', - params: [ 7.4, -0.5, -62.8, 3.80, 0.00, 0.00, 0.26 ], - rates: [ 0.1, -0.5, -3.3, 0.12, 0.00, 0.00, 0.02 ] }, - 'ITRF2014→ITRF94': { epoch: '2010.0', - params: [ 7.4, -0.5, -62.8, 3.80, 0.00, 0.00, 0.26 ], - rates: [ 0.1, -0.5, -3.3, 0.12, 0.00, 0.00, 0.02 ] }, - 'ITRF2014→ITRF93': { epoch: '2010.0', - params: [ -50.4, 3.3, -60.2, 4.29, -2.81, -3.38, 0.40 ], - rates: [ -2.8, -0.1, -2.5, 0.12, -0.11, -0.19, 0.07 ] }, - 'ITRF2014→ITRF92': { epoch: '2010.0', - params: [ 15.4, 1.5, -70.8, 3.09, 0.00, 0.00, 0.26 ], - rates: [ 0.1, -0.5, -3.3, 0.12, 0.00, 0.00, 0.02 ] }, - 'ITRF2014→ITRF91': { epoch: '2010.0', - params: [ 27.4, 15.5, -76.8, 4.49, 0.00, 0.00, 0.26 ], - rates: [ 0.1, -0.5, -3.3, 0.12, 0.00, 0.00, 0.02 ] }, - 'ITRF2014→ITRF90': { epoch: '2010.0', - params: [ 25.4, 11.5, -92.8, 4.79, 0.00, 0.00, 0.26 ], - rates: [ 0.1, -0.5, -3.3, 0.12, 0.00, 0.00, 0.02 ] }, - 'ITRF2014→ITRF89': { epoch: '2010.0', - params: [ 30.4, 35.5, -130.8, 8.19, 0.00, 0.00, 0.26 ], - rates: [ 0.1, -0.5, -3.3, 0.12, 0.00, 0.00, 0.02 ] }, - 'ITRF2014→ITRF88': { epoch: '2010.0', - params: [ 25.4, -0.5, -154.8, 11.29, 0.10, 0.00, 0.26 ], - rates: [ 0.1, -0.5, -3.3, 0.12, 0.00, 0.00, 0.02 ] }, - - 'ITRF2008→ITRF2005': { epoch: '2000.0', - params: [ -2.0, -0.9, -4.7, 0.94, 0.00, 0.00, 0.00 ], - rates: [ 0.3, 0.0, 0.0, 0.00, 0.00, 0.00, 0.00 ] }, - 'ITRF2008→ITRF2000': { epoch: '2000.0', - params: [ -1.9, -1.7, -10.5, 1.34, 0.00, 0.00, 0.00 ], - rates: [ 0.1, 0.1, -1.8, 0.08, 0.00, 0.00, 0.00 ] }, - 'ITRF2008→ITRF97': { epoch: '2000.0', - params: [ 4.8, 2.6, -33.2, 2.92, 0.00, 0.00, 0.06 ], - rates: [ 0.1, -0.5, -3.2, 0.09, 0.00, 0.00, 0.02 ] }, - 'ITRF2008→ITRF96': { epoch: '2000.0', - params: [ 4.8, 2.6, -33.2, 2.92, 0.00, 0.00, 0.06 ], - rates: [ 0.1, -0.5, -3.2, 0.09, 0.00, 0.00, 0.02 ] }, - 'ITRF2008→ITRF94': { epoch: '2000.0', - params: [ 4.8, 2.6, -33.2, 2.92, 0.00, 0.00, 0.06 ], - rates: [ 0.1, -0.5, -3.2, 0.09, 0.00, 0.00, 0.02 ] }, - 'ITRF2008→ITRF93': { epoch: '2000.0', - params: [ -24.0, 2.4, -38.6, 3.41, -1.71, -1.48, -0.30 ], - rates: [ -2.8, -0.1, -2.4, 0.09, -0.11, -0.19, 0.07 ] }, - 'ITRF2008→ITRF92': { epoch: '2000.0', - params: [ 12.8, 4.6, -41.2, 2.21, 0.00, 0.00, 0.06 ], - rates: [ 0.1, -0.5, -3.2, 0.09, 0.00, 0.00, 0.02 ] }, - 'ITRF2008→ITRF91': { epoch: '2000.0', - params: [ 24.8, 18.6, -47.2, 3.61, 0.00, 0.00, 0.06 ], - rates: [ 0.1, -0.5, -3.2, 0.09, 0.00, 0.00, 0.02 ] }, - 'ITRF2008→ITRF90': { epoch: '2000.0', - params: [ 22.8, 14.6, -63.2, 3.91, 0.00, 0.00, 0.06 ], - rates: [ 0.1, -0.5, -3.2, 0.09, 0.00, 0.00, 0.02 ] }, - 'ITRF2008→ITRF89': { epoch: '2000.0', - params: [ 27.8, 38.6, -101.2, 7.31, 0.00, 0.00, 0.06 ], - rates: [ 0.1, -0.5, -3.2, 0.09, 0.00, 0.00, 0.02 ] }, - 'ITRF2008→ITRF88': { epoch: '2000.0', - params: [ 22.8, 2.6, -125.2, 10.41, 0.10, 0.00, 0.06 ], - rates: [ 0.1, -0.5, -3.2, 0.09, 0.00, 0.00, 0.02 ] }, - - 'ITRF2005→ITRF2000': { epoch: '2000.0', - params: [ 0.1, -0.8, -5.8, 0.40, 0.000, 0.000, 0.000 ], - rates: [ -0.2, 0.1, -1.8, 0.08, 0.000, 0.000, 0.000 ] }, - - 'ITRF2000→ITRF97': { epoch: '1997.0', - params: [ 0.67, 0.61, -1.85, 1.55, 0.00, 0.00, 0.00 ], - rates: [ 0.00, -0.06, -0.14, 0.01, 0.00, 0.00, 0.02 ] }, - 'ITRF2000→ITRF96': { epoch: '1997.0', - params: [ 0.67, 0.61, -1.85, 1.55, 0.00, 0.00, 0.00 ], - rates: [ 0.00, -0.06, -0.14, 0.01, 0.00, 0.00, 0.02 ] }, - 'ITRF2000→ITRF94': { epoch: '1997.0', - params: [ 0.67, 0.61, -1.85, 1.55, 0.00, 0.00, 0.00 ], - rates: [ 0.00, -0.06, -0.14, 0.01, 0.00, 0.00, 0.02 ] }, - 'ITRF2000→ITRF93': { epoch: '1988.0', - params: [ 12.7, 6.5, -20.9, 1.95, -0.39, 0.80, -1.14 ], - rates: [ -2.9, -0.2, -0.6, 0.01, -0.11, -0.19, 0.07 ] }, - 'ITRF2000→ITRF92': { epoch: '1988.0', - params: [ 1.47, 1.35, -1.39, 0.75, 0.00, 0.00, -0.18 ], - rates: [ 0.00, -0.06, -0.14, 0.01, 0.00, 0.00, 0.02 ] }, - 'ITRF2000→ITRF91': { epoch: '1988.0', - params: [ 26.7, 27.5, -19.9, 2.15, 0.00, 0.00, -0.18 ], - rates: [ 0.0, -0.6, -1.4, 0.01, 0.00, 0.00, 0.02 ] }, - 'ITRF2000→ITRF90': { epoch: '1988.0', - params: [ 2.47, 2.35, -3.59, 2.45, 0.00, 0.00, -0.18 ], - rates: [ 0.00, -0.06, -0.14, 0.01, 0.00, 0.00, 0.02 ] }, - 'ITRF2000→ITRF89': { epoch: '1988.0', - params: [ 2.97, 4.75, -7.39, 5.85, 0.00, 0.00, -0.18 ], - rates: [ 0.00, -0.06, -0.14, 0.01, 0.00, 0.00, 0.02 ] }, - 'ITRF2000→ITRF88': { epoch: '1988.0', - params: [ 2.47, 1.15, -9.79, 8.95, 0.10, 0.00, -0.18 ], - rates: [ 0.00, -0.06, -0.14, 0.01, 0.00, 0.00, 0.02 ] }, - - 'ITRF2000→NAD83': { epoch: '1997.0', // note NAD83(CORS96) - params: [ 995.6, -1901.3, -521.5, 0.62, 25.915, 9.426, 11.599 ], - rates: [ 0.7, -0.7, 0.5, -0.18, 0.067, -0.757, -0.051 ] }, - - 'ITRF2014→ETRF2000': { epoch: '2000.0', - params: [ 53.7, 51.2, -55.1, 1.02, 0.891, 5.390, -8.712 ], - rates: [ 0.1, 0.1, -1.9, 0.11, 0.081, 0.490, -0.792 ] }, - 'ITRF2008→ETRF2000': { epoch: '2000.0', - params: [ 52.1, 49.3, -58.5, 1.34, 0.891, 5.390, -8.712 ], - rates: [ 0.1, 0.1, -1.8, 0.08, 0.081, 0.490, -0.792 ] }, - 'ITRF2005→ETRF2000': { epoch: '2000.0', - params: [ 54.1, 50.2, -53.8, 0.40, 0.891, 5.390, -8.712 ], - rates: [ -0.2, 0.1, -1.8, 0.08, 0.081, 0.490, -0.792 ] }, - 'ITRF2000→ETRF2000': { epoch: '2000.0', - params: [ 54.0, 51.0, -48.0, 0.00, 0.891, 5.390, -8.712 ], - rates: [ 0.0, 0.0, 0.0, 0.00, 0.081, 0.490, -0.792 ] }, - - 'ITRF2008→GDA94': { epoch: '1994.0', - params: [ -84.68, -19.42, 32.01, 9.710, -0.4254, 2.2578, 2.4015 ], - rates: [ 1.42, 1.34, 0.90, 0.109, 1.5461, 1.1820, 1.1551 ] }, - 'ITRF2005→GDA94': { epoch: '1994.0', - params: [ -79.73, -6.86, 38.03, 6.636, 0.0351, -2.1211, -2.1411 ], - rates: [ 2.25, -0.62, -0.56, 0.294, -1.4707, -1.1443, -1.1701 ] }, - 'ITRF2000→GDA94': { epoch: '1994.0', - params: [ -45.91, -29.85, -20.37, 7.070, -1.6705, 0.4594, 1.9356 ], - rates: [ -4.66, 3.55, 11.24, 0.249, 1.7454, 1.4868, 1.2240 ] }, -}; -/* Note WGS84(G730/G873/G1150) are coincident with ITRF at 10-centimetre level; WGS84(G1674) and - * ITRF20014 / ITRF2008 ‘are likely to agree at the centimeter level’ (QPS). - * - * sources: - * - ITRS: itrf.ensg.ign.fr/trans_para.php - * - NAD83: Transforming Positions and Velocities between the International Terrestrial Reference - * Frame of 2000 and North American Datum of 1983, Soler & Snay, 2004; - * www.ngs.noaa.gov/CORS/Articles/SolerSnayASCE.pdf - * - ETRS: etrs89.ensg.ign.fr/memo-V8.pdf / www.euref.eu/symposia/2016SanSebastian/01-02-Altamimi.pdf - * - GDA: ITRF to GDA94 coordinate transformations, Dawson & Woods, 2010 - * (note sign of rotations for GDA94 reversed from Dawson & Woods 2010 as “Australia assumes rotation - * to be of coordinate axes” rather than the more conventional “position around the coordinate axes”) - * more are available at: - * confluence.qps.nl/qinsy/files/en/29856813/45482834/2/1453459502000/ITRF_Transformation_Parameters.xlsx - */ diff --git a/src/js/geo/latlon-ellipsoidal-referenceframe.js b/src/js/geo/latlon-ellipsoidal-referenceframe.js deleted file mode 100644 index 1aa2773..0000000 --- a/src/js/geo/latlon-ellipsoidal-referenceframe.js +++ /dev/null @@ -1,533 +0,0 @@ -/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ -/* Geodesy tools for conversions between reference frames (c) Chris Veness 2016-2019 */ -/* MIT Licence */ -/* www.movable-type.co.uk/scripts/latlong-convert-coords.html */ -/* www.movable-type.co.uk/scripts/geodesy-library.html#latlon-ellipsoidal-referenceframe */ -/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - -import LatLonEllipsoidal, { Cartesian, Dms } from './latlon-ellipsoidal.js'; - - -/** - * Modern geodetic reference frames: a latitude/longitude point defines a geographic location on or - * above/below the earth’s surface, measured in degrees from the equator and the International - * Reference Meridian and metres above the ellipsoid within a given terrestrial reference frame at a - * given epoch. - * - * This module extends the core latlon-ellipsoidal module to include methods for converting between - * different reference frames. - * - * This is scratching the surface of complexities involved in high precision geodesy, but may be of - * interest and/or value to those with less demanding requirements. - * - * Note that ITRF solutions do not directly use an ellipsoid, but are specified by cartesian - * coordinates; the GRS80 ellipsoid is recommended for transformations to geographical coordinates - * (itrf.ensg.ign.fr). - * - * @module latlon-ellipsoidal-referenceframe - */ - - -/* - * Sources: - * - * - Soler & Snay, “Transforming Positions and Velocities between the International Terrestrial Refer- - * ence Frame of 2000 and North American Datum of 1983”, Journal of Surveying Engineering May 2004; - * www.ngs.noaa.gov/CORS/Articles/SolerSnayASCE.pdf. - * - * - Dawson & Woods, “ITRF to GDA94 coordinate transformations”, Journal of Applied Geodesy 4 (2010); - * www.ga.gov.au/webtemp/image_cache/GA19050.pdf. - */ - -/* eslint-disable key-spacing, indent */ - -/* - * Ellipsoid parameters; exposed through static getter below. - */ -const ellipsoids = { - WGS84: { a: 6378137, b: 6356752.314245, f: 1/298.257223563 }, - GRS80: { a: 6378137, b: 6356752.314140, f: 1/298.257222101 }, -}; - -/* - * Reference frames; exposed through static getter below. - */ -const referenceFrames = { - ITRF2014: { name: 'ITRF2014', epoch: 2010.0, ellipsoid: ellipsoids.GRS80 }, - ITRF2008: { name: 'ITRF2008', epoch: 2005.0, ellipsoid: ellipsoids.GRS80 }, - ITRF2005: { name: 'ITRF2005', epoch: 2000.0, ellipsoid: ellipsoids.GRS80 }, - ITRF2000: { name: 'ITRF2000', epoch: 1997.0, ellipsoid: ellipsoids.GRS80 }, - ITRF93: { name: 'ITRF93', epoch: 1988.0, ellipsoid: ellipsoids.GRS80 }, - ITRF91: { name: 'ITRF91', epoch: 1988.0, ellipsoid: ellipsoids.GRS80 }, - WGS84g1762: { name: 'WGS84g1762', epoch: 2005.0, ellipsoid: ellipsoids.WGS84 }, - WGS84g1674: { name: 'WGS84g1674', epoch: 2005.0, ellipsoid: ellipsoids.WGS84 }, - WGS84g1150: { name: 'WGS84g1150', epoch: 2001.0, ellipsoid: ellipsoids.WGS84 }, - ETRF2000: { name: 'ETRF2000', epoch: 2005.0, ellipsoid: ellipsoids.GRS80 }, // ETRF2000(R08) - NAD83: { name: 'NAD83', epoch: 1997.0, ellipsoid: ellipsoids.GRS80 }, // CORS96 - GDA94: { name: 'GDA94', epoch: 1994.0, ellipsoid: ellipsoids.GRS80 }, -}; - -/* - * Transform parameters; exposed through static getter below. - */ -import txParams from './latlon-ellipsoidal-referenceframe-txparams.js'; - - -// freeze static properties -Object.keys(ellipsoids).forEach(e => Object.freeze(ellipsoids[e])); -Object.keys(referenceFrames).forEach(trf => Object.freeze(referenceFrames[trf])); -Object.keys(txParams).forEach(tx => { Object.freeze(txParams[tx]); Object.freeze(txParams[tx].params); Object.freeze(txParams[tx].rates); }); - -/* eslint-enable key-spacing, indent */ - - -/* LatLon - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - - -/** - * Latitude/longitude points on an ellipsoidal model earth, with ellipsoid parameters and methods - * for converting between reference frames and to geocentric (ECEF) cartesian coordinates. - * - * @extends LatLonEllipsoidal - */ -class LatLonEllipsoidal_ReferenceFrame extends LatLonEllipsoidal { - - /** - * Creates geodetic latitude/longitude point on an ellipsoidal model earth using using a - * specified reference frame. - * - * Note that while the epoch defaults to the frame reference epoch, the accuracy of ITRF - * realisations is meaningless without knowing the observation epoch. - * - * @param {number} lat - Geodetic latitude in degrees. - * @param {number} lon - Geodetic longitude in degrees. - * @param {number} [height=0] - Height above ellipsoid in metres. - * @param {LatLon.referenceFrames} [referenceFrame=ITRF2014] - Reference frame this point is defined within. - * @param {number} [epoch=referenceFrame.epoch] - date of observation of coordinate (decimal year). - * defaults to reference epoch t₀ of reference frame. - * @throws {TypeError} Unrecognised reference frame. - * - * @example - * import LatLon from '/js/geodesy/latlon-ellipsoidal-referenceframe.js'; - * const p = new LatLon(51.47788, -0.00147, 0, LatLon.referenceFrames.ITRF2000); - */ - constructor(lat, lon, height=0, referenceFrame=referenceFrames.ITRF2014, epoch=undefined) { - if (!referenceFrame || referenceFrame.epoch==undefined) throw new TypeError('unrecognised reference frame'); - if (epoch != undefined && isNaN(Number(epoch))) throw new TypeError(`invalid epoch ’${epoch}’`); - - super(lat, lon, height); - - this._referenceFrame = referenceFrame; - if (epoch) this._epoch = Number(epoch); - } - - - /** - * Reference frame this point is defined within. - */ - get referenceFrame() { - return this._referenceFrame; - } - - - /** - * Point’s observed epoch. - */ - get epoch() { - return this._epoch || this.referenceFrame.epoch; - } - - - /** - * Ellipsoid parameters; semi-major axis (a), semi-minor axis (b), and flattening (f). - * - * The only ellipsoids used in modern geodesy are WGS-84 and GRS-80 (while based on differing - * defining parameters, the only effective difference is a 0.1mm variation in the minor axis b). - * - * @example - * const availableEllipsoids = Object.keys(LatLon.ellipsoids).join(); // WGS84,GRS80 - * const a = LatLon.ellipsoids.Airy1830.a; // 6377563.396 - */ - static get ellipsoids() { - return ellipsoids; - } - - - /** - * Reference frames, with their base ellipsoids and reference epochs. - * - * @example - * const availableReferenceFrames = Object.keys(LatLon.referenceFrames).join(); // ITRF2014,ITRF2008, ... - */ - static get referenceFrames() { - return referenceFrames; - } - - - /** - * 14-parameter Helmert transformation parameters between (dynamic) ITRS frames, and from ITRS - * frames to (static) regional TRFs NAD83, ETRF2000, and GDA94. - * - * This is a limited set of transformations; e.g. ITRF frames prior to ITRF2000 are not included. - * More transformations could be added on request. - * - * Many conversions are direct; for NAD83, successive ITRF transformations are chained back to - * ITRF2000. - */ - static get transformParameters() { - return txParams; - } - - - /** - * Parses a latitude/longitude point from a variety of formats. - * - * Latitude & longitude (in degrees) can be supplied as two separate parameters, as a single - * comma-separated lat/lon string, or as a single object with { lat, lon } or GeoJSON properties. - * - * The latitude/longitude values may be numeric or strings; they may be signed decimal or - * deg-min-sec (hexagesimal) suffixed by compass direction (NSEW); a variety of separators are - * accepted. Examples -3.62, '3 37 12W', '3°37′12″W'. - * - * Thousands/decimal separators must be comma/dot; use Dms.fromLocale to convert locale-specific - * thousands/decimal separators. - * - * @param {number|string|Object} lat|latlon - Geodetic Latitude (in degrees) or comma-separated lat/lon or lat/lon object. - * @param {number} [lon] - Longitude in degrees. - * @param {number} height - Height above ellipsoid in metres. - * @param {LatLon.referenceFrames} referenceFrame - Reference frame this point is defined within. - * @param {number} [epoch=referenceFrame.epoch] - date of observation of coordinate (decimal year). - * @returns {LatLon} Latitude/longitude point on ellipsoidal model earth using given reference frame. - * @throws {TypeError} Unrecognised reference frame. - * - * @example - * const p1 = LatLon.parse(51.47788, -0.00147, 17, LatLon.referenceFrames.ETRF2000); // numeric pair - * const p2 = LatLon.parse('51°28′40″N, 000°00′05″W', 17, LatLon.referenceFrames.ETRF2000); // dms string + height - * const p3 = LatLon.parse({ lat: 52.205, lon: 0.119 }, 17, LatLon.referenceFrames.ETRF2000); // { lat, lon } object numeric - */ - static parse(...args) { - if (args.length == 0) throw new TypeError('invalid (empty) point'); - - let referenceFrame = null, epoch = null; - - if (!isNaN(args[1]) && typeof args[2] == 'object') { // latlon, height, referenceFrame, [epoch] - [ referenceFrame ] = args.splice(2, 1); - [ epoch ] = args.splice(2, 1); - } - - if (!isNaN(args[2]) && typeof args[3] == 'object') { // lat, lon, height, referenceFrame, [epoch] - [ referenceFrame ] = args.splice(3, 1); - [ epoch ] = args.splice(3, 1); - } - - if (!referenceFrame || referenceFrame.epoch==undefined) throw new TypeError('unrecognised reference frame'); - - // args is now lat, lon, height or latlon, height as taken by LatLonEllipsoidal .parse() - - const point = super.parse(...args); // note super.parse() also invokes this.constructor() - - point._referenceFrame = referenceFrame; - if (epoch) point._epoch = Number(epoch); - - return point; - } - - - /** - * Converts ‘this’ lat/lon coordinate to new coordinate system. - * - * @param {LatLon.referenceFrames} toReferenceFrame - Reference frame this coordinate is to be converted to. - * @returns {LatLon} This point converted to new reference frame. - * @throws {Error} Undefined reference frame, Transformation not available. - * - * @example - * const pEtrf = new LatLon(51.47788000, -0.00147000, 0, LatLon.referenceFrames.ITRF2000); - * const pItrf = pEtrf.convertReferenceFrame(LatLon.referenceFrames.ETRF2000); // 51.47787826°N, 000.00147125°W - */ - convertReferenceFrame(toReferenceFrame) { - if (!toReferenceFrame || toReferenceFrame.epoch == undefined) throw new TypeError('unrecognised reference frame'); - - const oldCartesian = this.toCartesian(); // convert geodetic to cartesian - const newCartesian = oldCartesian.convertReferenceFrame(toReferenceFrame); // convert TRF - const newLatLon = newCartesian.toLatLon(); // convert cartesian back to to geodetic - - return newLatLon; - } - - - /** - * Converts ‘this’ point from (geodetic) latitude/longitude coordinates to (geocentric) cartesian - * (x/y/z) coordinates, based on same reference frame. - * - * Shadow of LatLonEllipsoidal.toCartesian(), returning Cartesian augmented with - * LatLonEllipsoidal_ReferenceFrame methods/properties. - * - * @returns {Cartesian} Cartesian point equivalent to lat/lon point, with x, y, z in metres from - * earth centre, augmented with reference frame conversion methods and properties. - */ - toCartesian() { - const cartesian = super.toCartesian(); - const cartesianReferenceFrame = new Cartesian_ReferenceFrame(cartesian.x, cartesian.y, cartesian.z, this.referenceFrame, this.epoch); - return cartesianReferenceFrame; - } - - - /** - * Returns a string representation of ‘this’ point, formatted as degrees, degrees+minutes, or - * degrees+minutes+seconds. - * - * @param {string} [format=d] - Format point as 'd', 'dm', 'dms'. - * @param {number} [dp=4|2|0] - Number of decimal places to use: default 4 for d, 2 for dm, 0 for dms. - * @param {number} [dpHeight=null] - Number of decimal places to use for height; default (null) is no height display. - * @param {boolean} [referenceFrame=false] - Whether to show reference frame point is defined on. - * @returns {string} Comma-separated formatted latitude/longitude. - * - * @example - * new LatLon(51.47788, -0.00147, 0, LatLon.referenceFrames.ITRF2014).toString(); // 51.4778°N, 000.0015°W - * new LatLon(51.47788, -0.00147, 0, LatLon.referenceFrames.ITRF2014).toString('dms'); // 51°28′40″N, 000°00′05″W - * new LatLon(51.47788, -0.00147, 42, LatLon.referenceFrames.ITRF2014).toString('dms', 0, 0); // 51°28′40″N, 000°00′05″W +42m - */ - toString(format='d', dp=undefined, dpHeight=null, referenceFrame=false) { - const ll = super.toString(format, dp, dpHeight); - - const epochFmt = { useGrouping: false, minimumFractionDigits: 1, maximumFractionDigits: 20 }; - const epoch = this.referenceFrame && this.epoch != this.referenceFrame.epoch ? this.epoch.toLocaleString('en', epochFmt) : ''; - - const trf = referenceFrame ? ` (${this.referenceFrame.name}${epoch?'@'+epoch:''})` : ''; - - return ll + trf; - } - -} - - -/* Cartesian - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - - -/** - * Augments Cartesian with reference frame and observation epoch the cooordinate is based on, and - * methods to convert between reference frames (using Helmert 14-parameter transforms) and to - * convert cartesian to geodetic latitude/longitude point. - * - * @extends Cartesian - */ -class Cartesian_ReferenceFrame extends Cartesian { - - /** - * Creates cartesian coordinate representing ECEF (earth-centric earth-fixed) point, on a given - * reference frame. The reference frame will identify the primary meridian (for the x-coordinate), - * and is also useful in transforming to/from geodetic (lat/lon) coordinates. - * - * @param {number} x - X coordinate in metres (=> 0°N,0°E). - * @param {number} y - Y coordinate in metres (=> 0°N,90°E). - * @param {number} z - Z coordinate in metres (=> 90°N). - * @param {LatLon.referenceFrames} [referenceFrame] - Reference frame this coordinate is defined within. - * @param {number} [epoch=referenceFrame.epoch] - date of observation of coordinate (decimal year). - * @throws {TypeError} Unrecognised reference frame, Invalid epoch. - * - * @example - * import { Cartesian } from '/js/geodesy/latlon-ellipsoidal-referenceframe.js'; - * const coord = new Cartesian(3980581.210, -111.159, 4966824.522); - */ - constructor(x, y, z, referenceFrame=undefined, epoch=undefined) { - if (referenceFrame!=undefined && referenceFrame.epoch==undefined) throw new TypeError('unrecognised reference frame'); - if (epoch!=undefined && isNaN(Number(epoch))) throw new TypeError(`invalid epoch ’${epoch}’`); - - super(x, y, z); - - if (referenceFrame) this._referenceFrame = referenceFrame; - if (epoch) this._epoch = epoch; - } - - - /** - * Reference frame this point is defined within. - */ - get referenceFrame() { - return this._referenceFrame; - } - set referenceFrame(referenceFrame) { - if (!referenceFrame || referenceFrame.epoch==undefined) throw new TypeError('unrecognised reference frame'); - this._referenceFrame = referenceFrame; - } - - /** - * Point’s observed epoch. - */ - get epoch() { - return this._epoch ? this._epoch : (this._referenceFrame ? this._referenceFrame.epoch : undefined); - } - set epoch(epoch) { - if (isNaN(Number(epoch))) throw new TypeError(`invalid epoch ’${epoch}’`); - if (this._epoch != this._referenceFrame.epoch) this._epoch = Number(epoch); - } - - - /** - * Converts ‘this’ (geocentric) cartesian (x/y/z) coordinate to (geodetic) latitude/longitude - * point (based on the same reference frame). - * - * Shadow of Cartesian.toLatLon(), returning LatLon augmented with LatLonEllipsoidal_ReferenceFrame - * methods convertReferenceFrame, toCartesian, etc. - * - * @returns {LatLon} Latitude/longitude point defined by cartesian coordinates, in given reference frame. - * @throws {Error} No reference frame defined. - * - * @example - * const c = new Cartesian(4027893.924, 307041.993, 4919474.294, LatLon.referenceFrames.ITRF2000); - * const p = c.toLatLon(); // 50.7978°N, 004.3592°E - */ - toLatLon() { - if (!this.referenceFrame) throw new Error('cartesian reference frame not defined'); - - const latLon = super.toLatLon(this.referenceFrame.ellipsoid); - const point = new LatLonEllipsoidal_ReferenceFrame(latLon.lat, latLon.lon, latLon.height, this.referenceFrame, this.epoch); - return point; - } - - - /** - * Converts ‘this’ cartesian coordinate to new reference frame using Helmert 14-parameter - * transformation. The observation epoch is unchanged. - * - * Note that different conversions have different tolerences; refer to the literature if - * tolerances are significant. - * - * @param {LatLon.referenceFrames} toReferenceFrame - Reference frame this coordinate is to be converted to. - * @returns {Cartesian} This point converted to new reference frame. - * @throws {Error} Undefined reference frame. - * - * @example - * const c = new Cartesian(3980574.247, -102.127, 4966830.065, LatLon.referenceFrames.ITRF2000); - * c.convertReferenceFrame(LatLon.referenceFrames.ETRF2000); // [3980574.395,-102.214,4966829.941](ETRF2000@1997.0) - */ - convertReferenceFrame(toReferenceFrame) { - if (!toReferenceFrame || toReferenceFrame.epoch == undefined) throw new TypeError('unrecognised reference frame'); - if (!this.referenceFrame) throw new TypeError('cartesian coordinate has no reference frame'); - - if (this.referenceFrame.name == toReferenceFrame.name) return this; // no-op! - - const oldTrf = this.referenceFrame; - const newTrf = toReferenceFrame; - - // WGS84(G730/G873/G1150) are coincident with ITRF at 10-centimetre level; WGS84(G1674) and - // ITRF20014 / ITRF2008 ‘are likely to agree at the centimeter level’ (QPS) - if (oldTrf.name.startsWith('ITRF') && newTrf.name.startsWith('WGS84')) return this; - if (oldTrf.name.startsWith('WGS84') && newTrf.name.startsWith('ITRF')) return this; - - const oldC = this; - let newC = null; - - // is requested transformation available in single step? - const txFwd = txParams[oldTrf.name+'→'+newTrf.name]; - const txRev = txParams[newTrf.name+'→'+oldTrf.name]; - - if (txFwd || txRev) { - // yes, single step available (either forward or reverse) - const tx = txFwd? txFwd : reverseTransform(txRev); - const t = this.epoch || this.referenceFrame.epoch; - const t0 = tx.epoch;//epoch || newTrf.epoch; - newC = oldC.applyTransform(tx.params, tx.rates, t-t0); // ...apply transform... - } else { - // find intermediate transform common to old & new to chain though; this is pretty yucky, - // but since with current transform params we can transform in no more than 2 steps, it works! - // TODO: find cleaner method! - const txAvailFromOld = Object.keys(txParams).filter(tx => tx.split('→')[0] == oldTrf.name).map(tx => tx.split('→')[1]); - const txAvailToNew = Object.keys(txParams).filter(tx => tx.split('→')[1] == newTrf.name).map(tx => tx.split('→')[0]); - const txIntermediateFwd = txAvailFromOld.filter(tx => txAvailToNew.includes(tx))[0]; - const txAvailFromNew = Object.keys(txParams).filter(tx => tx.split('→')[0] == newTrf.name).map(tx => tx.split('→')[1]); - const txAvailToOld = Object.keys(txParams).filter(tx => tx.split('→')[1] == oldTrf.name).map(tx => tx.split('→')[0]); - const txIntermediateRev = txAvailFromNew.filter(tx => txAvailToOld.includes(tx))[0]; - const txFwd1 = txParams[oldTrf.name+'→'+txIntermediateFwd]; - const txFwd2 = txParams[txIntermediateFwd+'→'+newTrf.name]; - const txRev1 = txParams[newTrf.name+'→'+txIntermediateRev]; - const txRev2 = txParams[txIntermediateRev+'→'+oldTrf.name]; - const tx1 = txIntermediateFwd ? txFwd1 : reverseTransform(txRev2); - const tx2 = txIntermediateFwd ? txFwd2 : reverseTransform(txRev1); - const t = this.epoch || this.referenceFrame.epoch; - newC = oldC.applyTransform(tx1.params, tx1.rates, t-tx1.epoch); // ...apply transform 1... - newC = newC.applyTransform(tx2.params, tx2.rates, t-tx2.epoch); // ...apply transform 2... - } - - newC.referenceFrame = toReferenceFrame; - newC.epoch = oldC.epoch; - - return newC; - - function reverseTransform(tx) { - return { epoch: tx.epoch, params: tx.params.map(p => -p), rates: tx.rates.map(r => -r) }; - } - } - - - /** - * Applies Helmert 14-parameter transformation to ‘this’ coordinate using supplied transform - * parameters and annual rates of change, with the secular variation given by the difference - * between the reference epoch t0 and the observation epoch tc. - * - * This is used in converting reference frames. - * - * See e.g. 3D Coordinate Transformations, Deakin, 1998. - * - * @private - * @param {number[]} params - Transform parameters tx, ty, tz, s, rx, ry, rz.. - * @param {number[]} rates - Rate of change of transform parameters ṫx, ṫy, ṫz, ṡ, ṙx, ṙy, ṙz. - * @param {number} δt - Period between reference and observed epochs, t − t₀. - * @returns {Cartesian} Transformed point (without reference frame). - */ - applyTransform(params, rates, δt) { - // this point - const x1 = this.x, y1 = this.y, z1 = this.z; - - // base parameters - const tx = params[0]/1000; // x-shift: normalise millimetres to metres - const ty = params[1]/1000; // y-shift: normalise millimetres to metres - const tz = params[2]/1000; // z-shift: normalise millimetres to metres - const s = params[3]/1e9; // scale: normalise parts-per-billion - const rx = (params[4]/3600/1000).toRadians(); // x-rotation: normalise milliarcseconds to radians - const ry = (params[5]/3600/1000).toRadians(); // y-rotation: normalise milliarcseconds to radians - const rz = (params[6]/3600/1000).toRadians(); // z-rotation: normalise milliarcseconds to radians - - // rate parameters - const ṫx = rates[0]/1000; // x-shift: normalise millimetres to metres - const ṫy = rates[1]/1000; // y-shift: normalise millimetres to metres - const ṫz = rates[2]/1000; // z-shift: normalise millimetres to metres - const ṡ = rates[3]/1e9; // scale: normalise parts-per-billion - const ṙx = (rates[4]/3600/1000).toRadians(); // x-rotation: normalise milliarcseconds to radians - const ṙy = (rates[5]/3600/1000).toRadians(); // y-rotation: normalise milliarcseconds to radians - const ṙz = (rates[6]/3600/1000).toRadians(); // z-rotation: normalise milliarcseconds to radians - - // combined (normalised) parameters - const T = { x: tx + ṫx*δt, y: ty + ṫy*δt, z: tz + ṫz*δt }; - const R = { x: rx + ṙx*δt, y: ry + ṙy*δt, z: rz + ṙz*δt }; - const S = 1 + s + ṡ*δt; - - // apply transform (shift, scale, rotate) - const x2 = T.x + x1*S - y1*R.z + z1*R.y; - const y2 = T.y + x1*R.z + y1*S - z1*R.x; - const z2 = T.z - x1*R.y + y1*R.x + z1*S; - - return new Cartesian_ReferenceFrame(x2, y2, z2); - } - - - /** - * Returns a string representation of ‘this’ cartesian point. TRF is shown if set, and - * observation epoch if different from reference epoch. - * - * @param {number} [dp=0] - Number of decimal places to use. - * @returns {string} Comma-separated latitude/longitude. - */ - toString(dp=0) { - const { x, y, z } = this; - const epochFmt = { useGrouping: false, minimumFractionDigits: 1, maximumFractionDigits: 20 }; - const epoch = this.referenceFrame && this.epoch != this.referenceFrame.epoch ? this.epoch.toLocaleString('en', epochFmt) : ''; - const trf = this.referenceFrame ? `(${this.referenceFrame.name}${epoch?'@'+epoch:''})` : ''; - return `[${x.toFixed(dp)},${y.toFixed(dp)},${z.toFixed(dp)}]${trf}`; - } -} - - -/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - -export { LatLonEllipsoidal_ReferenceFrame as default, Cartesian_ReferenceFrame as Cartesian, Dms }; diff --git a/src/js/geo/latlon-ellipsoidal-vincenty.js b/src/js/geo/latlon-ellipsoidal-vincenty.js deleted file mode 100644 index cec22d6..0000000 --- a/src/js/geo/latlon-ellipsoidal-vincenty.js +++ /dev/null @@ -1,331 +0,0 @@ -/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ -/* Vincenty Direct and Inverse Solution of Geodesics on the Ellipsoid (c) Chris Veness 2002-2022 */ -/* MIT Licence */ -/* www.ngs.noaa.gov/PUBS_LIB/inverse.pdf */ -/* 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; - if (fraction == 1) return point; - - 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 = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq))); - const B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq))); - - 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); // TV: 'iterate until negligible change in λ' (≈0.006mm) - 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σ) < 1e-24) break; // co-incident/antipodal points (σ < ≈0.006mm) - 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); // TV: 'iterate until negligible change in λ' (≈0.006mm) - if (iterations >= 1000) throw new EvalError('Vincenty formula failed to converge'); - - const uSq = cosSqα * (a*a - b*b) / (b*b); - const A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq))); - const B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq))); - 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 }; diff --git a/src/js/geo/latlon-ellipsoidal.js b/src/js/geo/latlon-ellipsoidal.js deleted file mode 100644 index 46873f6..0000000 --- a/src/js/geo/latlon-ellipsoidal.js +++ /dev/null @@ -1,429 +0,0 @@ -/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ -/* Geodesy tools for an ellipsoidal earth model (c) Chris Veness 2005-2022 */ -/* MIT Licence */ -/* Core class for latlon-ellipsoidal-datum & latlon-ellipsoidal-referenceframe. */ -/* */ -/* www.movable-type.co.uk/scripts/latlong-convert-coords.html */ -/* www.movable-type.co.uk/scripts/geodesy-library.html#latlon-ellipsoidal */ -/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - -import Dms from './dms.js'; -import Vector3d from './vector3d.js'; - - -/** - * A latitude/longitude point defines a geographic location on or above/below the earth’s surface, - * measured in degrees from the equator & the International Reference Meridian and in metres above - * the ellipsoid, and based on a given datum. - * - * As so much modern geodesy is based on WGS-84 (as used by GPS), this module includes WGS-84 - * ellipsoid parameters, and it has methods for converting geodetic (latitude/longitude) points to/from - * geocentric cartesian points; the latlon-ellipsoidal-datum and latlon-ellipsoidal-referenceframe - * modules provide transformation parameters for converting between historical datums and between - * modern reference frames. - * - * This module is used for both trigonometric geodesy (eg latlon-ellipsoidal-vincenty) and n-vector - * geodesy (eg latlon-nvector-ellipsoidal), and also for UTM/MGRS mapping. - * - * @module latlon-ellipsoidal - */ - - -/* - * Ellipsoid parameters; exposed through static getter below. - * - * The only ellipsoid defined is WGS84, for use in utm/mgrs, vincenty, nvector. - */ -const ellipsoids = { - WGS84: { a: 6378137, b: 6356752.314245, f: 1/298.257223563 }, -}; - - -/* - * Datums; exposed through static getter below. - * - * The only datum defined is WGS84, for use in utm/mgrs, vincenty, nvector. - */ -const datums = { - WGS84: { ellipsoid: ellipsoids.WGS84 }, -}; - - -// freeze static properties -Object.freeze(ellipsoids.WGS84); -Object.freeze(datums.WGS84); - - -/* LatLonEllipsoidal - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - - -/** - * Latitude/longitude points on an ellipsoidal model earth, with ellipsoid parameters and methods - * for converting points to/from cartesian (ECEF) coordinates. - * - * This is the core class, which will usually be used via LatLonEllipsoidal_Datum or - * LatLonEllipsoidal_ReferenceFrame. - */ -class LatLonEllipsoidal { - - /** - * Creates a geodetic latitude/longitude point on a (WGS84) ellipsoidal model earth. - * - * @param {number} lat - Latitude (in degrees). - * @param {number} lon - Longitude (in degrees). - * @param {number} [height=0] - Height above ellipsoid in metres. - * @throws {TypeError} Invalid lat/lon/height. - * - * @example - * import LatLon from '/js/geodesy/latlon-ellipsoidal.js'; - * const p = new LatLon(51.47788, -0.00147, 17); - */ - constructor(lat, lon, height=0) { - if (isNaN(lat) || lat == null) throw new TypeError(`invalid lat ‘${lat}’`); - if (isNaN(lon) || lon == null) throw new TypeError(`invalid lon ‘${lon}’`); - if (isNaN(height) || height == null) throw new TypeError(`invalid height ‘${height}’`); - - this._lat = Dms.wrap90(Number(lat)); - this._lon = Dms.wrap180(Number(lon)); - this._height = Number(height); - } - - - /** - * Latitude in degrees north from equator (including aliases lat, latitude): can be set as - * numeric or hexagesimal (deg-min-sec); returned as numeric. - */ - get lat() { return this._lat; } - get latitude() { return this._lat; } - set lat(lat) { - this._lat = isNaN(lat) ? Dms.wrap90(Dms.parse(lat)) : Dms.wrap90(Number(lat)); - if (isNaN(this._lat)) throw new TypeError(`invalid lat ‘${lat}’`); - } - set latitude(lat) { - this._lat = isNaN(lat) ? Dms.wrap90(Dms.parse(lat)) : Dms.wrap90(Number(lat)); - if (isNaN(this._lat)) throw new TypeError(`invalid latitude ‘${lat}’`); - } - - /** - * Longitude in degrees east from international reference meridian (including aliases lon, lng, - * longitude): can be set as numeric or hexagesimal (deg-min-sec); returned as numeric. - */ - get lon() { return this._lon; } - get lng() { return this._lon; } - get longitude() { return this._lon; } - set lon(lon) { - this._lon = isNaN(lon) ? Dms.wrap180(Dms.parse(lon)) : Dms.wrap180(Number(lon)); - if (isNaN(this._lon)) throw new TypeError(`invalid lon ‘${lon}’`); - } - set lng(lon) { - this._lon = isNaN(lon) ? Dms.wrap180(Dms.parse(lon)) : Dms.wrap180(Number(lon)); - if (isNaN(this._lon)) throw new TypeError(`invalid lng ‘${lon}’`); - } - set longitude(lon) { - this._lon = isNaN(lon) ? Dms.wrap180(Dms.parse(lon)) : Dms.wrap180(Number(lon)); - if (isNaN(this._lon)) throw new TypeError(`invalid longitude ‘${lon}’`); - } - - /** - * Height in metres above ellipsoid. - */ - get height() { return this._height; } - set height(height) { this._height = Number(height); if (isNaN(this._height)) throw new TypeError(`invalid height ‘${height}’`); } - - - /** - * Datum. - * - * Note this is replicated within LatLonEllipsoidal in order that a LatLonEllipsoidal object can - * be monkey-patched to look like a LatLonEllipsoidal_Datum, for Vincenty calculations on - * different ellipsoids. - * - * @private - */ - get datum() { return this._datum; } - set datum(datum) { this._datum = datum; } - - - /** - * Ellipsoids with their parameters; this module only defines WGS84 parameters a = 6378137, b = - * 6356752.314245, f = 1/298.257223563. - * - * @example - * const a = LatLon.ellipsoids.WGS84.a; // 6378137 - */ - static get ellipsoids() { - return ellipsoids; - } - - /** - * Datums; this module only defines WGS84 datum, hence no datum transformations. - * - * @example - * const a = LatLon.datums.WGS84.ellipsoid.a; // 6377563.396 - */ - static get datums() { - return datums; - } - - - /** - * Parses a latitude/longitude point from a variety of formats. - * - * Latitude & longitude (in degrees) can be supplied as two separate parameters, as a single - * comma-separated lat/lon string, or as a single object with { lat, lon } or GeoJSON properties. - * - * The latitude/longitude values may be numeric or strings; they may be signed decimal or - * deg-min-sec (hexagesimal) suffixed by compass direction (NSEW); a variety of separators are - * accepted. Examples -3.62, '3 37 12W', '3°37′12″W'. - * - * Thousands/decimal separators must be comma/dot; use Dms.fromLocale to convert locale-specific - * thousands/decimal separators. - * - * @param {number|string|Object} lat|latlon - Latitude (in degrees), or comma-separated lat/lon, or lat/lon object. - * @param {number} [lon] - Longitude (in degrees). - * @param {number} [height=0] - Height above ellipsoid in metres. - * @returns {LatLon} Latitude/longitude point on WGS84 ellipsoidal model earth. - * @throws {TypeError} Invalid coordinate. - * - * @example - * const p1 = LatLon.parse(51.47788, -0.00147); // numeric pair - * const p2 = LatLon.parse('51°28′40″N, 000°00′05″W', 17); // dms string + height - * const p3 = LatLon.parse({ lat: 52.205, lon: 0.119 }, 17); // { lat, lon } object numeric + height - */ - static parse(...args) { - if (args.length == 0) throw new TypeError('invalid (empty) point'); - - let lat=undefined, lon=undefined, height=undefined; - - // single { lat, lon } object - if (typeof args[0]=='object' && (args.length==1 || !isNaN(parseFloat(args[1])))) { - const ll = args[0]; - if (ll.type == 'Point' && Array.isArray(ll.coordinates)) { // GeoJSON - [ lon, lat, height ] = ll.coordinates; - height = height || 0; - } else { // regular { lat, lon } object - if (ll.latitude != undefined) lat = ll.latitude; - if (ll.lat != undefined) lat = ll.lat; - if (ll.longitude != undefined) lon = ll.longitude; - if (ll.lng != undefined) lon = ll.lng; - if (ll.lon != undefined) lon = ll.lon; - if (ll.height != undefined) height = ll.height; - lat = Dms.wrap90(Dms.parse(lat)); - lon = Dms.wrap180(Dms.parse(lon)); - } - if (args[1] != undefined) height = args[1]; - if (isNaN(lat) || isNaN(lon)) throw new TypeError(`invalid point ‘${JSON.stringify(args[0])}’`); - } - - // single comma-separated lat/lon - if (typeof args[0] == 'string' && args[0].split(',').length == 2) { - [ lat, lon ] = args[0].split(','); - lat = Dms.wrap90(Dms.parse(lat)); - lon = Dms.wrap180(Dms.parse(lon)); - height = args[1] || 0; - if (isNaN(lat) || isNaN(lon)) throw new TypeError(`invalid point ‘${args[0]}’`); - } - - // regular (lat, lon) arguments - if (lat==undefined && lon==undefined) { - [ lat, lon ] = args; - lat = Dms.wrap90(Dms.parse(lat)); - lon = Dms.wrap180(Dms.parse(lon)); - height = args[2] || 0; - if (isNaN(lat) || isNaN(lon)) throw new TypeError(`invalid point ‘${args.toString()}’`); - } - - return new this(lat, lon, height); // 'new this' as may return subclassed types - } - - - /** - * 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. - */ - toCartesian() { - // x = (ν+h)⋅cosφ⋅cosλ, y = (ν+h)⋅cosφ⋅sinλ, z = (ν⋅(1-e²)+h)⋅sinφ - // where ν = a/√(1−e²⋅sinφ⋅sinφ), e² = (a²-b²)/a² or (better conditioned) 2⋅f-f² - const ellipsoid = this.datum - ? this.datum.ellipsoid - : this.referenceFrame ? this.referenceFrame.ellipsoid : ellipsoids.WGS84; - - const φ = this.lat.toRadians(); - const λ = this.lon.toRadians(); - const h = this.height; - const { a, f } = ellipsoid; - - const sinφ = Math.sin(φ), cosφ = Math.cos(φ); - const sinλ = Math.sin(λ), cosλ = Math.cos(λ); - - const eSq = 2*f - f*f; // 1st eccentricity squared ≡ (a²-b²)/a² - const ν = a / Math.sqrt(1 - eSq*sinφ*sinφ); // radius of curvature in prime vertical - - const x = (ν+h) * cosφ * cosλ; - const y = (ν+h) * cosφ * sinλ; - const z = (ν*(1-eSq)+h) * sinφ; - - return new Cartesian(x, y, z); - } - - - /** - * Checks if another point is equal to ‘this’ point. - * - * @param {LatLon} point - Point to be compared against this point. - * @returns {bool} True if points have identical latitude, longitude, height, and datum/referenceFrame. - * @throws {TypeError} Invalid point. - * - * @example - * const p1 = new LatLon(52.205, 0.119); - * const p2 = new LatLon(52.205, 0.119); - * const equal = p1.equals(p2); // true - */ - equals(point) { - if (!(point instanceof LatLonEllipsoidal)) throw new TypeError(`invalid point ‘${point}’`); - - if (Math.abs(this.lat - point.lat) > Number.EPSILON) return false; - if (Math.abs(this.lon - point.lon) > Number.EPSILON) return false; - if (Math.abs(this.height - point.height) > Number.EPSILON) return false; - if (this.datum != point.datum) return false; - if (this.referenceFrame != point.referenceFrame) return false; - if (this.epoch != point.epoch) return false; - - return true; - } - - - /** - * Returns a string representation of ‘this’ point, formatted as degrees, degrees+minutes, or - * degrees+minutes+seconds. - * - * @param {string} [format=d] - Format point as 'd', 'dm', 'dms', or 'n' for signed numeric. - * @param {number} [dp=4|2|0] - Number of decimal places to use: default 4 for d, 2 for dm, 0 for dms. - * @param {number} [dpHeight=null] - Number of decimal places to use for height; default is no height display. - * @returns {string} Comma-separated formatted latitude/longitude. - * @throws {RangeError} Invalid format. - * - * @example - * const greenwich = new LatLon(51.47788, -0.00147, 46); - * const d = greenwich.toString(); // 51.4779°N, 000.0015°W - * const dms = greenwich.toString('dms', 2); // 51°28′40″N, 000°00′05″W - * const [lat, lon] = greenwich.toString('n').split(','); // 51.4779, -0.0015 - * const dmsh = greenwich.toString('dms', 0, 0); // 51°28′40″N, 000°00′06″W +46m - */ - toString(format='d', dp=undefined, dpHeight=null) { - // note: explicitly set dp to undefined for passing through to toLat/toLon - if (![ 'd', 'dm', 'dms', 'n' ].includes(format)) throw new RangeError(`invalid format ‘${format}’`); - - const height = (this.height>=0 ? ' +' : ' ') + this.height.toFixed(dpHeight) + 'm'; - if (format == 'n') { // signed numeric degrees - if (dp == undefined) dp = 4; - const lat = this.lat.toFixed(dp); - const lon = this.lon.toFixed(dp); - return `${lat}, ${lon}${dpHeight==null ? '' : height}`; - } - - const lat = Dms.toLat(this.lat, format, dp); - const lon = Dms.toLon(this.lon, format, dp); - - return `${lat}, ${lon}${dpHeight==null ? '' : height}`; - } - -} - - -/* Cartesian - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - - -/** - * ECEF (earth-centered earth-fixed) geocentric cartesian coordinates. - * - * @extends Vector3d - */ -class Cartesian extends Vector3d { - - /** - * Creates cartesian coordinate representing ECEF (earth-centric earth-fixed) point. - * - * @param {number} x - X coordinate in metres (=> 0°N,0°E). - * @param {number} y - Y coordinate in metres (=> 0°N,90°E). - * @param {number} z - Z coordinate in metres (=> 90°N). - * - * @example - * import { Cartesian } from '/js/geodesy/latlon-ellipsoidal.js'; - * const coord = new Cartesian(3980581.210, -111.159, 4966824.522); - */ - constructor(x, y, z) { - super(x, y, z); // arguably redundant constructor, but specifies units & axes - } - - - /** - * Converts ‘this’ (geocentric) cartesian (x/y/z) coordinate to (geodetic) latitude/longitude - * point on specified ellipsoid. - * - * Uses Bowring’s (1985) formulation for μm precision in concise form; ‘The accuracy of geodetic - * latitude and height equations’, B R Bowring, Survey Review vol 28, 218, Oct 1985. - * - * @param {LatLon.ellipsoids} [ellipsoid=WGS84] - Ellipsoid to use when converting point. - * @returns {LatLon} Latitude/longitude point defined by cartesian coordinates, on given ellipsoid. - * @throws {TypeError} Invalid ellipsoid. - * - * @example - * const c = new Cartesian(4027893.924, 307041.993, 4919474.294); - * const p = c.toLatLon(); // 50.7978°N, 004.3592°E - */ - toLatLon(ellipsoid=ellipsoids.WGS84) { - // note ellipsoid is available as a parameter for when toLatLon gets subclassed to - // Ellipsoidal_Datum / Ellipsoidal_Referenceframe. - if (!ellipsoid || !ellipsoid.a) throw new TypeError(`invalid ellipsoid ‘${ellipsoid}’`); - - const { x, y, z } = this; - const { a, b, f } = ellipsoid; - - const e2 = 2*f - f*f; // 1st eccentricity squared ≡ (a²−b²)/a² - const ε2 = e2 / (1-e2); // 2nd eccentricity squared ≡ (a²−b²)/b² - const p = Math.sqrt(x*x + y*y); // distance from minor axis - const R = Math.sqrt(p*p + z*z); // polar radius - - // parametric latitude (Bowring eqn.17, replacing tanβ = z·a / p·b) - const tanβ = (b*z)/(a*p) * (1+ε2*b/R); - const sinβ = tanβ / Math.sqrt(1+tanβ*tanβ); - const cosβ = sinβ / tanβ; - - // geodetic latitude (Bowring eqn.18: tanφ = z+ε²⋅b⋅sin³β / p−e²⋅cos³β) - const φ = isNaN(cosβ) ? 0 : Math.atan2(z + ε2*b*sinβ*sinβ*sinβ, p - e2*a*cosβ*cosβ*cosβ); - - // longitude - const λ = Math.atan2(y, x); - - // height above ellipsoid (Bowring eqn.7) - const sinφ = Math.sin(φ), cosφ = Math.cos(φ); - const ν = a / Math.sqrt(1-e2*sinφ*sinφ); // length of the normal terminated by the minor axis - const h = p*cosφ + z*sinφ - (a*a/ν); - - const point = new LatLonEllipsoidal(φ.toDegrees(), λ.toDegrees(), h); - - return point; - } - - - /** - * Returns a string representation of ‘this’ cartesian point. - * - * @param {number} [dp=0] - Number of decimal places to use. - * @returns {string} Comma-separated latitude/longitude. - */ - toString(dp=0) { - const x = this.x.toFixed(dp), y = this.y.toFixed(dp), z = this.z.toFixed(dp); - return `[${x},${y},${z}]`; - } - -} - - -/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - -export { LatLonEllipsoidal as default, Cartesian, Vector3d, Dms }; diff --git a/src/js/geo/latlon-nvector-ellipsoidal.js b/src/js/geo/latlon-nvector-ellipsoidal.js deleted file mode 100644 index c31378a..0000000 --- a/src/js/geo/latlon-nvector-ellipsoidal.js +++ /dev/null @@ -1,445 +0,0 @@ -/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ -/* Vector-based ellipsoidal geodetic (latitude/longitude) functions (c) Chris Veness 2015-2021 */ -/* MIT Licence */ -/* www.movable-type.co.uk/scripts/latlong-vectors.html */ -/* www.movable-type.co.uk/scripts/geodesy-library.html#latlon-nvector-ellipsoidal */ -/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - -import LatLonEllipsoidal, { Cartesian, Vector3d, Dms } from './latlon-ellipsoidal.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 (in order that n-vector - * = cartesian vector at 0°N,0°E); Gade uses x => 90°N, y => 0°N,90°E, z => 0°N,0°E. - * - * @module latlon-nvector-ellipsoidal - */ - - -/* LatLon_NvectorEllipsoidal - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - - -/** - * Latitude/longitude points on an ellipsoidal model earth augmented with methods for calculating - * delta vectors between points, and converting to n-vectors. - * - * @extends LatLonEllipsoidal - */ -class LatLon_NvectorEllipsoidal extends LatLonEllipsoidal { - - /** - * 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. - * @throws {TypeError} Invalid point. - * - * @example - * const a = new LatLon(49.66618, 3.45063, 99); - * const b = new LatLon(48.88667, 2.37472, 64); - * const delta = a.deltaTo(b); // [N:-86127,E:-78901,D:1104] - * const dist = delta.length; // 116809.178 m - * const brng = delta.bearing; // 222.493° - * const elev = delta.elevation; // -0.5416° - */ - deltaTo(point) { - if (!(point instanceof LatLonEllipsoidal)) throw new TypeError(`invalid point ‘${point}’`); - - // get delta in cartesian frame - const c1 = this.toCartesian(); - const c2 = point.toCartesian(); - const δc = c2.minus(c1); - - // get local (n-vector) coordinate frame - const n1 = this.toNvector(); - const a = new Vector3d(0, 0, 1); // axis vector pointing to 90°N - const d = n1.negate(); // down (pointing opposite to n-vector) - const e = a.cross(n1).unit(); // east (pointing perpendicular to the plane) - const n = e.cross(d); // north (by right hand rule) - - // rotation matrix is built from n-vector coordinate frame axes (using row vectors) - const 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 - const δn = new Cartesian( - 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. - * - * 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. - * - * @param {Ned} delta - Delta from ‘this’ point to supplied point in local tangent plane of this point. - * @returns {LatLon} Destination point. - * - * @example - * const a = new LatLon(49.66618, 3.45063, 99); - * const delta = Ned.fromDistanceBearingElevation(116809.178, 222.493, -0.5416); // [N:-86127,E:-78901,D:1104] - * const b = a.destinationPoint(delta); // 48.8867°N, 002.3747°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 - const δn = new Vector3d(delta.north, delta.east, delta.down); - - // get local (n-vector) coordinate frame - const n1 = this.toNvector(); - const a = new Vector3d(0, 0, 1); // axis vector pointing to 90°N - const d = n1.negate(); // down (pointing opposite to n-vector) - const e = a.cross(n1).unit(); // east (pointing perpendicular to the plane) - const n = e.cross(d); // north (by right hand rule) - - // rotation matrix is built from n-vector coordinate frame axes (using column vectors) - const 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 - const δc = new Cartesian( - 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 - const c1 = this.toCartesian(); // convert this LatLon to Cartesian - const v2 = c1.plus(δc); // the plus() gives us a plain vector,.. - const c2 = new Cartesian(v2.x, v2.y, v2.z); // ... need to convert it to Cartesian to get LatLon - - // return destination cartesian coordinate as latitude/longitude - return c2.toLatLon(); - } - - - /** - * Converts ‘this’ lat/lon point to n-vector (normal to the earth's surface). - * - * @returns {Nvector} N-vector representing lat/lon point. - * - * @example - * const p = new LatLon(45, 45); - * const n = p.toNvector(); // [0.5000,0.5000,0.7071] - */ - toNvector() { // note: replicated in LatLonNvectorSpherical - const φ = this.lat.toRadians(); - const λ = this.lon.toRadians(); - - const sinφ = Math.sin(φ), cosφ = Math.cos(φ); - const sinλ = Math.sin(λ), cosλ = Math.cos(λ); - - // right-handed vector: x -> 0°E,0°N; y -> 90°E,0°N, z -> 90°N - const x = cosφ * cosλ; - const y = cosφ * sinλ; - const z = sinφ; - - return new NvectorEllipsoidal(x, y, z, this.h, 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. - */ - toCartesian() { - const c = super.toCartesian(); // c is 'Cartesian' - - // return Cartesian_Nvector to have toNvector() available as method of exported LatLon - return new Cartesian_Nvector(c.x, c.y, c.z); - } - -} - - -/* Nvector - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - - -/** - * 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 commonality with latlon-nvector-spherical - - /** - * Creates a 3d n-vector normal to the Earth's surface. - * - * @param {number} x - X component of n-vector (towards 0°N, 0°E). - * @param {number} y - Y component of n-vector (towards 0°N, 90°E). - * @param {number} z - Z component of n-vector (towards 90°N). - * @param {number} [h=0] - Height above ellipsoid surface in metres. - * @param {LatLon.datums} [datum=WGS84] - Datum this n-vector is defined within. - */ - constructor(x, y, z, h=0, datum=LatLonEllipsoidal.datums.WGS84) { - const 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 - * const p = new Nvector(0.500000, 0.500000, 0.707107).toLatLon(); // 45.0000°N, 045.0000°E - */ - toLatLon() { - // tanφ = z / √(x²+y²), tanλ = y / x (same as spherical calculation) - - const { x, y, z } = this; - - const φ = Math.atan2(z, Math.sqrt(x*x + y*y)); - const λ = Math.atan2(y, 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 - * const c = new Nvector(0.500000, 0.500000, 0.707107).toCartesian(); // [3194419,3194419,4487349] - * const p = c.toLatLon(); // 45.0000°N, 045.0000°E - */ - toCartesian() { - const { b, f } = this.datum.ellipsoid; - const { x, y, z, h } = this; - - const m = (1-f) * (1-f); // (1−f)² = b²/a² - const n = b / Math.sqrt(x*x/m + y*y/m + z*z); - - const xʹ = n * x / m + x*h; - const yʹ = n * y / m + y*h; - const zʹ = n * z + z*h; - - return new Cartesian_Nvector(xʹ, yʹ, zʹ); - } - - - /** - * Returns a string representation of ‘this’ (unit) n-vector. Height component is only shown if - * dpHeight is specified. - * - * @param {number} [dp=3] - Number of decimal places to display. - * @param {number} [dpHeight=null] - Number of decimal places to use for height; default is no height display. - * @returns {string} Comma-separated x, y, z, h values. - * - * @example - * new Nvector(0.5000, 0.5000, 0.7071).toString(); // [0.500,0.500,0.707] - * new Nvector(0.5000, 0.5000, 0.7071, 1).toString(6, 0); // [0.500002,0.500002,0.707103+1m] - */ - toString(dp=3, dpHeight=null) { - const { x, y, z } = this; - const h = `${this.h>=0 ? '+' : ''}${this.h.toFixed(dpHeight)}m`; - - return `[${x.toFixed(dp)},${y.toFixed(dp)},${z.toFixed(dp)}${dpHeight==null ? '' : h}]`; - } - -} - - -/* Cartesian - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - - -/** - * Cartesian_Nvector extends Cartesian with method to convert cartesian coordinates to n-vectors. - * - * @extends Cartesian - */ -class Cartesian_Nvector extends Cartesian { - - - /** - * 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 - * const c = new Cartesian(3980581, 97, 4966825); - * const n = c.toNvector(); // { x: 0.6228, y: 0.0000, z: 0.7824, h: 0.0000 } - */ - toNvector(datum=LatLonEllipsoidal.datums.WGS84) { - const { a, f } = datum.ellipsoid; - const { x, y, z } = this; - - const e2 = 2*f - f*f; // e² = 1st eccentricity squared ≡ (a²-b²)/a² - const e4 = e2*e2; // e⁴ - - const p = (x*x + y*y) / (a*a); - const q = z*z * (1-e2) / (a*a); - const r = (p + q - e4) / 6; - const s = (e4*p*q) / (4*r*r*r); - const t = Math.cbrt(1 + s + Math.sqrt(2*s+s*s)); - const u = r * (1 + t + 1/t); - const v = Math.sqrt(u*u + e4*q); - const w = e2 * (u + v - q) / (2*v); - const k = Math.sqrt(u + v + w*w) - w; - const d = k * Math.sqrt(x*x + y*y) / (k + e2); - - const tmp = 1 / Math.sqrt(d*d + z*z); - const xʹ = tmp * k/(k+e2) * x; - const yʹ = tmp * k/(k+e2) * y; - const zʹ = tmp * z; - const h = (k + e2 - 1)/k * Math.sqrt(d*d + z*z); - - const 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. - * - * @example - * import { Ned } from '/js/geodesy/latlon-nvector-ellipsoidal.js'; - * const delta = new Ned(110569, 111297, 1936); // [N:110569,E:111297,D:1936] - */ - 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() { - const { north, east, down } = this; - - return Math.sqrt(north*north + east*east + down*down); - } - - - /** - * Bearing of NED vector. - * - * @returns {number} Bearing of NED vector in degrees from north. - */ - get bearing() { - const θ = Math.atan2(this.east, this.north); - - return Dms.wrap360(θ.toDegrees()); // normalise to range 0..360° - } - - - /** - * Elevation of NED vector. - * - * @returns {number} Elevation of NED vector in degrees from horizontal (ie tangent to ellipsoid surface). - */ - get elevation() { - const α = 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 in metres. - * @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. - * - * @example - * const delta = Ned.fromDistanceBearingElevation(116809.178, 222.493, -0.5416); // [N:-86127,E:-78901,D:1104] - */ - static fromDistanceBearingElevation(dist, brng, elev) { - const θ = Number(brng).toRadians(); - const α = Number(elev).toRadians(); - dist = Number(dist); - - const sinθ = Math.sin(θ), cosθ = Math.cos(θ); - const sinα = Math.sin(α), cosα = Math.cos(α); - - const n = cosθ * dist*cosα; - const e = sinθ * dist*cosα; - const d = -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)}]`; - } - -} - - -/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - -export { LatLon_NvectorEllipsoidal as default, NvectorEllipsoidal as Nvector, Cartesian_Nvector as Cartesian, Ned, Dms }; diff --git a/src/js/geo/latlon-nvector-spherical.js b/src/js/geo/latlon-nvector-spherical.js deleted file mode 100644 index 4d94c43..0000000 --- a/src/js/geo/latlon-nvector-spherical.js +++ /dev/null @@ -1,1021 +0,0 @@ -/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ -/* Vector-based spherical geodetic (latitude/longitude) functions (c) Chris Veness 2011-2022 */ -/* MIT Licence */ -/* www.movable-type.co.uk/scripts/latlong-vectors.html */ -/* www.movable-type.co.uk/scripts/geodesy-library.html#latlon-nvector-spherical */ -/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - -import Vector3d from './vector3d.js'; -import Dms from './dms.js'; - -const π = Math.PI; - - -/** - * Tools for working with points and paths on (a spherical model of) the earth’s surface using a - * vector-based approach using ‘n-vectors’. In contrast to the more common spherical trigonometry, - * a vector-based approach makes many calculations much simpler, and easier to follow. - * - * 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. - * - * Note also that on a spherical model earth, an n-vector is equivalent to a normalised version of - * an (ECEF) cartesian coordinate. - * - * @module latlon-nvector-spherical - */ - - -/* LatLonNvectorSpherical - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - - -/** - * Latitude/longitude points on an spherical model earth, and methods for calculating distances, - * bearings, destinations, etc on great circle paths. - */ -class LatLonNvectorSpherical { - - /** - * Creates a latitude/longitude point on the earth’s surface, using a spherical model earth. - * - * @param {number} lat - Latitude (in degrees). - * @param {number} lon - Longitude (in degrees). - * @throws {TypeError} Invalid lat/lon. - * - * @example - * import LatLon from '/js/geodesy/latlon-nvector-spherical.js'; - * const p = new LatLon(52.205, 0.119); - */ - constructor(lat, lon) { - if (isNaN(lat)) throw new TypeError(`invalid lat ‘${lat}’`); - if (isNaN(lon)) throw new TypeError(`invalid lon ‘${lon}’`); - - this._lat = Dms.wrap90(Number(lat)); - this._lon = Dms.wrap180(Number(lon)); - } - - - /** - * Latitude in degrees north from equator (including aliases lat, latitude): can be set as - * numeric or hexagesimal (deg-min-sec); returned as numeric. - */ - get lat() { return this._lat; } - get latitude() { return this._lat; } - set lat(lat) { - this._lat = isNaN(lat) ? Dms.wrap90(Dms.parse(lat)) : Dms.wrap90(Number(lat)); - if (isNaN(this._lat)) throw new TypeError(`invalid lat ‘${lat}’`); - } - set latitude(lat) { - this._lat = isNaN(lat) ? Dms.wrap90(Dms.parse(lat)) : Dms.wrap90(Number(lat)); - if (isNaN(this._lat)) throw new TypeError(`invalid latitude ‘${lat}’`); - } - - /** - * Longitude in degrees east from international reference meridian (including aliases lon, lng, - * longitude): can be set as numeric or hexagesimal (deg-min-sec); returned as numeric. - */ - get lon() { return this._lon; } - get lng() { return this._lon; } - get longitude() { return this._lon; } - set lon(lon) { - this._lon = isNaN(lon) ? Dms.wrap180(Dms.parse(lon)) : Dms.wrap180(Number(lon)); - if (isNaN(this._lon)) throw new TypeError(`invalid lon ‘${lon}’`); - } - set lng(lon) { - this._lon = isNaN(lon) ? Dms.wrap180(Dms.parse(lon)) : Dms.wrap180(Number(lon)); - if (isNaN(this._lon)) throw new TypeError(`invalid lng ‘${lon}’`); - } - set longitude(lon) { - this._lon = isNaN(lon) ? Dms.wrap180(Dms.parse(lon)) : Dms.wrap180(Number(lon)); - if (isNaN(this._lon)) throw new TypeError(`invalid longitude ‘${lon}’`); - } - - - /** Conversion factors; 1000 * LatLon.metresToKm gives 1. */ - static get metresToKm() { return 1/1000; } - /** Conversion factors; 1000 * LatLon.metresToMiles gives 0.621371192237334. */ - static get metresToMiles() { return 1/1609.344; } - /** Conversion factors; 1000 * LatLon.metresToMiles gives 0.5399568034557236. */ - static get metresToNauticalMiles() { return 1/1852; } - - - // TODO: is it worth LatLon.parse() for the n-vector version? - - - /** - * Converts ‘this’ latitude/longitude point to an n-vector (normal to earth's surface). - * - * @returns {Nvector} Normalised n-vector representing lat/lon point. - * - * @example - * const p = new LatLon(45, 45); - * const v = p.toNvector(); // [0.5000,0.5000,0.7071] - */ - toNvector() { // note: replicated in LatLon_NvectorEllipsoidal - const φ = this.lat.toRadians(); - const λ = this.lon.toRadians(); - - const sinφ = Math.sin(φ), cosφ = Math.cos(φ); - const sinλ = Math.sin(λ), cosλ = Math.cos(λ); - - // right-handed vector: x -> 0°E,0°N; y -> 90°E,0°N, z -> 90°N - const x = cosφ * cosλ; - const y = cosφ * sinλ; - const z = sinφ; - - return new NvectorSpherical(x, y, z); - } - - - /** - * Vector normal to great circle obtained by heading on given bearing from ‘this’ point. - * - * Direction of vector is such that initial bearing vector b = c × n, where n is an n-vector - * representing ‘this’ (start) point. - * - * @private - * @param {number} bearing - Compass bearing in degrees. - * @returns {Vector3d} Normalised vector representing great circle. - * - * @example - * const p1 = new LatLon(53.3206, -1.7297); - * const gc = p1.greatCircle(96.0); // [-0.794,0.129,0.594] - */ - greatCircle(bearing) { - const φ = this.lat.toRadians(); - const λ = this.lon.toRadians(); - const θ = Number(bearing).toRadians(); - - const x = Math.sin(λ) * Math.cos(θ) - Math.sin(φ) * Math.cos(λ) * Math.sin(θ); - const y = -Math.cos(λ) * Math.cos(θ) - Math.sin(φ) * Math.sin(λ) * Math.sin(θ); - const z = Math.cos(φ) * Math.sin(θ); - - return new Vector3d(x, y, z); - } - - - /** - * Returns the distance on the surface of the sphere from ‘this’ point to destination point. - * - * @param {LatLon} point - Latitude/longitude of destination point. - * @param {number} [radius=6371e3] - Radius of earth (defaults to mean radius in metres). - * @returns {number} Distance between this point and destination point, in same units as radius. - * @throws {TypeError} Invalid point/radius. - * - * @example - * const p1 = new LatLon(52.205, 0.119); - * const p2 = new LatLon(48.857, 2.351); - * const d = p1.distanceTo(p2); // 404.3 km - */ - distanceTo(point, radius=6371e3) { - if (!(point instanceof LatLonNvectorSpherical)) throw new TypeError(`invalid point ‘${point}’`); - if (isNaN(radius)) throw new TypeError(`invalid radius ‘${radius}’`); - - const R = Number(radius); - - const n1 = this.toNvector(); - const n2 = point.toNvector(); - - const sinθ = n1.cross(n2).length; - const cosθ = n1.dot(n2); - const δ = Math.atan2(sinθ, cosθ); // tanδ = |n₁×n₂| / n₁⋅n₂ - - return δ * R; - } - - - /** - * Returns the initial bearing from ‘this’ point to destination point. - * - * @param {LatLon} point - Latitude/longitude of destination point. - * @returns {number} Initial bearing in degrees from north (0°..360°). - * @throws {TypeError} Invalid point. - * - * @example - * const p1 = new LatLon(52.205, 0.119); - * const p2 = new LatLon(48.857, 2.351); - * const b1 = p1.initialBearingTo(p2); // 156.2° - */ - initialBearingTo(point) { - if (!(point instanceof LatLonNvectorSpherical)) throw new TypeError(`invalid point ‘${point}’`); - if (this.equals(point)) return NaN; // coincident points - - const p1 = this.toNvector(); - const p2 = point.toNvector(); - - const N = new NvectorSpherical(0, 0, 1); // n-vector representing north pole - - const c1 = p1.cross(p2); // great circle through p1 & p2 - const c2 = p1.cross(N); // great circle through p1 & north pole - - const θ = c1.angleTo(c2, p1); // bearing is (signed) angle between c1 & c2 - - return Dms.wrap360(θ.toDegrees()); // normalise to range 0..360° - } - - - /** - * Returns final bearing arriving at destination point from ‘this’ point; the final bearing will - * differ from the initial bearing by varying degrees according to distance and latitude. - * - * @param {LatLon} point - Latitude/longitude of destination point. - * @returns {number} Final bearing in degrees from north (0°..360°). - * @throws {TypeError} Invalid point. - * - * @example - * const p1 = new LatLon(52.205, 0.119); - * const p2 = new LatLon(48.857, 2.351); - * const b2 = p1.finalBearingTo(p2); // 157.9° - */ - finalBearingTo(point) { - if (!(point instanceof LatLonNvectorSpherical)) throw new TypeError(`invalid point ‘${point}’`); - - // get initial bearing from destination point to this point & reverse it by adding 180° - return Dms.wrap360(point.initialBearingTo(this) + 180); - } - - - /** - * Returns the midpoint between ‘this’ point and destination point. - * - * @param {LatLon} point - Latitude/longitude of destination point. - * @returns {LatLon} Midpoint between this point and destination point. - * @throws {TypeError} Invalid point. - * - * @example - * const p1 = new LatLon(52.205, 0.119); - * const p2 = new LatLon(48.857, 2.351); - * const pMid = p1.midpointTo(p2); // 50.5363°N, 001.2746°E - */ - midpointTo(point) { - if (!(point instanceof LatLonNvectorSpherical)) throw new TypeError(`invalid point ‘${point}’`); - - const n1 = this.toNvector(); - const n2 = point.toNvector(); - - const mid = n1.plus(n2); - - return new NvectorSpherical(mid.x, mid.y, mid.z).toLatLon(); - } - - - /** - * 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. - * @throws {TypeError} Invalid point/fraction. - * - * @example - * const p1 = new LatLon(52.205, 0.119); - * const p2 = new LatLon(48.857, 2.351); - * const pInt = p1.intermediatePointTo(p2, 0.25); // 51.3721°N, 000.7072°E - */ - intermediatePointTo(point, fraction) { - if (!(point instanceof LatLonNvectorSpherical)) throw new TypeError(`invalid point ‘${point}’`); - if (isNaN(fraction)) throw new TypeError(`invalid fraction ‘${fraction}’`); - - // angular distance between points; tanδ = |n₁×n₂| / n₁⋅n₂ - const n1 = this.toNvector(); - const n2 = point.toNvector(); - const sinθ = n1.cross(n2).length; - const cosθ = n1.dot(n2); - const δ = Math.atan2(sinθ, cosθ); - - // interpolated angular distance on straight line between points - const δi = δ * Number(fraction); - const sinδi = Math.sin(δi); - const cosδi = Math.cos(δi); - - // direction vector (perpendicular to n1 in plane of n2) - const d = n1.cross(n2).unit().cross(n1); // unit(n₁×n₂) × n₁ - - // interpolated position - const int = n1.times(cosδi).plus(d.times(sinδi)); // n₁⋅cosδᵢ + d⋅sinδᵢ - - return new NvectorSpherical(int.x, int.y, int.z).toLatLon(); - } - - - /** - * Returns the latitude/longitude point projected from the point at given fraction on a straight - * line between 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. - * @throws {TypeError} Invalid point. - * - * @example - * const p1 = new LatLon(52.205, 0.119); - * const p2 = new LatLon(48.857, 2.351); - * const pInt = p1.intermediatePointTo(p2, 0.25); // 51.3723°N, 000.7072°E - */ - intermediatePointOnChordTo(point, fraction) { - if (!(point instanceof LatLonNvectorSpherical)) throw new TypeError(`invalid point ‘${point}’`); - - const n1 = this.toNvector(); - const n2 = point.toNvector(); - - const int = n1.plus(n2.minus(n1).times(Number(fraction))); // n₁ + (n₂−n₁)·f ≡ n₁·(1-f) + n₂·f - - const n = new NvectorSpherical(int.x, int.y, int.z); - - return n.toLatLon(); - } - - - /** - * Returns the destination point from ‘this’ point having travelled the given distance on the - * given initial bearing (bearing normally varies around path followed). - * - * @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 - * const p1 = new LatLon(51.47788, -0.00147); - * const p2 = p1.destinationPoint(7794, 300.7); // 51.5136°N, 000.0983°W - */ - destinationPoint(distance, bearing, radius=6371e3) { - const n1 = this.toNvector(); // Gade's n_EA_E - const δ = distance / radius; // angular distance in radians - const θ = Number(bearing).toRadians(); // initial bearing in radians - - const N = new NvectorSpherical(0, 0, 1); // north pole - - const de = N.cross(n1).unit(); // east direction vector @ n1 (Gade's k_e_E) - const dn = n1.cross(de); // north direction vector @ n1 (Gade's (k_n_E) - - const deSinθ = de.times(Math.sin(θ)); - const dnCosθ = dn.times(Math.cos(θ)); - - const d = dnCosθ.plus(deSinθ); // direction vector @ n1 (≡ C×n1; C = great circle) - - const x = n1.times(Math.cos(δ)); // component of n2 parallel to n1 - const y = d.times(Math.sin(δ)); // component of n2 perpendicular to n1 - - const n2 = x.plus(y); // Gade's n_EB_E - - return new NvectorSpherical(n2.x, n2.y, n2.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) - * @throws {TypeError} Invalid parameter. - * - * @example - * const p1 = new LatLon(51.8853, 0.2545), brng1 = 108.55; - * const p2 = new LatLon(49.0034, 2.5735), brng2 = 32.44; - * const 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(`invalid path1start ‘${path1start}’`); - if (!(path2start instanceof LatLonNvectorSpherical)) throw new TypeError(`invalid path2start ‘${path2start}’`); - if (!(path1brngEnd instanceof LatLonNvectorSpherical) && isNaN(path1brngEnd)) throw new TypeError(`invalid path1brngEnd ‘${path1brngEnd}’`); - if (!(path2brngEnd instanceof LatLonNvectorSpherical) && isNaN(path2brngEnd)) throw new TypeError(`invalid path2brngEnd ‘${path2brngEnd}’`); - - if (path1start.equals(path2start)) return new LatLonNvectorSpherical(path1start.lat, path2start.lon); // coincident points - - // 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 - - const p1 = path1start.toNvector(); - const p2 = path2start.toNvector(); - - let c1 = null, c2 = null, path1def = null, path2def = null; - // 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(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(path2brngEnd); - path2def = 'bearing'; - } - - // there are two (antipodal) candidate intersection points; we have to choose which to return - const i1 = c1.cross(c2); - const i2 = c2.cross(c1); - - // TODO 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) - let 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 - const mid = p1.plus(p2).plus(path1brngEnd.toNvector()).plus(path2brngEnd.toNvector()); // eslint-disable-line no-case-declarations - 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). - * @throws {TypeError} Invalid parameter. - * - * @example - * const pCurrent = new LatLon(53.2611, -0.7972); - * - * const p1 = new LatLon(53.3206, -1.7297), brng = 96.0; - * const d = pCurrent.crossTrackDistanceTo(p1, brng); // Number(d.toPrecision(4)): -305.7 - * - * const p1 = new LatLon(53.3206, -1.7297), p2 = new LatLon(53.1887, 0.1334); - * const d = pCurrent.crossTrackDistanceTo(p1, p2); // Number(d.toPrecision(4)): -307.5 - */ - crossTrackDistanceTo(pathStart, pathBrngEnd, radius=6371e3) { - if (!(pathStart instanceof LatLonNvectorSpherical)) throw new TypeError(`invalid pathStart ‘${pathStart}’`); - if (!(pathBrngEnd instanceof LatLonNvectorSpherical || !isNaN(pathBrngEnd))) throw new TypeError(`invalid pathBrngEnd ‘${pathBrngEnd}’`); - - if (this.equals(pathStart)) return 0; - - const p = this.toNvector(); - const R = Number(radius); - - const gc = pathBrngEnd instanceof LatLonNvectorSpherical // (note JavaScript is not good at method overloading) - ? pathStart.toNvector().cross(pathBrngEnd.toNvector()) // great circle defined by two points - : pathStart.greatCircle(pathBrngEnd); // great circle defined by point + bearing - - const α = gc.angleTo(p) - π/2; // angle between point & great-circle - - return α * R; - } - - - /** - * Returns how far ‘this’ point is along a path from from start-point, heading on bearing or towards - * end-point. That is, if a perpendicular is drawn from ‘this’ point to the (great circle) path, the - * along-track distance is the distance from the start point to where the perpendicular crosses the - * path. - * - * @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 along great circle to point nearest ‘this’ point. - * - * @example - * const pCurrent = new LatLon(53.2611, -0.7972); - * const p1 = new LatLon(53.3206, -1.7297); - * const p2 = new LatLon(53.1887, 0.1334); - * const d = pCurrent.alongTrackDistanceTo(p1, p2); // 62.331 km - */ - alongTrackDistanceTo(pathStart, pathBrngEnd, radius=6371e3) { - if (!(pathStart instanceof LatLonNvectorSpherical)) throw new TypeError(`invalid pathStart ‘${pathStart}’`); - if (!(pathBrngEnd instanceof LatLonNvectorSpherical || !isNaN(pathBrngEnd))) throw new TypeError(`invalid pathBrngEnd ‘${pathBrngEnd}’`); - - const p = this.toNvector(); - const R = Number(radius); - - const gc = pathBrngEnd instanceof LatLonNvectorSpherical // (note JavaScript is not good at method overloading) - ? pathStart.toNvector().cross(pathBrngEnd.toNvector()) // great circle defined by two points - : pathStart.greatCircle(pathBrngEnd); // great circle defined by point + bearing - - const pat = gc.cross(p).cross(gc); // along-track point c × p × c - - const α = pathStart.toNvector().angleTo(pat, gc); // angle between start point and along-track point - - return α * R; - } - - - /** - * 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 {LatLon} Closest point on segment. - * - * @example - * const p1 = new LatLon(51.0, 1.0); - * const p2 = new LatLon(51.0, 2.0); - * - * const p0 = new LatLon(51.0, 1.9); - * const p = p0.nearestPointOnSegment(p1, p2); // 51.0004°N, 001.9000°E - * const d = p.distanceTo(p); // 42.71 m - * - * const p0 = new LatLon(51.0, 2.1); - * const p = p0.nearestPointOnSegment(p1, p2); // 51.0000°N, 002.0000°E - */ - nearestPointOnSegment(point1, point2) { - let p = null; - - if (this.isWithinExtent(point1, point2) && !point1.equals(point2)) { - // closer to segment than to its endpoints, find closest point on segment - const n0 = this.toNvector(), n1 = point1.toNvector(), n2 = point2.toNvector(); - const c1 = n1.cross(n2); // n1×n2 = vector representing great circle through p1, p2 - const c2 = n0.cross(c1); // n0×c1 = vector representing great circle through p0 normal to c1 - const n = c1.cross(c2); // c2×c1 = nearest point on c1 to n0 - p = new NvectorSpherical(n.x, n.y, n.z).toLatLon(); - } else { - // beyond segment extent, take closer endpoint - const d1 = this.distanceTo(point1); - const d2 = this.distanceTo(point2); - const pCloser = d1p1, p0->p2, p1->p2, p2->p1 - const δ10 = n0.minus(n1), δ12 = n2.minus(n1); - const δ20 = n0.minus(n2), δ21 = n1.minus(n2); - - // dot product δ10⋅δ12 tells us if p0 is on p2 side of p1, similarly for δ20⋅δ21 - const extent1 = δ10.dot(δ12); - const extent2 = δ20.dot(δ21); - - const isSameHemisphere = n0.dot(n1)>=0 && n0.dot(n2)>=0; - - return extent1>=0 && extent2>=0 && isSameHemisphere; - } - - - /** - * Locates a point given two known locations and bearings from those locations. - * - * @param {LatLon} point1 - First reference point. - * @param {number} bearing1 - Bearing (in degrees from north) from first reference point. - * @param {LatLon} point2 - Second reference point. - * @param {number} bearing2 - Bearing (in degrees from north) from second reference point. - * @returns {LatLon} Triangulated point. - * - * @example - * const p1 = new LatLon(50.7175,1.65139), p2 = new LatLon(50.9250,1.7094); - * const p = LatLon.triangulate(p1, 333.3508, p2, 310.1414); // 51.1297°N, 001.3214°E - */ - static triangulate(point1, bearing1, point2, bearing2) { - const n1 = point1.toNvector(), θ1 = Number(bearing1).toRadians(); - const n2 = point2.toNvector(), θ2 = Number(bearing2).toRadians(); - - const N = new NvectorSpherical(0, 0, 1); // north pole - - const de1 = N.cross(n1).unit(); // east vector @ n1 - const dn1 = n1.cross(de1); // north vector @ n1 - const de1Sinθ = de1.times(Math.sin(θ1)); - const dn1Cosθ = dn1.times(Math.cos(θ1)); - const d1 = dn1Cosθ.plus(de1Sinθ); // direction vector @ n1 - - const c1 = n1.cross(d1); // great circle p1 + bearing1 - - const de2 = N.cross(n2).unit(); // east vector @ n2 - const dn2 = n2.cross(de2); // north vector @ n2 - const de2Sinθ = de2.times(Math.sin(θ2)); - const dn2Cosθ = dn2.times(Math.cos(θ2)); - const d2 = dn2Cosθ.plus(de2Sinθ); // direction vector @ n2 - - const c2 = n2.cross(d2); // great circle p2 + bearing2 - - const ni = c1.cross(c2); // n-vector of intersection point - - return new NvectorSpherical(ni.x, ni.y, ni.z).toLatLon(); - } - - - /** - * Locates a latitude/longitude point at given distances from three other points. - * - * @param {LatLon} point1 - First reference point. - * @param {number} distance1 - Distance to first reference point (same units as radius). - * @param {LatLon} point2 - Second reference point. - * @param {number} distance2 - Distance to second reference point (same units as radius). - * @param {LatLon} point3 - Third reference point. - * @param {number} distance3 - Distance to third reference point (same units as radius). - * @param {number} [radius=6371e3] - (Mean) radius of earth (defaults to radius in metres). - * @returns {LatLon} Trilaterated point. - * - * @example - * LatLon.trilaterate(new LatLon(0, 0), 157e3, new LatLon(0, 1), 111e3, new LatLon(1, 0), 111e3); // 00.9985°N, 000.9986°E - */ - static trilaterate(point1, distance1, point2, distance2, point3, distance3, radius=6371e3) { - // from en.wikipedia.org/wiki/Trilateration - - const n1 = point1.toNvector(), δ1 = Number(distance1)/Number(radius); - const n2 = point2.toNvector(), δ2 = Number(distance2)/Number(radius); - const n3 = point3.toNvector(), δ3 = Number(distance3)/Number(radius); - - // the following uses x,y coordinate system with origin at n1, x axis n1->n2 - const eX = n2.minus(n1).unit(); // unit vector in x direction n1->n2 - const i = eX.dot(n3.minus(n1)); // signed magnitude of x component of n1->n3 - const eY = n3.minus(n1).minus(eX.times(i)).unit(); // unit vector in y direction - const d = n2.minus(n1).length; // distance n1->n2 - const j = eY.dot(n3.minus(n1)); // signed magnitude of y component of n1->n3 - const x = (δ1*δ1 - δ2*δ2 + d*d) / (2*d); // x component of n1 -> intersection - const y = (δ1*δ1 - δ3*δ3 + i*i + j*j) / (2*j) - x*i/j; // y component of n1 -> intersection - // const eZ = eX.cross(eY); // unit vector perpendicular to plane - // const z = Math.sqrt(δ1*δ1 - x*x - y*y); // z will be NaN for no intersections - - if (!isFinite(x) || !isFinite(y)) return null; // coincident points? - - const n = n1.plus(eX.times(x)).plus(eY.times(y)); // note don't use z component; assume points at same height - - return new NvectorSpherical(n.x, n.y, n.z).toLatLon(); - } - - - - /** - * Tests whether ‘this’ point is enclosed by the polygon defined by a set of points. - * - * @param {LatLon[]} polygon - Ordered array of points defining vertices of polygon. - * @returns {bool} Whether this point is enclosed by polygon. - * - * @example - * const bounds = [ new LatLon(45,1), new LatLon(45,2), new LatLon(46,2), new LatLon(46,1) ]; - * const p = new LatLon(45.1, 1.1); - * const inside = p.isEnclosedBy(bounds); // true - */ - isEnclosedBy(polygon) { - // this method uses angle summation test; on a plane, angles for an enclosed point will sum - // to 360°, angles for an exterior point will sum to 0°. On a sphere, enclosed point angles - // will sum to less than 360° (due to spherical excess), exterior point angles will be small - // but non-zero. TODO: are any winding number optimisations applicable to spherical surface? - - if (!(polygon instanceof Array)) throw new TypeError(`isEnclosedBy: polygon must be Array (not ${classOf(polygon)})`); - if (!(polygon[0] instanceof LatLonNvectorSpherical)) throw new TypeError(`isEnclosedBy: polygon must be Array of LatLon (not ${classOf(polygon[0])})`); - if (polygon.length < 3) return false; // or throw? - - const nVertices = polygon.length; - - const p = this.toNvector(); - - // get vectors from p to each vertex - const vectorToVertex = []; - for (let v=0; v π; - } - - - /** - * Calculates the area of a spherical polygon where the sides of the polygon are great circle - * arcs joining the vertices. - * - * Uses Girard’s theorem: A = [Σθᵢ − (n−2)·π]·R² - * - * @param {LatLon[]} polygon - Array of points defining vertices of the polygon. - * @param {number} [radius=6371e3] - (Mean) radius of earth (defaults to radius in metres). - * @returns {number} The area of the polygon in the same units as radius. - * - * @example - * const polygon = [ new LatLon(0,0), new LatLon(1,0), new LatLon(0,1) ]; - * const area = LatLon.areaOf(polygon); // 6.18e9 m² - */ - static areaOf(polygon, radius=6371e3) { - const R = Number(radius); - - // get great-circle vector representing each segment - const c = []; - for (let v=0; v π/2) centreV = centreV.negate(); - - const centreP = new NvectorSpherical(centreV.x, centreV.y, centreV.z).toLatLon(); - - return centreP; - } - static centerOf(polygon) { return LatLonNvectorSpherical.centreOf(polygon); } // for en-us American English - - - /** - * 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. - * - * @example - * const p = LatLon.meanOf([ new LatLon(1, 1), new LatLon(4, 2), new LatLon(1, 3) ]); // 02.0001°N, 002.0000°E - */ - static meanOf(points) { - let m = new NvectorSpherical(0, 0, 0); // null vector - - // add all vectors - for (let 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. - * - * @param {LatLon} point - Point to be compared against this point. - * @returns {bool} True if points have identical latitude and longitude values. - * @throws {TypeError} Invalid point. - * - * @example - * const p1 = new LatLon(52.205, 0.119); - * const p2 = new LatLon(52.205, 0.119); - * const equal = p1.equals(p2); // true - */ - equals(point) { - if (!(point instanceof LatLonNvectorSpherical)) throw new TypeError(`invalid point ‘${point}’`); - - if (Math.abs(this.lat - point.lat) > Number.EPSILON) return false; - if (Math.abs(this.lon - point.lon) > Number.EPSILON) return false; - - return true; - } - - - /** - * Converts ‘this’ point to a GeoJSON object. - * - * @returns {Object} this point as a GeoJSON ‘Point’ object. - */ - toGeoJSON() { - return { type: 'Point', coordinates: [ this.lon, this.lat ] }; - } - - - /** - * Returns a string representation of ‘this’ point, formatted as degrees, degrees+minutes, or - * degrees+minutes+seconds. - * - * @param {string} [format=d] - Format point as 'd', 'dm', 'dms', or 'n' for signed numeric. - * @param {number} [dp=4|2|0] - Number of decimal places to use: default 4 for d, 2 for dm, 0 for dms. - * @returns {string} Comma-separated formatted latitude/longitude. - * - * @example - * const greenwich = new LatLon(51.47788, -0.00147); - * const d = greenwich.toString(); // 51.4778°N, 000.0015°W - * const dms = greenwich.toString('dms', 2); // 51°28′40.37″N, 000°00′05.29″W - * const [lat, lon] = greenwich.toString('n').split(','); // 51.4778, -0.0015 - */ - toString(format='d', dp=undefined) { - // note: explicitly set dp to undefined for passing through to toLat/toLon - if (![ 'd', 'dm', 'dms', 'n' ].includes(format)) throw new RangeError(`invalid format ‘${format}’`); - - if (format == 'n') { // signed numeric degrees - if (dp == undefined) dp = 4; - return `${this.lat.toFixed(dp)},${this.lon.toFixed(dp)}`; - } - const lat = Dms.toLat(this.lat, format, dp); - const lon = Dms.toLon(this.lon, format, dp); - return `${lat}, ${lon}`; - } - -} - - -/* Nvector - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - - -/** - * An n-vector is a (unit) 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 a (normalised) earth-centred earth-fixed - * (ECEF) vector. - * - * @extends Vector3d - */ -class NvectorSpherical extends Vector3d { - - // note commonality with latlon-nvector-ellipsoidal - - /** - * Creates a 3d n-vector normal to the Earth’s surface. - * - * @param {number} x - X component of n-vector (towards 0°N, 0°E). - * @param {number} y - Y component of n-vector (towards 0°N, 90°E). - * @param {number} z - Z component of n-vector (towards 90°N). - * - * @example - * import { Nvector } from '/js/geodesy/latlon-nvector-spherical.js'; - * const n = new Nvector(0.5000, 0.5000, 0.7071); - */ - constructor(x, y, z) { - const u = new Vector3d(x, y, z).unit(); // n-vectors are always normalised - - super(u.x, u.y, u.z); - } - - - /** - * Converts ‘this’ n-vector to latitude/longitude point. - * - * @returns {LatLon} Latitude/longitude point vector points to. - * - * @example - * const n = new Nvector(0.5000, 0.5000, 0.7071); - * const p = n.toLatLon(); // 45.0°N, 045.0°E - */ - toLatLon() { - // tanφ = z / √(x²+y²), tanλ = y / x (same as ellipsoidal calculation) - - const x = this.x, y = this.y, z = this.z; - - const φ = Math.atan2(z, Math.sqrt(x*x + y*y)); - const λ = Math.atan2(y, x); - - return new LatLonNvectorSpherical(φ.toDegrees(), λ.toDegrees()); - } - - - /** - * Vector normal to great circle obtained by heading on given bearing from point given by - * ‘this’ n-vector. - * - * Direction of vector is such that initial bearing vector b = c × n, where n is an n-vector - * representing ‘this’ (start) point. - * - * @private - * @param {number} bearing - Compass bearing in degrees. - * @returns {Vector3d} Normalised vector representing great circle. - * - * @example - * const n1 = new LatLon(53.3206, -1.7297).toNvector(); - * const gc = n1.greatCircle(96.0); // [-0.794,0.129,0.594] - */ - greatCircle(bearing) { - const θ = Number(bearing).toRadians(); - - const N = new Vector3d(0, 0, 1); // n-vector representing north pole - const e = N.cross(this); // easting - const n = this.cross(e); // northing - const eʹ = e.times(Math.cos(θ)/e.length); - const nʹ = n.times(Math.sin(θ)/n.length); - const c = nʹ.minus(eʹ); - - return c; - } - - - /** - * 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. - * - * @example - * const v = new Nvector(0.5000, 0.5000, 0.7071).toString(); // [0.500,0.500,0.707] - */ - toString(dp=3) { - const x = this.x.toFixed(dp); - const y = this.y.toFixed(dp); - const z = this.z.toFixed(dp); - - return `[${x},${y},${z}]`; - } - -} - - -/** - * Return class of supplied argument; javascriptweblog.wordpress.com/2011/08/08. - * - * @param {any} thing - Object whose class is to be determined. - * @returns {string} Class of supplied object. - */ -function classOf(thing) { - return ({}).toString.call(thing).match(/\s([a-zA-Z0-9]+)/)[1]; -} - - -/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - -export { LatLonNvectorSpherical as default, NvectorSpherical as Nvector, Dms }; diff --git a/src/js/geo/latlon-spherical.js b/src/js/geo/latlon-spherical.js deleted file mode 100644 index 67ba244..0000000 --- a/src/js/geo/latlon-spherical.js +++ /dev/null @@ -1,865 +0,0 @@ -/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ -/* Latitude/longitude spherical geodesy tools (c) Chris Veness 2002-2022 */ -/* MIT Licence */ -/* www.movable-type.co.uk/scripts/latlong.html */ -/* www.movable-type.co.uk/scripts/geodesy-library.html#latlon-spherical */ -/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - -import Dms from './dms.js'; - -const π = Math.PI; - - -/** - * Library of geodesy functions for operations on a spherical earth model. - * - * Includes distances, bearings, destinations, etc, for both great circle paths and rhumb lines, - * and other related functions. - * - * All calculations are done using simple spherical trigonometric formulae. - * - * @module latlon-spherical - */ - -// note greek letters (e.g. φ, λ, θ) are used for angles in radians to distinguish from angles in -// degrees (e.g. lat, lon, brng) - - -/* LatLonSpherical - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - - -/** - * Latitude/longitude points on a spherical model earth, and methods for calculating distances, - * bearings, destinations, etc on (orthodromic) great-circle paths and (loxodromic) rhumb lines. - */ -class LatLonSpherical { - - /** - * Creates a latitude/longitude point on the earth’s surface, using a spherical model earth. - * - * @param {number} lat - Latitude (in degrees). - * @param {number} lon - Longitude (in degrees). - * @throws {TypeError} Invalid lat/lon. - * - * @example - * import LatLon from '/js/geodesy/latlon-spherical.js'; - * const p = new LatLon(52.205, 0.119); - */ - constructor(lat, lon) { - if (isNaN(lat)) throw new TypeError(`invalid lat ‘${lat}’`); - if (isNaN(lon)) throw new TypeError(`invalid lon ‘${lon}’`); - - this._lat = Dms.wrap90(Number(lat)); - this._lon = Dms.wrap180(Number(lon)); - } - - - /** - * Latitude in degrees north from equator (including aliases lat, latitude): can be set as - * numeric or hexagesimal (deg-min-sec); returned as numeric. - */ - get lat() { return this._lat; } - get latitude() { return this._lat; } - set lat(lat) { - this._lat = isNaN(lat) ? Dms.wrap90(Dms.parse(lat)) : Dms.wrap90(Number(lat)); - if (isNaN(this._lat)) throw new TypeError(`invalid lat ‘${lat}’`); - } - set latitude(lat) { - this._lat = isNaN(lat) ? Dms.wrap90(Dms.parse(lat)) : Dms.wrap90(Number(lat)); - if (isNaN(this._lat)) throw new TypeError(`invalid latitude ‘${lat}’`); - } - - /** - * Longitude in degrees east from international reference meridian (including aliases lon, lng, - * longitude): can be set as numeric or hexagesimal (deg-min-sec); returned as numeric. - */ - get lon() { return this._lon; } - get lng() { return this._lon; } - get longitude() { return this._lon; } - set lon(lon) { - this._lon = isNaN(lon) ? Dms.wrap180(Dms.parse(lon)) : Dms.wrap180(Number(lon)); - if (isNaN(this._lon)) throw new TypeError(`invalid lon ‘${lon}’`); - } - set lng(lon) { - this._lon = isNaN(lon) ? Dms.wrap180(Dms.parse(lon)) : Dms.wrap180(Number(lon)); - if (isNaN(this._lon)) throw new TypeError(`invalid lng ‘${lon}’`); - } - set longitude(lon) { - this._lon = isNaN(lon) ? Dms.wrap180(Dms.parse(lon)) : Dms.wrap180(Number(lon)); - if (isNaN(this._lon)) throw new TypeError(`invalid longitude ‘${lon}’`); - } - - - /** Conversion factors; 1000 * LatLon.metresToKm gives 1. */ - static get metresToKm() { return 1/1000; } - /** Conversion factors; 1000 * LatLon.metresToMiles gives 0.621371192237334. */ - static get metresToMiles() { return 1/1609.344; } - /** Conversion factors; 1000 * LatLon.metresToMiles gives 0.5399568034557236. */ - static get metresToNauticalMiles() { return 1/1852; } - - - /** - * Parses a latitude/longitude point from a variety of formats. - * - * Latitude & longitude (in degrees) can be supplied as two separate parameters, as a single - * comma-separated lat/lon string, or as a single object with { lat, lon } or GeoJSON properties. - * - * The latitude/longitude values may be numeric or strings; they may be signed decimal or - * deg-min-sec (hexagesimal) suffixed by compass direction (NSEW); a variety of separators are - * accepted. Examples -3.62, '3 37 12W', '3°37′12″W'. - * - * Thousands/decimal separators must be comma/dot; use Dms.fromLocale to convert locale-specific - * thousands/decimal separators. - * - * @param {number|string|Object} lat|latlon - Latitude (in degrees) or comma-separated lat/lon or lat/lon object. - * @param {number|string} [lon] - Longitude (in degrees). - * @returns {LatLon} Latitude/longitude point. - * @throws {TypeError} Invalid point. - * - * @example - * const p1 = LatLon.parse(52.205, 0.119); // numeric pair (≡ new LatLon) - * const p2 = LatLon.parse('52.205', '0.119'); // numeric string pair (≡ new LatLon) - * const p3 = LatLon.parse('52.205, 0.119'); // single string numerics - * const p4 = LatLon.parse('52°12′18.0″N', '000°07′08.4″E'); // DMS pair - * const p5 = LatLon.parse('52°12′18.0″N, 000°07′08.4″E'); // single string DMS - * const p6 = LatLon.parse({ lat: 52.205, lon: 0.119 }); // { lat, lon } object numeric - * const p7 = LatLon.parse({ lat: '52°12′18.0″N', lng: '000°07′08.4″E' }); // { lat, lng } object DMS - * const p8 = LatLon.parse({ type: 'Point', coordinates: [ 0.119, 52.205] }); // GeoJSON - */ - static parse(...args) { - if (args.length == 0) throw new TypeError('invalid (empty) point'); - if (args[0]===null || args[1]===null) throw new TypeError('invalid (null) point'); - - let lat=undefined, lon=undefined; - - if (args.length == 2) { // regular (lat, lon) arguments - [ lat, lon ] = args; - lat = Dms.wrap90(Dms.parse(lat)); - lon = Dms.wrap180(Dms.parse(lon)); - if (isNaN(lat) || isNaN(lon)) throw new TypeError(`invalid point ‘${args.toString()}’`); - } - - if (args.length == 1 && typeof args[0] == 'string') { // single comma-separated lat,lon string - [ lat, lon ] = args[0].split(','); - lat = Dms.wrap90(Dms.parse(lat)); - lon = Dms.wrap180(Dms.parse(lon)); - if (isNaN(lat) || isNaN(lon)) throw new TypeError(`invalid point ‘${args[0]}’`); - } - - if (args.length == 1 && typeof args[0] == 'object') { // single { lat, lon } object - const ll = args[0]; - if (ll.type == 'Point' && Array.isArray(ll.coordinates)) { // GeoJSON - [ lon, lat ] = ll.coordinates; - } else { // regular { lat, lon } object - if (ll.latitude != undefined) lat = ll.latitude; - if (ll.lat != undefined) lat = ll.lat; - if (ll.longitude != undefined) lon = ll.longitude; - if (ll.lng != undefined) lon = ll.lng; - if (ll.lon != undefined) lon = ll.lon; - lat = Dms.wrap90(Dms.parse(lat)); - lon = Dms.wrap180(Dms.parse(lon)); - } - if (isNaN(lat) || isNaN(lon)) throw new TypeError(`invalid point ‘${JSON.stringify(args[0])}’`); - } - - if (isNaN(lat) || isNaN(lon)) throw new TypeError(`invalid point ‘${args.toString()}’`); - - return new LatLonSpherical(lat, lon); - } - - - /** - * Returns the distance along the surface of the earth from ‘this’ point to destination point. - * - * Uses haversine formula: a = sin²(Δφ/2) + cosφ1·cosφ2 · sin²(Δλ/2); d = 2 · atan2(√a, √(a-1)). - * - * @param {LatLon} point - Latitude/longitude of destination point. - * @param {number} [radius=6371e3] - Radius of earth (defaults to mean radius in metres). - * @returns {number} Distance between this point and destination point, in same units as radius. - * @throws {TypeError} Invalid radius. - * - * @example - * const p1 = new LatLon(52.205, 0.119); - * const p2 = new LatLon(48.857, 2.351); - * const d = p1.distanceTo(p2); // 404.3×10³ m - * const m = p1.distanceTo(p2, 3959); // 251.2 miles - */ - distanceTo(point, radius=6371e3) { - if (!(point instanceof LatLonSpherical)) point = LatLonSpherical.parse(point); // allow literal forms - if (isNaN(radius)) throw new TypeError(`invalid radius ‘${radius}’`); - - // a = sin²(Δφ/2) + cos(φ1)⋅cos(φ2)⋅sin²(Δλ/2) - // δ = 2·atan2(√(a), √(1−a)) - // see mathforum.org/library/drmath/view/51879.html for derivation - - const R = radius; - const φ1 = this.lat.toRadians(), λ1 = this.lon.toRadians(); - const φ2 = point.lat.toRadians(), λ2 = point.lon.toRadians(); - const Δφ = φ2 - φ1; - const Δλ = λ2 - λ1; - - const a = Math.sin(Δφ/2)*Math.sin(Δφ/2) + Math.cos(φ1)*Math.cos(φ2) * Math.sin(Δλ/2)*Math.sin(Δλ/2); - const c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a)); - const d = R * c; - - return d; - } - - - /** - * Returns the initial bearing from ‘this’ point to destination point. - * - * @param {LatLon} point - Latitude/longitude of destination point. - * @returns {number} Initial bearing in degrees from north (0°..360°). - * - * @example - * const p1 = new LatLon(52.205, 0.119); - * const p2 = new LatLon(48.857, 2.351); - * const b1 = p1.initialBearingTo(p2); // 156.2° - */ - initialBearingTo(point) { - if (!(point instanceof LatLonSpherical)) point = LatLonSpherical.parse(point); // allow literal forms - if (this.equals(point)) return NaN; // coincident points - - // tanθ = sinΔλ⋅cosφ2 / cosφ1⋅sinφ2 − sinφ1⋅cosφ2⋅cosΔλ - // see mathforum.org/library/drmath/view/55417.html for derivation - - const φ1 = this.lat.toRadians(); - const φ2 = point.lat.toRadians(); - const Δλ = (point.lon - this.lon).toRadians(); - - const x = Math.cos(φ1) * Math.sin(φ2) - Math.sin(φ1) * Math.cos(φ2) * Math.cos(Δλ); - const y = Math.sin(Δλ) * Math.cos(φ2); - const θ = Math.atan2(y, x); - - const bearing = θ.toDegrees(); - - return Dms.wrap360(bearing); - } - - - /** - * Returns final bearing arriving at destination point from ‘this’ point; the final bearing will - * differ from the initial bearing by varying degrees according to distance and latitude. - * - * @param {LatLon} point - Latitude/longitude of destination point. - * @returns {number} Final bearing in degrees from north (0°..360°). - * - * @example - * const p1 = new LatLon(52.205, 0.119); - * const p2 = new LatLon(48.857, 2.351); - * const b2 = p1.finalBearingTo(p2); // 157.9° - */ - finalBearingTo(point) { - if (!(point instanceof LatLonSpherical)) point = LatLonSpherical.parse(point); // allow literal forms - - // get initial bearing from destination point to this point & reverse it by adding 180° - - const bearing = point.initialBearingTo(this) + 180; - - return Dms.wrap360(bearing); - } - - - /** - * Returns the midpoint between ‘this’ point and destination point. - * - * @param {LatLon} point - Latitude/longitude of destination point. - * @returns {LatLon} Midpoint between this point and destination point. - * - * @example - * const p1 = new LatLon(52.205, 0.119); - * const p2 = new LatLon(48.857, 2.351); - * const pMid = p1.midpointTo(p2); // 50.5363°N, 001.2746°E - */ - midpointTo(point) { - if (!(point instanceof LatLonSpherical)) point = LatLonSpherical.parse(point); // allow literal forms - - // φm = atan2( sinφ1 + sinφ2, √( (cosφ1 + cosφ2⋅cosΔλ)² + cos²φ2⋅sin²Δλ ) ) - // λm = λ1 + atan2(cosφ2⋅sinΔλ, cosφ1 + cosφ2⋅cosΔλ) - // midpoint is sum of vectors to two points: mathforum.org/library/drmath/view/51822.html - - const φ1 = this.lat.toRadians(); - const λ1 = this.lon.toRadians(); - const φ2 = point.lat.toRadians(); - const Δλ = (point.lon - this.lon).toRadians(); - - // get cartesian coordinates for the two points - const A = { x: Math.cos(φ1), y: 0, z: Math.sin(φ1) }; // place point A on prime meridian y=0 - const B = { x: Math.cos(φ2)*Math.cos(Δλ), y: Math.cos(φ2)*Math.sin(Δλ), z: Math.sin(φ2) }; - - // vector to midpoint is sum of vectors to two points (no need to normalise) - const C = { x: A.x + B.x, y: A.y + B.y, z: A.z + B.z }; - - const φm = Math.atan2(C.z, Math.sqrt(C.x*C.x + C.y*C.y)); - const λm = λ1 + Math.atan2(C.y, C.x); - - const lat = φm.toDegrees(); - const lon = λm.toDegrees(); - - return new LatLonSpherical(lat, lon); - } - - - /** - * 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(52.205, 0.119); - * const p2 = new LatLon(48.857, 2.351); - * const pInt = p1.intermediatePointTo(p2, 0.25); // 51.3721°N, 000.7073°E - */ - intermediatePointTo(point, fraction) { - if (!(point instanceof LatLonSpherical)) point = LatLonSpherical.parse(point); // allow literal forms - if (this.equals(point)) return new LatLonSpherical(this.lat, this.lon); // coincident points - - const φ1 = this.lat.toRadians(), λ1 = this.lon.toRadians(); - const φ2 = point.lat.toRadians(), λ2 = point.lon.toRadians(); - - // distance between points - const Δφ = φ2 - φ1; - const Δλ = λ2 - λ1; - const a = Math.sin(Δφ/2) * Math.sin(Δφ/2) - + Math.cos(φ1) * Math.cos(φ2) * Math.sin(Δλ/2) * Math.sin(Δλ/2); - const δ = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a)); - - const A = Math.sin((1-fraction)*δ) / Math.sin(δ); - const B = Math.sin(fraction*δ) / Math.sin(δ); - - const x = A * Math.cos(φ1) * Math.cos(λ1) + B * Math.cos(φ2) * Math.cos(λ2); - const y = A * Math.cos(φ1) * Math.sin(λ1) + B * Math.cos(φ2) * Math.sin(λ2); - const z = A * Math.sin(φ1) + B * Math.sin(φ2); - - const φ3 = Math.atan2(z, Math.sqrt(x*x + y*y)); - const λ3 = Math.atan2(y, x); - - const lat = φ3.toDegrees(); - const lon = λ3.toDegrees(); - - return new LatLonSpherical(lat, lon); - } - - - /** - * Returns the destination point from ‘this’ point having travelled the given distance on the - * given initial bearing (bearing normally varies around path followed). - * - * @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. - * @throws {TypeError} Invalid distance/bearing/radius. - * - * @example - * const p1 = new LatLon(51.47788, -0.00147); - * const p2 = p1.destinationPoint(7794, 300.7); // 51.5136°N, 000.0983°W - */ - destinationPoint(distance, bearing, radius=6371e3) { - if (isNaN(distance)) throw new TypeError(`invalid distance ‘${distance}’`); - if (isNaN(bearing)) throw new TypeError(`invalid bearing ‘${bearing}’`); - if (isNaN(radius)) throw new TypeError(`invalid radius ‘${radius}’`); - - // sinφ2 = sinφ1⋅cosδ + cosφ1⋅sinδ⋅cosθ - // tanΔλ = sinθ⋅sinδ⋅cosφ1 / cosδ−sinφ1⋅sinφ2 - // see mathforum.org/library/drmath/view/52049.html for derivation - - const δ = distance / radius; // angular distance in radians - const θ = Number(bearing).toRadians(); - - const φ1 = this.lat.toRadians(), λ1 = this.lon.toRadians(); - - const sinφ2 = Math.sin(φ1) * Math.cos(δ) + Math.cos(φ1) * Math.sin(δ) * Math.cos(θ); - const φ2 = Math.asin(sinφ2); - const y = Math.sin(θ) * Math.sin(δ) * Math.cos(φ1); - const x = Math.cos(δ) - Math.sin(φ1) * sinφ2; - const λ2 = λ1 + Math.atan2(y, x); - - const lat = φ2.toDegrees(); - const lon = λ2.toDegrees(); - - return new LatLonSpherical(lat, lon); - } - - - /** - * Returns the point of intersection of two paths defined by point and bearing. - * - * @param {LatLon} p1 - First point. - * @param {number} brng1 - Initial bearing from first point. - * @param {LatLon} p2 - Second point. - * @param {number} brng2 - Initial bearing from second point. - * @returns {LatLon|null} Destination point (null if no unique intersection defined). - * - * @example - * const p1 = new LatLon(51.8853, 0.2545), brng1 = 108.547; - * const p2 = new LatLon(49.0034, 2.5735), brng2 = 32.435; - * const pInt = LatLon.intersection(p1, brng1, p2, brng2); // 50.9078°N, 004.5084°E - */ - static intersection(p1, brng1, p2, brng2) { - if (!(p1 instanceof LatLonSpherical)) p1 = LatLonSpherical.parse(p1); // allow literal forms - if (!(p2 instanceof LatLonSpherical)) p2 = LatLonSpherical.parse(p2); // allow literal forms - if (isNaN(brng1)) throw new TypeError(`invalid brng1 ‘${brng1}’`); - if (isNaN(brng2)) throw new TypeError(`invalid brng2 ‘${brng2}’`); - - // see www.edwilliams.org/avform.htm#Intersection - - const φ1 = p1.lat.toRadians(), λ1 = p1.lon.toRadians(); - const φ2 = p2.lat.toRadians(), λ2 = p2.lon.toRadians(); - const θ13 = Number(brng1).toRadians(), θ23 = Number(brng2).toRadians(); - const Δφ = φ2 - φ1, Δλ = λ2 - λ1; - - // angular distance p1-p2 - const δ12 = 2 * Math.asin(Math.sqrt(Math.sin(Δφ/2) * Math.sin(Δφ/2) - + Math.cos(φ1) * Math.cos(φ2) * Math.sin(Δλ/2) * Math.sin(Δλ/2))); - if (Math.abs(δ12) < Number.EPSILON) return new LatLonSpherical(p1.lat, p1.lon); // coincident points - - // initial/final bearings between points - const cosθa = (Math.sin(φ2) - Math.sin(φ1)*Math.cos(δ12)) / (Math.sin(δ12)*Math.cos(φ1)); - const cosθb = (Math.sin(φ1) - Math.sin(φ2)*Math.cos(δ12)) / (Math.sin(δ12)*Math.cos(φ2)); - const θa = Math.acos(Math.min(Math.max(cosθa, -1), 1)); // protect against rounding errors - const θb = Math.acos(Math.min(Math.max(cosθb, -1), 1)); // protect against rounding errors - - const θ12 = Math.sin(λ2-λ1)>0 ? θa : 2*π-θa; - const θ21 = Math.sin(λ2-λ1)>0 ? 2*π-θb : θb; - - const α1 = θ13 - θ12; // angle 2-1-3 - const α2 = θ21 - θ23; // angle 1-2-3 - - if (Math.sin(α1) == 0 && Math.sin(α2) == 0) return null; // infinite intersections - if (Math.sin(α1) * Math.sin(α2) < 0) return null; // ambiguous intersection (antipodal/360°) - - const cosα3 = -Math.cos(α1)*Math.cos(α2) + Math.sin(α1)*Math.sin(α2)*Math.cos(δ12); - - const δ13 = Math.atan2(Math.sin(δ12)*Math.sin(α1)*Math.sin(α2), Math.cos(α2) + Math.cos(α1)*cosα3); - - const φ3 = Math.asin(Math.min(Math.max(Math.sin(φ1)*Math.cos(δ13) + Math.cos(φ1)*Math.sin(δ13)*Math.cos(θ13), -1), 1)); - - const Δλ13 = Math.atan2(Math.sin(θ13)*Math.sin(δ13)*Math.cos(φ1), Math.cos(δ13) - Math.sin(φ1)*Math.sin(φ3)); - const λ3 = λ1 + Δλ13; - - const lat = φ3.toDegrees(); - const lon = λ3.toDegrees(); - - return new LatLonSpherical(lat, lon); - } - - - /** - * Returns (signed) distance from ‘this’ point to great circle defined by start-point and - * end-point. - * - * @param {LatLon} pathStart - Start point of great circle path. - * @param {LatLon} pathEnd - End point of great circle path. - * @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 - * const pCurrent = new LatLon(53.2611, -0.7972); - * const p1 = new LatLon(53.3206, -1.7297); - * const p2 = new LatLon(53.1887, 0.1334); - * const d = pCurrent.crossTrackDistanceTo(p1, p2); // -307.5 m - */ - crossTrackDistanceTo(pathStart, pathEnd, radius=6371e3) { - if (!(pathStart instanceof LatLonSpherical)) pathStart = LatLonSpherical.parse(pathStart); // allow literal forms - if (!(pathEnd instanceof LatLonSpherical)) pathEnd = LatLonSpherical.parse(pathEnd); // allow literal forms - const R = radius; - - if (this.equals(pathStart)) return 0; - - const δ13 = pathStart.distanceTo(this, R) / R; - const θ13 = pathStart.initialBearingTo(this).toRadians(); - const θ12 = pathStart.initialBearingTo(pathEnd).toRadians(); - - const δxt = Math.asin(Math.sin(δ13) * Math.sin(θ13 - θ12)); - - return δxt * R; - } - - - /** - * Returns how far ‘this’ point is along a path from from start-point, heading towards end-point. - * That is, if a perpendicular is drawn from ‘this’ point to the (great circle) path, the - * along-track distance is the distance from the start point to where the perpendicular crosses - * the path. - * - * @param {LatLon} pathStart - Start point of great circle path. - * @param {LatLon} pathEnd - End point of great circle path. - * @param {number} [radius=6371e3] - (Mean) radius of earth (defaults to radius in metres). - * @returns {number} Distance along great circle to point nearest ‘this’ point. - * - * @example - * const pCurrent = new LatLon(53.2611, -0.7972); - * const p1 = new LatLon(53.3206, -1.7297); - * const p2 = new LatLon(53.1887, 0.1334); - * const d = pCurrent.alongTrackDistanceTo(p1, p2); // 62.331 km - */ - alongTrackDistanceTo(pathStart, pathEnd, radius=6371e3) { - if (!(pathStart instanceof LatLonSpherical)) pathStart = LatLonSpherical.parse(pathStart); // allow literal forms - if (!(pathEnd instanceof LatLonSpherical)) pathEnd = LatLonSpherical.parse(pathEnd); // allow literal forms - const R = radius; - - if (this.equals(pathStart)) return 0; - - const δ13 = pathStart.distanceTo(this, R) / R; - const θ13 = pathStart.initialBearingTo(this).toRadians(); - const θ12 = pathStart.initialBearingTo(pathEnd).toRadians(); - - const δxt = Math.asin(Math.sin(δ13) * Math.sin(θ13-θ12)); - - const δat = Math.acos(Math.cos(δ13) / Math.abs(Math.cos(δxt))); - - return δat*Math.sign(Math.cos(θ12-θ13)) * R; - } - - - /** - * Returns maximum latitude reached when travelling on a great circle on given bearing from - * ‘this’ point (‘Clairaut’s formula’). Negate the result for the minimum latitude (in the - * southern hemisphere). - * - * The maximum latitude is independent of longitude; it will be the same for all points on a - * given latitude. - * - * @param {number} bearing - Initial bearing. - * @returns {number} Maximum latitude reached. - */ - maxLatitude(bearing) { - const θ = Number(bearing).toRadians(); - - const φ = this.lat.toRadians(); - - const φMax = Math.acos(Math.abs(Math.sin(θ) * Math.cos(φ))); - - return φMax.toDegrees(); - } - - - /** - * Returns the pair of meridians at which a great circle defined by two points crosses the given - * latitude. If the great circle doesn't reach the given latitude, null is returned. - * - * @param {LatLon} point1 - First point defining great circle. - * @param {LatLon} point2 - Second point defining great circle. - * @param {number} latitude - Latitude crossings are to be determined for. - * @returns {Object|null} Object containing { lon1, lon2 } or null if given latitude not reached. - */ - static crossingParallels(point1, point2, latitude) { - if (point1.equals(point2)) return null; // coincident points - - const φ = Number(latitude).toRadians(); - - const φ1 = point1.lat.toRadians(); - const λ1 = point1.lon.toRadians(); - const φ2 = point2.lat.toRadians(); - const λ2 = point2.lon.toRadians(); - - const Δλ = λ2 - λ1; - - const x = Math.sin(φ1) * Math.cos(φ2) * Math.cos(φ) * Math.sin(Δλ); - const y = Math.sin(φ1) * Math.cos(φ2) * Math.cos(φ) * Math.cos(Δλ) - Math.cos(φ1) * Math.sin(φ2) * Math.cos(φ); - const z = Math.cos(φ1) * Math.cos(φ2) * Math.sin(φ) * Math.sin(Δλ); - - if (z * z > x * x + y * y) return null; // great circle doesn't reach latitude - - const λm = Math.atan2(-y, x); // longitude at max latitude - const Δλi = Math.acos(z / Math.sqrt(x*x + y*y)); // Δλ from λm to intersection points - - const λi1 = λ1 + λm - Δλi; - const λi2 = λ1 + λm + Δλi; - - const lon1 = λi1.toDegrees(); - const lon2 = λi2.toDegrees(); - - return { - lon1: Dms.wrap180(lon1), - lon2: Dms.wrap180(lon2), - }; - } - - - /* Rhumb - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - - - /** - * Returns the distance travelling from ‘this’ point to destination point along a rhumb line. - * - * @param {LatLon} point - Latitude/longitude of destination point. - * @param {number} [radius=6371e3] - (Mean) radius of earth (defaults to radius in metres). - * @returns {number} Distance in km between this point and destination point (same units as radius). - * - * @example - * const p1 = new LatLon(51.127, 1.338); - * const p2 = new LatLon(50.964, 1.853); - * const d = p1.distanceTo(p2); // 40.31 km - */ - rhumbDistanceTo(point, radius=6371e3) { - if (!(point instanceof LatLonSpherical)) point = LatLonSpherical.parse(point); // allow literal forms - - // see www.edwilliams.org/avform.htm#Rhumb - - const R = radius; - const φ1 = this.lat.toRadians(); - const φ2 = point.lat.toRadians(); - const Δφ = φ2 - φ1; - let Δλ = Math.abs(point.lon - this.lon).toRadians(); - // if dLon over 180° take shorter rhumb line across the anti-meridian: - if (Math.abs(Δλ) > π) Δλ = Δλ > 0 ? -(2 * π - Δλ) : (2 * π + Δλ); - - // on Mercator projection, longitude distances shrink by latitude; q is the 'stretch factor' - // q becomes ill-conditioned along E-W line (0/0); use empirical tolerance to avoid it (note ε is too small) - const Δψ = Math.log(Math.tan(φ2 / 2 + π / 4) / Math.tan(φ1 / 2 + π / 4)); - const q = Math.abs(Δψ) > 10e-12 ? Δφ / Δψ : Math.cos(φ1); - - // distance is pythagoras on 'stretched' Mercator projection, √(Δφ² + q²·Δλ²) - const δ = Math.sqrt(Δφ*Δφ + q*q * Δλ*Δλ); // angular distance in radians - const d = δ * R; - - return d; - } - - - /** - * Returns the bearing from ‘this’ point to destination point along a rhumb line. - * - * @param {LatLon} point - Latitude/longitude of destination point. - * @returns {number} Bearing in degrees from north. - * - * @example - * const p1 = new LatLon(51.127, 1.338); - * const p2 = new LatLon(50.964, 1.853); - * const d = p1.rhumbBearingTo(p2); // 116.7° - */ - rhumbBearingTo(point) { - if (!(point instanceof LatLonSpherical)) point = LatLonSpherical.parse(point); // allow literal forms - if (this.equals(point)) return NaN; // coincident points - - const φ1 = this.lat.toRadians(); - const φ2 = point.lat.toRadians(); - let Δλ = (point.lon - this.lon).toRadians(); - // if dLon over 180° take shorter rhumb line across the anti-meridian: - if (Math.abs(Δλ) > π) Δλ = Δλ > 0 ? -(2 * π - Δλ) : (2 * π + Δλ); - - const Δψ = Math.log(Math.tan(φ2 / 2 + π / 4) / Math.tan(φ1 / 2 + π / 4)); - - const θ = Math.atan2(Δλ, Δψ); - - const bearing = θ.toDegrees(); - - return Dms.wrap360(bearing); - } - - - /** - * Returns the destination point having travelled along a rhumb line from ‘this’ point the given - * distance on the given bearing. - * - * @param {number} distance - Distance travelled, in same units as earth radius (default: metres). - * @param {number} bearing - Bearing in degrees from north. - * @param {number} [radius=6371e3] - (Mean) radius of earth (defaults to radius in metres). - * @returns {LatLon} Destination point. - * - * @example - * const p1 = new LatLon(51.127, 1.338); - * const p2 = p1.rhumbDestinationPoint(40300, 116.7); // 50.9642°N, 001.8530°E - */ - rhumbDestinationPoint(distance, bearing, radius=6371e3) { - const φ1 = this.lat.toRadians(), λ1 = this.lon.toRadians(); - const θ = Number(bearing).toRadians(); - - const δ = distance / radius; // angular distance in radians - - const Δφ = δ * Math.cos(θ); - let φ2 = φ1 + Δφ; - - // check for some daft bugger going past the pole, normalise latitude if so - if (Math.abs(φ2) > π / 2) φ2 = φ2 > 0 ? π - φ2 : -π - φ2; - - const Δψ = Math.log(Math.tan(φ2 / 2 + π / 4) / Math.tan(φ1 / 2 + π / 4)); - const q = Math.abs(Δψ) > 10e-12 ? Δφ / Δψ : Math.cos(φ1); // E-W course becomes ill-conditioned with 0/0 - - const Δλ = δ * Math.sin(θ) / q; - const λ2 = λ1 + Δλ; - - const lat = φ2.toDegrees(); - const lon = λ2.toDegrees(); - - return new LatLonSpherical(lat, lon); - } - - - /** - * Returns the loxodromic midpoint (along a rhumb line) between ‘this’ point and second point. - * - * @param {LatLon} point - Latitude/longitude of second point. - * @returns {LatLon} Midpoint between this point and second point. - * - * @example - * const p1 = new LatLon(51.127, 1.338); - * const p2 = new LatLon(50.964, 1.853); - * const pMid = p1.rhumbMidpointTo(p2); // 51.0455°N, 001.5957°E - */ - rhumbMidpointTo(point) { - if (!(point instanceof LatLonSpherical)) point = LatLonSpherical.parse(point); // allow literal forms - - // see mathforum.org/kb/message.jspa?messageID=148837 - - const φ1 = this.lat.toRadians(); let λ1 = this.lon.toRadians(); - const φ2 = point.lat.toRadians(), λ2 = point.lon.toRadians(); - - if (Math.abs(λ2 - λ1) > π) λ1 += 2 * π; // crossing anti-meridian - - const φ3 = (φ1 + φ2) / 2; - const f1 = Math.tan(π / 4 + φ1 / 2); - const f2 = Math.tan(π / 4 + φ2 / 2); - const f3 = Math.tan(π / 4 + φ3 / 2); - let λ3 = ((λ2 - λ1) * Math.log(f3) + λ1 * Math.log(f2) - λ2 * Math.log(f1)) / Math.log(f2 / f1); - - if (!isFinite(λ3)) λ3 = (λ1 + λ2) / 2; // parallel of latitude - - const lat = φ3.toDegrees(); - const lon = λ3.toDegrees(); - - return new LatLonSpherical(lat, lon); - } - - - /* Area - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - - - /** - * Calculates the area of a spherical polygon where the sides of the polygon are great circle - * arcs joining the vertices. - * - * @param {LatLon[]} polygon - Array of points defining vertices of the polygon. - * @param {number} [radius=6371e3] - (Mean) radius of earth (defaults to radius in metres). - * @returns {number} The area of the polygon in the same units as radius. - * - * @example - * const polygon = [new LatLon(0,0), new LatLon(1,0), new LatLon(0,1)]; - * const area = LatLon.areaOf(polygon); // 6.18e9 m² - */ - static areaOf(polygon, radius=6371e3) { - // uses method due to Karney: osgeo-org.1560.x6.nabble.com/Area-of-a-spherical-polygon-td3841625.html; - // for each edge of the polygon, tan(E/2) = tan(Δλ/2)·(tan(φ₁/2)+tan(φ₂/2)) / (1+tan(φ₁/2)·tan(φ₂/2)) - // where E is the spherical excess of the trapezium obtained by extending the edge to the equator - // (Karney's method is probably more efficient than the more widely known L’Huilier’s Theorem) - - const R = radius; - - // close polygon so that last point equals first point - const closed = polygon[0].equals(polygon[polygon.length-1]); - if (!closed) polygon.push(polygon[0]); - - const nVertices = polygon.length - 1; - - let S = 0; // spherical excess in steradians - for (let v=0; v Number.EPSILON) return false; - if (Math.abs(this.lon - point.lon) > Number.EPSILON) return false; - - return true; - } - - - /** - * Converts ‘this’ point to a GeoJSON object. - * - * @returns {Object} this point as a GeoJSON ‘Point’ object. - */ - toGeoJSON() { - return { type: 'Point', coordinates: [ this.lon, this.lat ] }; - } - - - /** - * Returns a string representation of ‘this’ point, formatted as degrees, degrees+minutes, or - * degrees+minutes+seconds. - * - * @param {string} [format=d] - Format point as 'd', 'dm', 'dms', or 'n' for signed numeric. - * @param {number} [dp=4|2|0] - Number of decimal places to use: default 4 for d, 2 for dm, 0 for dms. - * @returns {string} Comma-separated formatted latitude/longitude. - * @throws {RangeError} Invalid format. - * - * @example - * const greenwich = new LatLon(51.47788, -0.00147); - * const d = greenwich.toString(); // 51.4779°N, 000.0015°W - * const dms = greenwich.toString('dms', 2); // 51°28′40.37″N, 000°00′05.29″W - * const [lat, lon] = greenwich.toString('n').split(','); // 51.4779, -0.0015 - */ - toString(format='d', dp=undefined) { - // note: explicitly set dp to undefined for passing through to toLat/toLon - if (![ 'd', 'dm', 'dms', 'n' ].includes(format)) throw new RangeError(`invalid format ‘${format}’`); - - if (format == 'n') { // signed numeric degrees - if (dp == undefined) dp = 4; - return `${this.lat.toFixed(dp)},${this.lon.toFixed(dp)}`; - } - const lat = Dms.toLat(this.lat, format, dp); - const lon = Dms.toLon(this.lon, format, dp); - return `${lat}, ${lon}`; - } - -} - - -/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - -export { LatLonSpherical as default, Dms }; diff --git a/src/js/geo/mgrs.js b/src/js/geo/mgrs.js deleted file mode 100644 index a0030f5..0000000 --- a/src/js/geo/mgrs.js +++ /dev/null @@ -1,305 +0,0 @@ -/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ -/* MGRS / UTM Conversion Functions (c) Chris Veness 2014-2022 */ -/* MIT Licence */ -/* www.movable-type.co.uk/scripts/latlong-utm-mgrs.html */ -/* www.movable-type.co.uk/scripts/geodesy-library.html#mgrs */ -/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - -import Utm, { LatLon as LatLonEllipsoidal, Dms } from './utm.js'; - - -/** - * Military Grid Reference System (MGRS/NATO) grid references provides geocoordinate references - * covering the entire globe, based on UTM projections. - * - * MGRS references comprise a grid zone designator, a 100km square identification, and an easting - * and northing (in metres); e.g. ‘31U DQ 48251 11932’. - * - * Depending on requirements, some parts of the reference may be omitted (implied), and - * eastings/northings may be given to varying resolution. - * - * qv www.fgdc.gov/standards/projects/FGDC-standards-projects/usng/fgdc_std_011_2001_usng.pdf - * - * @module mgrs - */ - - -/* - * Latitude bands C..X 8° each, covering 80°S to 84°N - */ -const latBands = 'CDEFGHJKLMNPQRSTUVWXX'; // X is repeated for 80-84°N - - -/* - * 100km grid square column (‘e’) letters repeat every third zone - */ -const e100kLetters = [ 'ABCDEFGH', 'JKLMNPQR', 'STUVWXYZ' ]; - - -/* - * 100km grid square row (‘n’) letters repeat every other zone - */ -const n100kLetters = [ 'ABCDEFGHJKLMNPQRSTUV', 'FGHJKLMNPQRSTUVABCDE' ]; - - -/* Mgrs - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - - -/** - * Military Grid Reference System (MGRS/NATO) grid references, with methods to parse references, and - * to convert to UTM coordinates. - */ -class Mgrs { - - /** - * Creates an Mgrs grid reference object. - * - * @param {number} zone - 6° longitudinal zone (1..60 covering 180°W..180°E). - * @param {string} band - 8° latitudinal band (C..X covering 80°S..84°N). - * @param {string} e100k - First letter (E) of 100km grid square. - * @param {string} n100k - Second letter (N) of 100km grid square. - * @param {number} easting - Easting in metres within 100km grid square. - * @param {number} northing - Northing in metres within 100km grid square. - * @param {LatLon.datums} [datum=WGS84] - Datum UTM coordinate is based on. - * @throws {RangeError} Invalid MGRS grid reference. - * - * @example - * import Mgrs from '/js/geodesy/mgrs.js'; - * const mgrsRef = new Mgrs(31, 'U', 'D', 'Q', 48251, 11932); // 31U DQ 48251 11932 - */ - constructor(zone, band, e100k, n100k, easting, northing, datum=LatLonEllipsoidal.datums.WGS84) { - if (!(1<=zone && zone<=60)) throw new RangeError(`invalid MGRS zone ‘${zone}’`); - if (zone != parseInt(zone)) throw new RangeError(`invalid MGRS zone ‘${zone}’`); - const errors = []; // check & report all other possible errors rather than reporting one-by-one - if (band.length!=1 || latBands.indexOf(band) == -1) errors.push(`invalid MGRS band ‘${band}’`); - if (e100k.length!=1 || e100kLetters[(zone-1)%3].indexOf(e100k) == -1) errors.push(`invalid MGRS 100km grid square column ‘${e100k}’ for zone ${zone}`); - if (n100k.length!=1 || n100kLetters[0].indexOf(n100k) == -1) errors.push(`invalid MGRS 100km grid square row ‘${n100k}’`); - if (isNaN(Number(easting))) errors.push(`invalid MGRS easting ‘${easting}’`); - if (isNaN(Number(northing))) errors.push(`invalid MGRS northing ‘${northing}’`); - if (Number(easting) < 0 || Number(easting) > 99999) errors.push(`invalid MGRS easting ‘${easting}’`); - if (Number(northing) < 0 || Number(northing) > 99999) errors.push(`invalid MGRS northing ‘${northing}’`); - if (!datum || datum.ellipsoid==undefined) errors.push(`unrecognised datum ‘${datum}’`); - if (errors.length > 0) throw new RangeError(errors.join(', ')); - - this.zone = Number(zone); - this.band = band; - this.e100k = e100k; - this.n100k = n100k; - this.easting = Math.floor(easting); - this.northing = Math.floor(northing); - this.datum = datum; - } - - - /** - * Converts MGRS grid reference to UTM coordinate. - * - * Grid references refer to squares rather than points (with the size of the square indicated - * by the precision of the reference); this conversion will return the UTM coordinate of the SW - * corner of the grid reference square. - * - * @returns {Utm} UTM coordinate of SW corner of this MGRS grid reference. - * - * @example - * const mgrsRef = Mgrs.parse('31U DQ 48251 11932'); - * const utmCoord = mgrsRef.toUtm(); // 31 N 448251 5411932 - */ - toUtm() { - const hemisphere = this.band>='N' ? 'N' : 'S'; - - // get easting specified by e100k (note +1 because eastings start at 166e3 due to 500km false origin) - const col = e100kLetters[(this.zone-1)%3].indexOf(this.e100k) + 1; - const e100kNum = col * 100e3; // e100k in metres - - // get northing specified by n100k - const row = n100kLetters[(this.zone-1)%2].indexOf(this.n100k); - const n100kNum = row * 100e3; // n100k in metres - - // get latitude of (bottom of) band (10 bands above the equator, 8°latitude each) - const latBand = (latBands.indexOf(this.band)-10)*8; - - // get southern-most northing of bottom of band, using floor() to extend to include entirety - // of bottom-most 100km square - note in northern hemisphere, centre of zone will be furthest - // south; in southern hemisphere extremity of zone will be furthest south, so use 3°E / 0°E - const lon = this.band >= 'N' ? 3 : 0; - const nBand = Math.floor(new LatLonEllipsoidal(latBand, lon).toUtm().northing/100e3)*100e3; - - // 100km grid square row letters repeat every 2,000km north; add enough 2,000km blocks to - // get into required band - let n2M = 0; // northing of 2,000km block - while (n2M + n100kNum + this.northing < nBand) n2M += 2000e3; - - return new Utm_Mgrs(this.zone, hemisphere, e100kNum+this.easting, n2M+n100kNum+this.northing, this.datum); - } - - - /** - * Parses string representation of MGRS grid reference. - * - * An MGRS grid reference comprises (space-separated) - * - grid zone designator (GZD) - * - 100km grid square letter-pair - * - easting - * - northing. - * - * @param {string} mgrsGridRef - String representation of MGRS grid reference. - * @returns {Mgrs} Mgrs grid reference object. - * @throws {Error} Invalid MGRS grid reference. - * - * @example - * const mgrsRef = Mgrs.parse('31U DQ 48251 11932'); - * const mgrsRef = Mgrs.parse('31UDQ4825111932'); // military style no separators - * // mgrsRef: { zone:31, band:'U', e100k:'D', n100k:'Q', easting:48251, northing:11932 } - */ - static parse(mgrsGridRef) { - if (!mgrsGridRef) throw new Error(`invalid MGRS grid reference ‘${mgrsGridRef}’`); - - // check for military-style grid reference with no separators - if (!mgrsGridRef.trim().match(/\s/)) { // convert mgrsGridRef to standard space-separated format - const ref = mgrsGridRef.match(/(\d\d?[A-Z])([A-Z]{2})([0-9]{2,10})/i); - if (!ref) throw new Error(`invalid MGRS grid reference ‘${mgrsGridRef}’`); - - const [ , gzd, en100k, en ] = ref; // split grid ref into gzd, en100k, en - const [ easting, northing ] = [ en.slice(0, en.length/2), en.slice(-en.length/2) ]; - mgrsGridRef = `${gzd} ${en100k} ${easting} ${northing}`; - } - - // match separate elements (separated by whitespace) - const ref = mgrsGridRef.match(/\S+/g); // returns [ gzd, e100k, easting, northing ] - if (ref==null || ref.length!=4) throw new Error(`invalid MGRS grid reference ‘${mgrsGridRef}’`); - - const [ gzd, en100k, e, n ] = ref; // split grid ref into gzd, en100k, e, n - const [ , zone, band ] = gzd.match(/(\d\d?)([A-Z])/i); // split gzd into zone, band - const [ e100k, n100k ] = en100k.split(''); // split 100km letter-pair into e, n - - // standardise to 10-digit refs - ie metres) (but only if < 10-digit refs, to allow decimals) - const easting = e.length>=5 ? e : e.padEnd(5, '0'); - const northing = n.length>=5 ? n : n.padEnd(5, '0'); - - return new Mgrs(zone, band, e100k, n100k, easting, northing); - } - - - /** - * Returns a string representation of an MGRS grid reference. - * - * To distinguish from civilian UTM coordinate representations, no space is included within the - * zone/band grid zone designator. - * - * Components are separated by spaces: for a military-style unseparated string, use - * Mgrs.toString().replace(/ /g, ''); - * - * Note that MGRS grid references get truncated, not rounded (unlike UTM coordinates); grid - * references indicate a bounding square, rather than a point, with the size of the square - * indicated by the precision - a precision of 10 indicates a 1-metre square, a precision of 4 - * indicates a 1,000-metre square (hence 31U DQ 48 11 indicates a 1km square with SW corner at - * 31 N 448000 5411000, which would include the 1m square 31U DQ 48251 11932). - * - * @param {number} [digits=10] - Precision of returned grid reference (eg 4 = km, 10 = m). - * @returns {string} This grid reference in standard format. - * @throws {RangeError} Invalid precision. - * - * @example - * const mgrsStr = new Mgrs(31, 'U', 'D', 'Q', 48251, 11932).toString(); // 31U DQ 48251 11932 - */ - toString(digits=10) { - if (![ 2, 4, 6, 8, 10 ].includes(Number(digits))) throw new RangeError(`invalid precision ‘${digits}’`); - - const { zone, band, e100k, n100k, easting, northing } = this; - - // truncate to required precision - const eRounded = Math.floor(easting/Math.pow(10, 5-digits/2)); - const nRounded = Math.floor(northing/Math.pow(10, 5-digits/2)); - - // ensure leading zeros - const zPadded = zone.toString().padStart(2, '0'); - const ePadded = eRounded.toString().padStart(digits/2, '0'); - const nPadded = nRounded.toString().padStart(digits/2, '0'); - - return `${zPadded}${band} ${e100k}${n100k} ${ePadded} ${nPadded}`; - } -} - - -/* Utm_Mgrs - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - - -/** - * Extends Utm with method to convert UTM coordinate to MGRS reference. - * - * @extends Utm - */ -class Utm_Mgrs extends Utm { - - /** - * Converts UTM coordinate to MGRS reference. - * - * @returns {Mgrs} - * @throws {TypeError} Invalid UTM coordinate. - * - * @example - * const utmCoord = new Utm(31, 'N', 448251, 5411932); - * const mgrsRef = utmCoord.toMgrs(); // 31U DQ 48251 11932 - */ - toMgrs() { - // MGRS zone is same as UTM zone - const zone = this.zone; - - // convert UTM to lat/long to get latitude to determine band - const latlong = this.toLatLon(); - // grid zones are 8° tall, 0°N is 10th band - const band = latBands.charAt(Math.floor(latlong.lat.toFixed(12)/8+10)); // latitude band - - // columns in zone 1 are A-H, zone 2 J-R, zone 3 S-Z, then repeating every 3rd zone - const col = Math.floor(this.easting / 100e3); - // (note -1 because eastings start at 166e3 due to 500km false origin) - const e100k = e100kLetters[(zone-1)%3].charAt(col-1); - - // rows in even zones are A-V, in odd zones are F-E - const row = Math.floor(this.northing / 100e3) % 20; - const n100k = n100kLetters[(zone-1)%2].charAt(row); - - // truncate easting/northing to within 100km grid square & round to 1-metre precision - const easting = Math.floor(this.easting % 100e3); - const northing = Math.floor(this.northing % 100e3); - - return new Mgrs(zone, band, e100k, n100k, easting, northing); - } - -} - - -/** - * Extends LatLonEllipsoidal adding toMgrs() method to the Utm object returned by LatLon.toUtm(). - * - * @extends LatLonEllipsoidal - */ -class Latlon_Utm_Mgrs extends LatLonEllipsoidal { - - /** - * Converts latitude/longitude to UTM coordinate. - * - * Shadow of LatLon.toUtm, returning Utm augmented with toMgrs() method. - * - * @param {number} [zoneOverride] - Use specified zone rather than zone within which point lies; - * note overriding the UTM zone has the potential to result in negative eastings, and - * perverse results within Norway/Svalbard exceptions (this is unlikely to be relevant - * for MGRS, but is needed as Mgrs passes through the Utm class). - * @returns {Utm} UTM coordinate. - * @throws {Error} If point not valid, if point outside latitude range. - * - * @example - * const latlong = new LatLon(48.8582, 2.2945); - * const utmCoord = latlong.toUtm(); // 31 N 448252 5411933 - */ - toUtm(zoneOverride=undefined) { - const utm = super.toUtm(zoneOverride); - return new Utm_Mgrs(utm.zone, utm.hemisphere, utm.easting, utm.northing, utm.datum, utm.convergence, utm.scale); - } - -} - - -/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - -export { Mgrs as default, Utm_Mgrs as Utm, Latlon_Utm_Mgrs as LatLon, Dms }; diff --git a/src/js/geo/osgridref.js b/src/js/geo/osgridref.js deleted file mode 100644 index f73c24e..0000000 --- a/src/js/geo/osgridref.js +++ /dev/null @@ -1,348 +0,0 @@ -/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ -/* Ordnance Survey Grid Reference functions (c) Chris Veness 2005-2021 */ -/* MIT Licence */ -/* www.movable-type.co.uk/scripts/latlong-gridref.html */ -/* www.movable-type.co.uk/scripts/geodesy-library.html#osgridref */ -/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - -import LatLonEllipsoidal, { Dms } from './latlon-ellipsoidal-datum.js'; - - -/** - * Ordnance Survey OSGB grid references provide geocoordinate references for UK mapping purposes. - * - * Formulation implemented here due to Thomas, Redfearn, etc is as published by OS, but is inferior - * to Krüger as used by e.g. Karney 2011. - * - * www.ordnancesurvey.co.uk/documents/resources/guide-coordinate-systems-great-britain.pdf. - * - * Note OSGB grid references cover Great Britain only; Ireland and the Channel Islands have their - * own references. - * - * Note that these formulae are based on ellipsoidal calculations, and according to the OS are - * accurate to about 4–5 metres – for greater accuracy, a geoid-based transformation (OSTN15) must - * be used. - */ - -/* - * Converted 2015 to work with WGS84 by default, OSGB36 as option; - * www.ordnancesurvey.co.uk/blog/2014/12/confirmation-on-changes-to-latitude-and-longitude - */ - - -/* OsGridRef - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - - -const nationalGrid = { - trueOrigin: { lat: 49, lon: -2 }, // true origin of grid 49°N,2°W on OSGB36 datum - falseOrigin: { easting: -400e3, northing: 100e3 }, // easting & northing of false origin, metres from true origin - scaleFactor: 0.9996012717, // scale factor on central meridian - ellipsoid: LatLonEllipsoidal.ellipsoids.Airy1830, -}; -// note Irish National Grid uses t/o 53°30′N, 8°W, f/o 200kmW, 250kmS, scale factor 1.000035, on Airy 1830 Modified ellipsoid - - -/** - * OS Grid References with methods to parse and convert them to latitude/longitude points. - */ -class OsGridRef { - - /** - * Creates an OsGridRef object. - * - * @param {number} easting - Easting in metres from OS Grid false origin. - * @param {number} northing - Northing in metres from OS Grid false origin. - * - * @example - * import OsGridRef from '/js/geodesy/osgridref.js'; - * const gridref = new OsGridRef(651409, 313177); - */ - constructor(easting, northing) { - this.easting = Number(easting); - this.northing = Number(northing); - - if (isNaN(easting) || this.easting<0 || this.easting>700e3) throw new RangeError(`invalid easting ‘${easting}’`); - if (isNaN(northing) || this.northing<0 || this.northing>1300e3) throw new RangeError(`invalid northing ‘${northing}’`); - } - - - /** - * Converts ‘this’ Ordnance Survey Grid Reference easting/northing coordinate to latitude/longitude - * (SW corner of grid square). - * - * While OS Grid References are based on OSGB-36, the Ordnance Survey have deprecated the use of - * OSGB-36 for latitude/longitude coordinates (in favour of WGS-84), hence this function returns - * WGS-84 by default, with OSGB-36 as an option. See www.ordnancesurvey.co.uk/blog/2014/12/2. - * - * Note formulation implemented here due to Thomas, Redfearn, etc is as published by OS, but is - * inferior to Krüger as used by e.g. Karney 2011. - * - * @param {LatLon.datum} [datum=WGS84] - Datum to convert grid reference into. - * @returns {LatLon} Latitude/longitude of supplied grid reference. - * - * @example - * const gridref = new OsGridRef(651409.903, 313177.270); - * const pWgs84 = gridref.toLatLon(); // 52°39′28.723″N, 001°42′57.787″E - * // to obtain (historical) OSGB36 lat/lon point: - * const pOsgb = gridref.toLatLon(LatLon.datums.OSGB36); // 52°39′27.253″N, 001°43′04.518″E - */ - toLatLon(datum=LatLonEllipsoidal.datums.WGS84) { - const { easting: E, northing: N } = this; - - const { a, b } = nationalGrid.ellipsoid; // a = 6377563.396, b = 6356256.909 - const φ0 = nationalGrid.trueOrigin.lat.toRadians(); // latitude of true origin, 49°N - const λ0 = nationalGrid.trueOrigin.lon.toRadians(); // longitude of true origin, 2°W - const E0 = -nationalGrid.falseOrigin.easting; // easting of true origin, 400km - const N0 = -nationalGrid.falseOrigin.northing; // northing of true origin, -100km - const F0 = nationalGrid.scaleFactor; // 0.9996012717 - - const e2 = 1 - (b*b)/(a*a); // eccentricity squared - const n = (a-b)/(a+b), n2 = n*n, n3 = n*n*n; // n, n², n³ - - let φ=φ0, M=0; - do { - φ = (N-N0-M)/(a*F0) + φ; - - const Ma = (1 + n + (5/4)*n2 + (5/4)*n3) * (φ-φ0); - const Mb = (3*n + 3*n*n + (21/8)*n3) * Math.sin(φ-φ0) * Math.cos(φ+φ0); - const Mc = ((15/8)*n2 + (15/8)*n3) * Math.sin(2*(φ-φ0)) * Math.cos(2*(φ+φ0)); - const Md = (35/24)*n3 * Math.sin(3*(φ-φ0)) * Math.cos(3*(φ+φ0)); - M = b * F0 * (Ma - Mb + Mc - Md); // meridional arc - - } while (Math.abs(N-N0-M) >= 0.00001); // ie until < 0.01mm - - const cosφ = Math.cos(φ), sinφ = Math.sin(φ); - const ν = a*F0/Math.sqrt(1-e2*sinφ*sinφ); // nu = transverse radius of curvature - const ρ = a*F0*(1-e2)/Math.pow(1-e2*sinφ*sinφ, 1.5); // rho = meridional radius of curvature - const η2 = ν/ρ-1; // eta = ? - - const tanφ = Math.tan(φ); - const tan2φ = tanφ*tanφ, tan4φ = tan2φ*tan2φ, tan6φ = tan4φ*tan2φ; - const secφ = 1/cosφ; - const ν3 = ν*ν*ν, ν5 = ν3*ν*ν, ν7 = ν5*ν*ν; - const VII = tanφ/(2*ρ*ν); - const VIII = tanφ/(24*ρ*ν3)*(5+3*tan2φ+η2-9*tan2φ*η2); - const IX = tanφ/(720*ρ*ν5)*(61+90*tan2φ+45*tan4φ); - const X = secφ/ν; - const XI = secφ/(6*ν3)*(ν/ρ+2*tan2φ); - const XII = secφ/(120*ν5)*(5+28*tan2φ+24*tan4φ); - const XIIA = secφ/(5040*ν7)*(61+662*tan2φ+1320*tan4φ+720*tan6φ); - - const dE = (E-E0), dE2 = dE*dE, dE3 = dE2*dE, dE4 = dE2*dE2, dE5 = dE3*dE2, dE6 = dE4*dE2, dE7 = dE5*dE2; - φ = φ - VII*dE2 + VIII*dE4 - IX*dE6; - const λ = λ0 + X*dE - XI*dE3 + XII*dE5 - XIIA*dE7; - - let point = new LatLon_OsGridRef(φ.toDegrees(), λ.toDegrees(), 0, LatLonEllipsoidal.datums.OSGB36); - - if (datum != LatLonEllipsoidal.datums.OSGB36) { - // if point is required in datum other than OSGB36, convert it - point = point.convertDatum(datum); - // convertDatum() gives us a LatLon: convert to LatLon_OsGridRef which includes toOsGrid() - point = new LatLon_OsGridRef(point.lat, point.lon, point.height, point.datum); - } - - return point; - } - - - /** - * Parses grid reference to OsGridRef object. - * - * Accepts standard grid references (eg 'SU 387 148'), with or without whitespace separators, from - * two-digit references up to 10-digit references (1m × 1m square), or fully numeric comma-separated - * references in metres (eg '438700,114800'). - * - * @param {string} gridref - Standard format OS Grid Reference. - * @returns {OsGridRef} Numeric version of grid reference in metres from false origin (SW corner of - * supplied grid square). - * @throws {Error} Invalid grid reference. - * - * @example - * const grid = OsGridRef.parse('TG 51409 13177'); // grid: { easting: 651409, northing: 313177 } - */ - static parse(gridref) { - gridref = String(gridref).trim(); - - // check for fully numeric comma-separated gridref format - let match = gridref.match(/^(\d+),\s*(\d+)$/); - if (match) return new OsGridRef(match[1], match[2]); - - // validate format - match = gridref.match(/^[HNST][ABCDEFGHJKLMNOPQRSTUVWXYZ]\s*[0-9]+\s*[0-9]+$/i); - if (!match) throw new Error(`invalid grid reference ‘${gridref}’`); - - // get numeric values of letter references, mapping A->0, B->1, C->2, etc: - let l1 = gridref.toUpperCase().charCodeAt(0) - 'A'.charCodeAt(0); // 500km square - let l2 = gridref.toUpperCase().charCodeAt(1) - 'A'.charCodeAt(0); // 100km square - // shuffle down letters after 'I' since 'I' is not used in grid: - if (l1 > 7) l1--; - if (l2 > 7) l2--; - - // convert grid letters into 100km-square indexes from false origin (grid square SV): - const e100km = ((l1 - 2) % 5) * 5 + (l2 % 5); - const n100km = (19 - Math.floor(l1 / 5) * 5) - Math.floor(l2 / 5); - - // skip grid letters to get numeric (easting/northing) part of ref - let en = gridref.slice(2).trim().split(/\s+/); - // if e/n not whitespace separated, split half way - if (en.length == 1) en = [ en[0].slice(0, en[0].length / 2), en[0].slice(en[0].length / 2) ]; - - // validation - if (en[0].length != en[1].length) throw new Error(`invalid grid reference ‘${gridref}’`); - - // standardise to 10-digit refs (metres) - en[0] = en[0].padEnd(5, '0'); - en[1] = en[1].padEnd(5, '0'); - - const e = e100km + en[0]; - const n = n100km + en[1]; - - return new OsGridRef(e, n); - } - - - /** - * Converts ‘this’ numeric grid reference to standard OS Grid Reference. - * - * @param {number} [digits=10] - Precision of returned grid reference (10 digits = metres); - * digits=0 will return grid reference in numeric format. - * @returns {string} This grid reference in standard format. - * - * @example - * const gridref = new OsGridRef(651409, 313177).toString(8); // 'TG 5140 1317' - * const gridref = new OsGridRef(651409, 313177).toString(0); // '651409,313177' - */ - toString(digits=10) { - if (![ 0,2,4,6,8,10,12,14,16 ].includes(Number(digits))) throw new RangeError(`invalid precision ‘${digits}’`); // eslint-disable-line comma-spacing - - let { easting: e, northing: n } = this; - - // use digits = 0 to return numeric format (in metres) - note northing may be >= 1e7 - if (digits == 0) { - const format = { useGrouping: false, minimumIntegerDigits: 6, maximumFractionDigits: 3 }; - const ePad = e.toLocaleString('en', format); - const nPad = n.toLocaleString('en', format); - return `${ePad},${nPad}`; - } - - // get the 100km-grid indices - const e100km = Math.floor(e / 100000), n100km = Math.floor(n / 100000); - - // translate those into numeric equivalents of the grid letters - let l1 = (19 - n100km) - (19 - n100km) % 5 + Math.floor((e100km + 10) / 5); - let l2 = (19 - n100km) * 5 % 25 + e100km % 5; - - // compensate for skipped 'I' and build grid letter-pairs - if (l1 > 7) l1++; - if (l2 > 7) l2++; - const letterPair = String.fromCharCode(l1 + 'A'.charCodeAt(0), l2 + 'A'.charCodeAt(0)); - - // strip 100km-grid indices from easting & northing, and reduce precision - e = Math.floor((e % 100000) / Math.pow(10, 5 - digits / 2)); - n = Math.floor((n % 100000) / Math.pow(10, 5 - digits / 2)); - - // pad eastings & northings with leading zeros - e = e.toString().padStart(digits/2, '0'); - n = n.toString().padStart(digits/2, '0'); - - return `${letterPair} ${e} ${n}`; - } - -} - - -/* LatLon_OsGridRef - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - - -/** - * Extends LatLon class with method to convert LatLon point to OS Grid Reference. - * - * @extends LatLonEllipsoidal - */ -class LatLon_OsGridRef extends LatLonEllipsoidal { - - /** - * Converts latitude/longitude to Ordnance Survey grid reference easting/northing coordinate. - * - * @returns {OsGridRef} OS Grid Reference easting/northing. - * - * @example - * const grid = new LatLon(52.65798, 1.71605).toOsGrid(); // TG 51409 13177 - * // for conversion of (historical) OSGB36 latitude/longitude point: - * const grid = new LatLon(52.65798, 1.71605).toOsGrid(LatLon.datums.OSGB36); - */ - toOsGrid() { - // if necessary convert to OSGB36 first - const point = this.datum == LatLonEllipsoidal.datums.OSGB36 - ? this - : this.convertDatum(LatLonEllipsoidal.datums.OSGB36); - - const φ = point.lat.toRadians(); - const λ = point.lon.toRadians(); - - const { a, b } = nationalGrid.ellipsoid; // a = 6377563.396, b = 6356256.909 - const φ0 = nationalGrid.trueOrigin.lat.toRadians(); // latitude of true origin, 49°N - const λ0 = nationalGrid.trueOrigin.lon.toRadians(); // longitude of true origin, 2°W - const E0 = -nationalGrid.falseOrigin.easting; // easting of true origin, 400km - const N0 = -nationalGrid.falseOrigin.northing; // northing of true origin, -100km - const F0 = nationalGrid.scaleFactor; // 0.9996012717 - - const e2 = 1 - (b*b)/(a*a); // eccentricity squared - const n = (a-b)/(a+b), n2 = n*n, n3 = n*n*n; // n, n², n³ - - const cosφ = Math.cos(φ), sinφ = Math.sin(φ); - const ν = a*F0/Math.sqrt(1-e2*sinφ*sinφ); // nu = transverse radius of curvature - const ρ = a*F0*(1-e2)/Math.pow(1-e2*sinφ*sinφ, 1.5); // rho = meridional radius of curvature - const η2 = ν/ρ-1; // eta = ? - - const Ma = (1 + n + (5/4)*n2 + (5/4)*n3) * (φ-φ0); - const Mb = (3*n + 3*n*n + (21/8)*n3) * Math.sin(φ-φ0) * Math.cos(φ+φ0); - const Mc = ((15/8)*n2 + (15/8)*n3) * Math.sin(2*(φ-φ0)) * Math.cos(2*(φ+φ0)); - const Md = (35/24)*n3 * Math.sin(3*(φ-φ0)) * Math.cos(3*(φ+φ0)); - const M = b * F0 * (Ma - Mb + Mc - Md); // meridional arc - - const cos3φ = cosφ*cosφ*cosφ; - const cos5φ = cos3φ*cosφ*cosφ; - const tan2φ = Math.tan(φ)*Math.tan(φ); - const tan4φ = tan2φ*tan2φ; - - const I = M + N0; - const II = (ν/2)*sinφ*cosφ; - const III = (ν/24)*sinφ*cos3φ*(5-tan2φ+9*η2); - const IIIA = (ν/720)*sinφ*cos5φ*(61-58*tan2φ+tan4φ); - const IV = ν*cosφ; - const V = (ν/6)*cos3φ*(ν/ρ-tan2φ); - const VI = (ν/120) * cos5φ * (5 - 18*tan2φ + tan4φ + 14*η2 - 58*tan2φ*η2); - - const Δλ = λ-λ0; - const Δλ2 = Δλ*Δλ, Δλ3 = Δλ2*Δλ, Δλ4 = Δλ3*Δλ, Δλ5 = Δλ4*Δλ, Δλ6 = Δλ5*Δλ; - - let N = I + II*Δλ2 + III*Δλ4 + IIIA*Δλ6; - let E = E0 + IV*Δλ + V*Δλ3 + VI*Δλ5; - - N = Number(N.toFixed(3)); // round to mm precision - E = Number(E.toFixed(3)); - - try { - return new OsGridRef(E, N); // note: gets truncated to SW corner of 1m grid square - } catch (e) { - throw new Error(`${e.message} from (${point.lat.toFixed(6)},${point.lon.toFixed(6)}).toOsGrid()`); - } - } - - - /** - * Override LatLonEllipsoidal.convertDatum() with version which returns LatLon_OsGridRef. - */ - convertDatum(toDatum) { - const osgbED = super.convertDatum(toDatum); // returns LatLonEllipsoidal_Datum - const osgbOSGR = new LatLon_OsGridRef(osgbED.lat, osgbED.lon, osgbED.height, osgbED.datum); - return osgbOSGR; - } - -} - - -/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - -export { OsGridRef as default, LatLon_OsGridRef as LatLon, Dms }; diff --git a/src/js/geo/utm.js b/src/js/geo/utm.js deleted file mode 100644 index 6a6c625..0000000 --- a/src/js/geo/utm.js +++ /dev/null @@ -1,379 +0,0 @@ -/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ -/* UTM / WGS-84 Conversion Functions (c) Chris Veness 2014-2022 */ -/* MIT Licence */ -/* www.movable-type.co.uk/scripts/latlong-utm-mgrs.html */ -/* www.movable-type.co.uk/scripts/geodesy-library.html#utm */ -/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - -/* eslint-disable indent */ - -import LatLonEllipsoidal, { Dms } from './latlon-ellipsoidal-datum.js'; - - -/** - * The Universal Transverse Mercator (UTM) system is a 2-dimensional Cartesian coordinate system - * providing locations on the surface of the Earth. - * - * UTM is a set of 60 transverse Mercator projections, normally based on the WGS-84 ellipsoid. - * Within each zone, coordinates are represented as eastings and northings, measures in metres; e.g. - * ‘31 N 448251 5411932’. - * - * This method based on Karney 2011 ‘Transverse Mercator with an accuracy of a few nanometers’, - * building on Krüger 1912 ‘Konforme Abbildung des Erdellipsoids in der Ebene’. - * - * @module utm - */ - - -/* Utm - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - - -/** - * UTM coordinates, with functions to parse them and convert them to LatLon points. - */ -class Utm { - - /** - * Creates a Utm coordinate object comprising zone, hemisphere, easting, northing on a given - * datum (normally WGS84). - * - * @param {number} zone - UTM 6° longitudinal zone (1..60 covering 180°W..180°E). - * @param {string} hemisphere - N for northern hemisphere, S for southern hemisphere. - * @param {number} easting - Easting in metres from false easting (-500km from central meridian). - * @param {number} northing - Northing in metres from equator (N) or from false northing -10,000km (S). - * @param {LatLon.datums} [datum=WGS84] - Datum UTM coordinate is based on. - * @param {number} [convergence=null] - Meridian convergence (bearing of grid north - * clockwise from true north), in degrees. - * @param {number} [scale=null] - Grid scale factor. - * @params {boolean=true} verifyEN - Check easting/northing is within 'normal' values (may be - * suppressed for extended coherent coordinates or alternative datums - * e.g. ED50 (epsg.io/23029). - * @throws {TypeError} Invalid UTM coordinate. - * - * @example - * import Utm from '/js/geodesy/utm.js'; - * const utmCoord = new Utm(31, 'N', 448251, 5411932); - */ - constructor(zone, hemisphere, easting, northing, datum=LatLonEllipsoidal.datums.WGS84, convergence=null, scale=null, verifyEN=true) { - if (!(1<=zone && zone<=60)) throw new RangeError(`invalid UTM zone ‘${zone}’`); - if (zone != parseInt(zone)) throw new RangeError(`invalid UTM zone ‘${zone}’`); - if (typeof hemisphere != 'string' || !hemisphere.match(/[NS]/i)) throw new RangeError(`invalid UTM hemisphere ‘${hemisphere}’`); - if (verifyEN) { // (rough) range-check of E/N values - if (!(0<=easting && easting<=1000e3)) throw new RangeError(`invalid UTM easting ‘${easting}’`); - if (hemisphere.toUpperCase()=='N' && !(0<=northing && northing<9329006)) throw new RangeError(`invalid UTM northing ‘${northing}’`); - if (hemisphere.toUpperCase()=='S' && !(1116914 1e-12); // using IEEE 754 δτi -> 0 after 2-3 iterations - // note relatively large convergence test as δτi toggles on ±1.12e-16 for eg 31 N 400000 5000000 - const τ = τi; - - const φ = Math.atan(τ); - - let λ = Math.atan2(sinhηʹ, cosξʹ); - - // ---- convergence: Karney 2011 Eq 26, 27 - - let p = 1; - for (let j=1; j<=6; j++) p -= 2*j*β[j] * Math.cos(2*j*ξ) * Math.cosh(2*j*η); - let q = 0; - for (let j=1; j<=6; j++) q += 2*j*β[j] * Math.sin(2*j*ξ) * Math.sinh(2*j*η); - - const γʹ = Math.atan(Math.tan(ξʹ) * Math.tanh(ηʹ)); - const γʺ = Math.atan2(q, p); - - const γ = γʹ + γʺ; - - // ---- scale: Karney 2011 Eq 28 - - const sinφ = Math.sin(φ); - const kʹ = Math.sqrt(1 - e*e*sinφ*sinφ) * Math.sqrt(1 + τ*τ) * Math.sqrt(sinhηʹ*sinhηʹ + cosξʹ*cosξʹ); - const kʺ = A / a / Math.sqrt(p*p + q*q); - - const k = k0 * kʹ * kʺ; - - // ------------ - - const λ0 = ((z-1)*6 - 180 + 3).toRadians(); // longitude of central meridian - λ += λ0; // move λ from zonal to global coordinates - - // round to reasonable precision - const lat = Number(φ.toDegrees().toFixed(14)); // nm precision (1nm = 10^-14°) - const lon = Number(λ.toDegrees().toFixed(14)); // (strictly lat rounding should be φ⋅cosφ!) - const convergence = Number(γ.toDegrees().toFixed(9)); - const scale = Number(k.toFixed(12)); - - const latLong = new LatLon_Utm(lat, lon, 0, this.datum); - // ... and add the convergence and scale into the LatLon object ... wonderful JavaScript! - latLong.convergence = convergence; - latLong.scale = scale; - - return latLong; - } - - - /** - * Parses string representation of UTM coordinate. - * - * A UTM coordinate comprises (space-separated) - * - zone - * - hemisphere - * - easting - * - northing. - * - * @param {string} utmCoord - UTM coordinate (WGS 84). - * @param {Datum} [datum=WGS84] - Datum coordinate is defined in (default WGS 84). - * @returns {Utm} Parsed UTM coordinate. - * @throws {TypeError} Invalid UTM coordinate. - * - * @example - * const utmCoord = Utm.parse('31 N 448251 5411932'); - * // utmCoord: {zone: 31, hemisphere: 'N', easting: 448251, northing: 5411932 } - */ - static parse(utmCoord, datum=LatLonEllipsoidal.datums.WGS84) { - // match separate elements (separated by whitespace) - utmCoord = utmCoord.trim().match(/\S+/g); - - if (utmCoord==null || utmCoord.length!=4) throw new Error(`invalid UTM coordinate ‘${utmCoord}’`); - - const zone = utmCoord[0], hemisphere = utmCoord[1], easting = utmCoord[2], northing = utmCoord[3]; - - return new this(zone, hemisphere, easting, northing, datum); // 'new this' as may return subclassed types - } - - - /** - * Returns a string representation of a UTM coordinate. - * - * To distinguish from MGRS grid zone designators, a space is left between the zone and the - * hemisphere. - * - * Note that UTM coordinates get rounded, not truncated (unlike MGRS grid references). - * - * @param {number} [digits=0] - Number of digits to appear after the decimal point (3 ≡ mm). - * @returns {string} A string representation of the coordinate. - * - * @example - * const utm = new Utm('31', 'N', 448251, 5411932).toString(4); // 31 N 448251.0000 5411932.0000 - */ - toString(digits=0) { - - const z = this.zone.toString().padStart(2, '0'); - const h = this.hemisphere; - const e = this.easting.toFixed(digits); - const n = this.northing.toFixed(digits); - - return `${z} ${h} ${e} ${n}`; - } - -} - - -/* LatLon_Utm - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - - -/** - * Extends LatLon with method to convert LatLon points to UTM coordinates. - * - * @extends LatLon - */ -class LatLon_Utm extends LatLonEllipsoidal { - - /** - * Converts latitude/longitude to UTM coordinate. - * - * Implements Karney’s method, using Krüger series to order n⁶, giving results accurate to 5nm - * for distances up to 3900km from the central meridian. - * - * @param {number} [zoneOverride] - Use specified zone rather than zone within which point lies; - * note overriding the UTM zone has the potential to result in negative eastings, and - * perverse results within Norway/Svalbard exceptions. - * @returns {Utm} UTM coordinate. - * @throws {TypeError} Latitude outside UTM limits. - * - * @example - * const latlong = new LatLon(48.8582, 2.2945); - * const utmCoord = latlong.toUtm(); // 31 N 448252 5411933 - */ - toUtm(zoneOverride=undefined) { - if (!(-80<=this.lat && this.lat<=84)) throw new RangeError(`latitude ‘${this.lat}’ outside UTM limits`); - - const falseEasting = 500e3, falseNorthing = 10000e3; - - let zone = zoneOverride || Math.floor((this.lon+180)/6) + 1; // longitudinal zone - let λ0 = ((zone-1)*6 - 180 + 3).toRadians(); // longitude of central meridian - - // ---- handle Norway/Svalbard exceptions - // grid zones are 8° tall; 0°N is offset 10 into latitude bands array - const mgrsLatBands = 'CDEFGHJKLMNPQRSTUVWXX'; // X is repeated for 80-84°N - const latBand = mgrsLatBands.charAt(Math.floor(this.lat/8+10)); - // adjust zone & central meridian for Norway - if (zone==31 && latBand=='V' && this.lon>= 3) { zone++; λ0 += (6).toRadians(); } - // adjust zone & central meridian for Svalbard - if (zone==32 && latBand=='X' && this.lon< 9) { zone--; λ0 -= (6).toRadians(); } - if (zone==32 && latBand=='X' && this.lon>= 9) { zone++; λ0 += (6).toRadians(); } - if (zone==34 && latBand=='X' && this.lon< 21) { zone--; λ0 -= (6).toRadians(); } - if (zone==34 && latBand=='X' && this.lon>=21) { zone++; λ0 += (6).toRadians(); } - if (zone==36 && latBand=='X' && this.lon< 33) { zone--; λ0 -= (6).toRadians(); } - if (zone==36 && latBand=='X' && this.lon>=33) { zone++; λ0 += (6).toRadians(); } - - const φ = this.lat.toRadians(); // latitude ± from equator - const λ = this.lon.toRadians() - λ0; // longitude ± from central meridian - - // allow alternative ellipsoid to be specified - const ellipsoid = this.datum ? this.datum.ellipsoid : LatLonEllipsoidal.ellipsoids.WGS84; - const { a, f } = ellipsoid; // WGS-84: a = 6378137, f = 1/298.257223563; - - const k0 = 0.9996; // UTM scale on the central meridian - - // ---- easting, northing: Karney 2011 Eq 7-14, 29, 35: - - const e = Math.sqrt(f*(2-f)); // eccentricity - const n = f / (2 - f); // 3rd flattening - const n2 = n*n, n3 = n*n2, n4 = n*n3, n5 = n*n4, n6 = n*n5; - - const cosλ = Math.cos(λ), sinλ = Math.sin(λ), tanλ = Math.tan(λ); - - const τ = Math.tan(φ); // τ ≡ tanφ, τʹ ≡ tanφʹ; prime (ʹ) indicates angles on the conformal sphere - const σ = Math.sinh(e*Math.atanh(e*τ/Math.sqrt(1+τ*τ))); - - const τʹ = τ*Math.sqrt(1+σ*σ) - σ*Math.sqrt(1+τ*τ); - - const ξʹ = Math.atan2(τʹ, cosλ); - const ηʹ = Math.asinh(sinλ / Math.sqrt(τʹ*τʹ + cosλ*cosλ)); - - const A = a/(1+n) * (1 + 1/4*n2 + 1/64*n4 + 1/256*n6); // 2πA is the circumference of a meridian - - const α = [ null, // note α is one-based array (6th order Krüger expressions) - 1/2*n - 2/3*n2 + 5/16*n3 + 41/180*n4 - 127/288*n5 + 7891/37800*n6, - 13/48*n2 - 3/5*n3 + 557/1440*n4 + 281/630*n5 - 1983433/1935360*n6, - 61/240*n3 - 103/140*n4 + 15061/26880*n5 + 167603/181440*n6, - 49561/161280*n4 - 179/168*n5 + 6601661/7257600*n6, - 34729/80640*n5 - 3418889/1995840*n6, - 212378941/319334400*n6 ]; - - let ξ = ξʹ; - for (let j=1; j<=6; j++) ξ += α[j] * Math.sin(2*j*ξʹ) * Math.cosh(2*j*ηʹ); - - let η = ηʹ; - for (let j=1; j<=6; j++) η += α[j] * Math.cos(2*j*ξʹ) * Math.sinh(2*j*ηʹ); - - let x = k0 * A * η; - let y = k0 * A * ξ; - - // ---- convergence: Karney 2011 Eq 23, 24 - - let pʹ = 1; - for (let j=1; j<=6; j++) pʹ += 2*j*α[j] * Math.cos(2*j*ξʹ) * Math.cosh(2*j*ηʹ); - let qʹ = 0; - for (let j=1; j<=6; j++) qʹ += 2*j*α[j] * Math.sin(2*j*ξʹ) * Math.sinh(2*j*ηʹ); - - const γʹ = Math.atan(τʹ / Math.sqrt(1+τʹ*τʹ)*tanλ); - const γʺ = Math.atan2(qʹ, pʹ); - - const γ = γʹ + γʺ; - - // ---- scale: Karney 2011 Eq 25 - - const sinφ = Math.sin(φ); - const kʹ = Math.sqrt(1 - e*e*sinφ*sinφ) * Math.sqrt(1 + τ*τ) / Math.sqrt(τʹ*τʹ + cosλ*cosλ); - const kʺ = A / a * Math.sqrt(pʹ*pʹ + qʹ*qʹ); - - const k = k0 * kʹ * kʺ; - - // ------------ - - // shift x/y to false origins - x = x + falseEasting; // make x relative to false easting - if (y < 0) y = y + falseNorthing; // make y in southern hemisphere relative to false northing - - // round to reasonable precision - x = Number(x.toFixed(9)); // nm precision - y = Number(y.toFixed(9)); // nm precision - const convergence = Number(γ.toDegrees().toFixed(9)); - const scale = Number(k.toFixed(12)); - - const h = this.lat>=0 ? 'N' : 'S'; // hemisphere - - return new Utm(zone, h, x, y, this.datum, convergence, scale, !!zoneOverride); - } -} - - -/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - -export { Utm as default, LatLon_Utm as LatLon, Dms }; diff --git a/src/js/geo/vector3d.js b/src/js/geo/vector3d.js deleted file mode 100644 index 08ddce6..0000000 --- a/src/js/geo/vector3d.js +++ /dev/null @@ -1,256 +0,0 @@ -/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ -/* Vector handling functions (c) Chris Veness 2011-2019 */ -/* MIT Licence */ -/* www.movable-type.co.uk/scripts/geodesy-library.html#vector3d */ -/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - - -/** - * Library of 3-d vector manipulation routines. - * - * @module vector3d - */ - - -/* Vector3d - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ - - -/** - * Functions for manipulating generic 3-d vectors. - * - * Functions return vectors as return results, so that operations can be chained. - * - * @example - * const v = v1.cross(v2).dot(v3) // ≡ v1×v2⋅v3 - */ -class Vector3d { - - /** - * Creates a 3-d vector. - * - * @param {number} x - X component of vector. - * @param {number} y - Y component of vector. - * @param {number} z - Z component of vector. - * - * @example - * import Vector3d from '/js/geodesy/vector3d.js'; - * const v = new Vector3d(0.267, 0.535, 0.802); - */ - constructor(x, y, z) { - if (isNaN(x) || isNaN(y) || isNaN(z)) throw new TypeError(`invalid vector [${x},${y},${z}]`); - - this.x = Number(x); - this.y = Number(y); - this.z = Number(z); - } - - - /** - * Length (magnitude or norm) of ‘this’ vector. - * - * @returns {number} Magnitude of this vector. - */ - get length() { - return Math.sqrt(this.x * this.x + this.y * this.y + this.z * this.z); - } - - - /** - * Adds supplied vector to ‘this’ vector. - * - * @param {Vector3d} v - Vector to be added to this vector. - * @returns {Vector3d} Vector representing sum of this and v. - */ - plus(v) { - if (!(v instanceof Vector3d)) throw new TypeError('v is not Vector3d object'); - - return new Vector3d(this.x + v.x, this.y + v.y, this.z + v.z); - } - - - /** - * Subtracts supplied vector from ‘this’ vector. - * - * @param {Vector3d} v - Vector to be subtracted from this vector. - * @returns {Vector3d} Vector representing difference between this and v. - */ - minus(v) { - if (!(v instanceof Vector3d)) throw new TypeError('v is not Vector3d object'); - - return new Vector3d(this.x - v.x, this.y - v.y, this.z - v.z); - } - - - /** - * Multiplies ‘this’ vector by a scalar value. - * - * @param {number} x - Factor to multiply this vector by. - * @returns {Vector3d} Vector scaled by x. - */ - times(x) { - if (isNaN(x)) throw new TypeError(`invalid scalar value ‘${x}’`); - - return new Vector3d(this.x * x, this.y * x, this.z * x); - } - - - /** - * Divides ‘this’ vector by a scalar value. - * - * @param {number} x - Factor to divide this vector by. - * @returns {Vector3d} Vector divided by x. - */ - dividedBy(x) { - if (isNaN(x)) throw new TypeError(`invalid scalar value ‘${x}’`); - - return new Vector3d(this.x / x, this.y / x, this.z / x); - } - - - /** - * Multiplies ‘this’ vector by the supplied vector using dot (scalar) product. - * - * @param {Vector3d} v - Vector to be dotted with this vector. - * @returns {number} Dot product of ‘this’ and v. - */ - dot(v) { - if (!(v instanceof Vector3d)) throw new TypeError('v is not Vector3d object'); - - return this.x * v.x + this.y * v.y + this.z * v.z; - } - - - /** - * Multiplies ‘this’ vector by the supplied vector using cross (vector) product. - * - * @param {Vector3d} v - Vector to be crossed with this vector. - * @returns {Vector3d} Cross product of ‘this’ and v. - */ - cross(v) { - if (!(v instanceof Vector3d)) throw new TypeError('v is not Vector3d object'); - - const x = this.y * v.z - this.z * v.y; - const y = this.z * v.x - this.x * v.z; - const z = this.x * v.y - this.y * v.x; - - return new Vector3d(x, y, z); - } - - - /** - * Negates a vector to point in the opposite direction. - * - * @returns {Vector3d} Negated vector. - */ - negate() { - return new Vector3d(-this.x, -this.y, -this.z); - } - - - /** - * Normalizes a vector to its unit vector - * – if the vector is already unit or is zero magnitude, this is a no-op. - * - * @returns {Vector3d} Normalised version of this vector. - */ - unit() { - const norm = this.length; - if (norm == 1) return this; - if (norm == 0) return this; - - const x = this.x / norm; - const y = this.y / norm; - const z = this.z / norm; - - return new Vector3d(x, y, z); - } - - - /** - * Calculates the angle between ‘this’ vector and supplied vector atan2(|p₁×p₂|, p₁·p₂) (or if - * (extra-planar) ‘n’ supplied then atan2(n·p₁×p₂, p₁·p₂). - * - * @param {Vector3d} v - Vector whose angle is to be determined from ‘this’ vector. - * @param {Vector3d} [n] - Plane normal: if supplied, angle is signed +ve if this->v is - * clockwise looking along n, -ve in opposite direction. - * @returns {number} Angle (in radians) between this vector and supplied vector (in range 0..π - * if n not supplied, range -π..+π if n supplied). - */ - angleTo(v, n=undefined) { - if (!(v instanceof Vector3d)) throw new TypeError('v is not Vector3d object'); - if (!(n instanceof Vector3d || n == undefined)) throw new TypeError('n is not Vector3d object'); - - // q.v. stackoverflow.com/questions/14066933#answer-16544330, but n·p₁×p₂ is numerically - // ill-conditioned, so just calculate sign to apply to |p₁×p₂| - - // if n·p₁×p₂ is -ve, negate |p₁×p₂| - const sign = n==undefined || this.cross(v).dot(n)>=0 ? 1 : -1; - - const sinθ = this.cross(v).length * sign; - const cosθ = this.dot(v); - - return Math.atan2(sinθ, cosθ); - } - - - /** - * Rotates ‘this’ point around an axis by a specified angle. - * - * @param {Vector3d} axis - The axis being rotated around. - * @param {number} angle - The angle of rotation (in degrees). - * @returns {Vector3d} The rotated point. - */ - rotateAround(axis, angle) { - if (!(axis instanceof Vector3d)) throw new TypeError('axis is not Vector3d object'); - - const θ = angle.toRadians(); - - // en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle - // en.wikipedia.org/wiki/Quaternions_and_spatial_rotation#Quaternion-derived_rotation_matrix - const p = this.unit(); - const a = axis.unit(); - - const s = Math.sin(θ); - const c = Math.cos(θ); - const t = 1-c; - const x = a.x, y = a.y, z = a.z; - - const r = [ // rotation matrix for rotation about supplied axis - [ t*x*x + c, t*x*y - s*z, t*x*z + s*y ], - [ t*x*y + s*z, t*y*y + c, t*y*z - s*x ], - [ t*x*z - s*y, t*y*z + s*x, t*z*z + c ], - ]; - - // multiply r × p - const rp = [ - r[0][0]*p.x + r[0][1]*p.y + r[0][2]*p.z, - r[1][0]*p.x + r[1][1]*p.y + r[1][2]*p.z, - r[2][0]*p.x + r[2][1]*p.y + r[2][2]*p.z, - ]; - const p2 = new Vector3d(rp[0], rp[1], rp[2]); - - return p2; - // qv en.wikipedia.org/wiki/Rodrigues'_rotation_formula... - } - - - /** - * String representation of vector. - * - * @param {number} [dp=3] - Number of decimal places to be used. - * @returns {string} Vector represented as [x,y,z]. - */ - toString(dp=3) { - return `[${this.x.toFixed(dp)},${this.y.toFixed(dp)},${this.z.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 default Vector3d; diff --git a/src/js/globalize.min.js b/src/js/globalize.min.js deleted file mode 100644 index d60d8db..0000000 --- a/src/js/globalize.min.js +++ /dev/null @@ -1,5 +0,0 @@ -/*! - * Globalize v1.7.0 2021-08-02T11:53Z Released under the MIT license - * http://git.io/TrdQbw - */ -!function(e,n){"function"==typeof define&&define.amd?define(["cldr","cldr/event"],n):"object"==typeof exports?module.exports=n(require("cldrjs")):e.Globalize=n(e.Cldr)}(this,(function(e){var n=function(e,n){return e=e.replace(/{[0-9a-zA-Z-_. ]+}/g,(function(e){return e=e.replace(/^{([^}]*)}$/,"$1"),"string"==typeof(t=n[e])?t:"number"==typeof t?""+t:JSON.stringify(t);var t}))},t=function(){var e=arguments[0],n=[].slice.call(arguments,1);return n.forEach((function(n){var t;for(t in n)e[t]=n[t]})),e},r=function(e,r,a){var i;return r=e+(r?": "+n(r,a):""),(i=new Error(r)).code=e,t(i,a),i},a=function(e,n,t){e.length&&e[e.length-1].type===n?e[e.length-1].value+=t:e.push({type:n,value:t})},i=function(e){return JSON.stringify(e,(function(e,n){return n&&n.runtimeKey?n.runtimeKey:n}))},u=function(e,n,t,a){if(!t)throw r(e,n,a)},o=function(e){return Array.isArray(e)?e:e?[e]:[]},c=function(e,n,t){var r;r=o((t=t||{}).skip).some((function(n){return n.test(e)})),u("E_MISSING_CLDR","Missing required CLDR content `{path}`.",n||r,{path:e})},l=function(e,n){u("E_MISSING_PARAMETER","Missing required parameter `{name}`.",void 0!==e,{name:n})},f=function(e,n,t,r){u("E_INVALID_PAR_TYPE","Invalid `{name}` parameter ({value}). {expected} expected.",t,{expected:r,name:n,value:e})},s=function(n,t){f(n,t,void 0===n||"string"==typeof n||n instanceof e,"String or Cldr instance")},d=function(e){return null!==e&&""+e=="[object Object]"},m=function(n){return n instanceof e?n:new e(n)};function v(e){e.once("get",c),e.get("supplemental/likelySubtags")}function _(e){if(!(this instanceof _))return new _(e);l(e,"locale"),s(e,"locale"),this.cldr=m(e),v(this.cldr)}return _.load=function(){e.load.apply(e,arguments)},_.locale=function(e){return s(e,"locale"),arguments.length&&(this.cldr=m(e),v(this.cldr)),this.cldr},_._alwaysArray=o,_._createError=r,_._formatMessage=n,_._formatMessageToParts=function(e,n){var t=0,r=[];return e.replace(/{[0-9a-zA-Z-_. ]+}/g,(function(i,u){var o=i.slice(1,-1);a(r,"literal",e.slice(t,u)),a(r,"variable",n[o]),r[r.length-1].name=o,t+=u+i.length})),r.filter((function(e){return""!==e.value}))},_._isPlainObject=d,_._objectExtend=t,_._partsJoin=function(e){return e.map((function(e){return e.value})).join("")},_._partsPush=a,_._regexpEscape=function(e){return e.replace(/([.*+?^=!:${}()|\[\]\/\\])/g,"\\$1")},_._runtimeBind=function(e,n,t,r){var a=i(e),u=function(e){if(void 0!==e.name)return e.name;var n=/^function\s+([\w\$]+)\s*\(/.exec(e.toString());return n&&n.length>0?n[1]:void 0}(t),o=n.locale;return u?(t.runtimeKey=function(e,n,t,r){var a,u;return r=r||i(t),u=e+n+r,(a=[].reduce.call(u,(function(e,n){return 0|(e=(e<<5)-e+n.charCodeAt(0))}),0))>0?"a"+a:"b"+Math.abs(a)}(u,o,null,a),t.generatorString=function(){return'Globalize("'+o+'").'+u+"("+a.slice(1,-1)+")"},t.runtimeArgs=r,t):t},_._stringPad=function(e,n,t){var r;for("string"!=typeof e&&(e=String(e)),r=e.length;r=t&&e<=r,{maximum:r,minimum:t,name:n,value:e})},_._validateParameterTypePlainObject=function(e,n){f(e,n,void 0===e||d(e),"Plain Object")},_._validateParameterType=f,_})); diff --git a/src/js/jquery.js b/src/js/jquery.js deleted file mode 100644 index b86de89..0000000 --- a/src/js/jquery.js +++ /dev/null @@ -1,10993 +0,0 @@ -/*! - * jQuery JavaScript Library v3.6.3 - * https://jquery.com/ - * - * Includes Sizzle.js - * https://sizzlejs.com/ - * - * Copyright OpenJS Foundation and other contributors - * Released under the MIT license - * https://jquery.org/license - * - * Date: 2022-12-20T21:28Z - */ -( function( global, factory ) { - - "use strict"; - - if ( typeof module === "object" && typeof module.exports === "object" ) { - - // For CommonJS and CommonJS-like environments where a proper `window` - // is present, execute the factory and get jQuery. - // For environments that do not have a `window` with a `document` - // (such as Node.js), expose a factory as module.exports. - // This accentuates the need for the creation of a real `window`. - // e.g. var jQuery = require("jquery")(window); - // See ticket trac-14549 for more info. - module.exports = global.document ? - factory( global, true ) : - function( w ) { - if ( !w.document ) { - throw new Error( "jQuery requires a window with a document" ); - } - return factory( w ); - }; - } else { - factory( global ); - } - -// Pass this if window is not defined yet -} )( typeof window !== "undefined" ? window : this, function( window, noGlobal ) { - -// Edge <= 12 - 13+, Firefox <=18 - 45+, IE 10 - 11, Safari 5.1 - 9+, iOS 6 - 9.1 -// throw exceptions when non-strict code (e.g., ASP.NET 4.5) accesses strict mode -// arguments.callee.caller (trac-13335). But as of jQuery 3.0 (2016), strict mode should be common -// enough that all such attempts are guarded in a try block. -"use strict"; - -var arr = []; - -var getProto = Object.getPrototypeOf; - -var slice = arr.slice; - -var flat = arr.flat ? function( array ) { - return arr.flat.call( array ); -} : function( array ) { - return arr.concat.apply( [], array ); -}; - - -var push = arr.push; - -var indexOf = arr.indexOf; - -var class2type = {}; - -var toString = class2type.toString; - -var hasOwn = class2type.hasOwnProperty; - -var fnToString = hasOwn.toString; - -var ObjectFunctionString = fnToString.call( Object ); - -var support = {}; - -var isFunction = function isFunction( obj ) { - - // Support: Chrome <=57, Firefox <=52 - // In some browsers, typeof returns "function" for HTML elements - // (i.e., `typeof document.createElement( "object" ) === "function"`). - // We don't want to classify *any* DOM node as a function. - // Support: QtWeb <=3.8.5, WebKit <=534.34, wkhtmltopdf tool <=0.12.5 - // Plus for old WebKit, typeof returns "function" for HTML collections - // (e.g., `typeof document.getElementsByTagName("div") === "function"`). (gh-4756) - return typeof obj === "function" && typeof obj.nodeType !== "number" && - typeof obj.item !== "function"; - }; - - -var isWindow = function isWindow( obj ) { - return obj != null && obj === obj.window; - }; - - -var document = window.document; - - - - var preservedScriptAttributes = { - type: true, - src: true, - nonce: true, - noModule: true - }; - - function DOMEval( code, node, doc ) { - doc = doc || document; - - var i, val, - script = doc.createElement( "script" ); - - script.text = code; - if ( node ) { - for ( i in preservedScriptAttributes ) { - - // Support: Firefox 64+, Edge 18+ - // Some browsers don't support the "nonce" property on scripts. - // On the other hand, just using `getAttribute` is not enough as - // the `nonce` attribute is reset to an empty string whenever it - // becomes browsing-context connected. - // See https://github.com/whatwg/html/issues/2369 - // See https://html.spec.whatwg.org/#nonce-attributes - // The `node.getAttribute` check was added for the sake of - // `jQuery.globalEval` so that it can fake a nonce-containing node - // via an object. - val = node[ i ] || node.getAttribute && node.getAttribute( i ); - if ( val ) { - script.setAttribute( i, val ); - } - } - } - doc.head.appendChild( script ).parentNode.removeChild( script ); - } - - -function toType( obj ) { - if ( obj == null ) { - return obj + ""; - } - - // Support: Android <=2.3 only (functionish RegExp) - return typeof obj === "object" || typeof obj === "function" ? - class2type[ toString.call( obj ) ] || "object" : - typeof obj; -} -/* global Symbol */ -// Defining this global in .eslintrc.json would create a danger of using the global -// unguarded in another place, it seems safer to define global only for this module - - - -var - version = "3.6.3", - - // Define a local copy of jQuery - jQuery = function( selector, context ) { - - // The jQuery object is actually just the init constructor 'enhanced' - // Need init if jQuery is called (just allow error to be thrown if not included) - return new jQuery.fn.init( selector, context ); - }; - -jQuery.fn = jQuery.prototype = { - - // The current version of jQuery being used - jquery: version, - - constructor: jQuery, - - // The default length of a jQuery object is 0 - length: 0, - - toArray: function() { - return slice.call( this ); - }, - - // Get the Nth element in the matched element set OR - // Get the whole matched element set as a clean array - get: function( num ) { - - // Return all the elements in a clean array - if ( num == null ) { - return slice.call( this ); - } - - // Return just the one element from the set - return num < 0 ? this[ num + this.length ] : this[ num ]; - }, - - // Take an array of elements and push it onto the stack - // (returning the new matched element set) - pushStack: function( elems ) { - - // Build a new jQuery matched element set - var ret = jQuery.merge( this.constructor(), elems ); - - // Add the old object onto the stack (as a reference) - ret.prevObject = this; - - // Return the newly-formed element set - return ret; - }, - - // Execute a callback for every element in the matched set. - each: function( callback ) { - return jQuery.each( this, callback ); - }, - - map: function( callback ) { - return this.pushStack( jQuery.map( this, function( elem, i ) { - return callback.call( elem, i, elem ); - } ) ); - }, - - slice: function() { - return this.pushStack( slice.apply( this, arguments ) ); - }, - - first: function() { - return this.eq( 0 ); - }, - - last: function() { - return this.eq( -1 ); - }, - - even: function() { - return this.pushStack( jQuery.grep( this, function( _elem, i ) { - return ( i + 1 ) % 2; - } ) ); - }, - - odd: function() { - return this.pushStack( jQuery.grep( this, function( _elem, i ) { - return i % 2; - } ) ); - }, - - eq: function( i ) { - var len = this.length, - j = +i + ( i < 0 ? len : 0 ); - return this.pushStack( j >= 0 && j < len ? [ this[ j ] ] : [] ); - }, - - end: function() { - return this.prevObject || this.constructor(); - }, - - // For internal use only. - // Behaves like an Array's method, not like a jQuery method. - push: push, - sort: arr.sort, - splice: arr.splice -}; - -jQuery.extend = jQuery.fn.extend = function() { - var options, name, src, copy, copyIsArray, clone, - target = arguments[ 0 ] || {}, - i = 1, - length = arguments.length, - deep = false; - - // Handle a deep copy situation - if ( typeof target === "boolean" ) { - deep = target; - - // Skip the boolean and the target - target = arguments[ i ] || {}; - i++; - } - - // Handle case when target is a string or something (possible in deep copy) - if ( typeof target !== "object" && !isFunction( target ) ) { - target = {}; - } - - // Extend jQuery itself if only one argument is passed - if ( i === length ) { - target = this; - i--; - } - - for ( ; i < length; i++ ) { - - // Only deal with non-null/undefined values - if ( ( options = arguments[ i ] ) != null ) { - - // Extend the base object - for ( name in options ) { - copy = options[ name ]; - - // Prevent Object.prototype pollution - // Prevent never-ending loop - if ( name === "__proto__" || target === copy ) { - continue; - } - - // Recurse if we're merging plain objects or arrays - if ( deep && copy && ( jQuery.isPlainObject( copy ) || - ( copyIsArray = Array.isArray( copy ) ) ) ) { - src = target[ name ]; - - // Ensure proper type for the source value - if ( copyIsArray && !Array.isArray( src ) ) { - clone = []; - } else if ( !copyIsArray && !jQuery.isPlainObject( src ) ) { - clone = {}; - } else { - clone = src; - } - copyIsArray = false; - - // Never move original objects, clone them - target[ name ] = jQuery.extend( deep, clone, copy ); - - // Don't bring in undefined values - } else if ( copy !== undefined ) { - target[ name ] = copy; - } - } - } - } - - // Return the modified object - return target; -}; - -jQuery.extend( { - - // Unique for each copy of jQuery on the page - expando: "jQuery" + ( version + Math.random() ).replace( /\D/g, "" ), - - // Assume jQuery is ready without the ready module - isReady: true, - - error: function( msg ) { - throw new Error( msg ); - }, - - noop: function() {}, - - isPlainObject: function( obj ) { - var proto, Ctor; - - // Detect obvious negatives - // Use toString instead of jQuery.type to catch host objects - if ( !obj || toString.call( obj ) !== "[object Object]" ) { - return false; - } - - proto = getProto( obj ); - - // Objects with no prototype (e.g., `Object.create( null )`) are plain - if ( !proto ) { - return true; - } - - // Objects with prototype are plain iff they were constructed by a global Object function - Ctor = hasOwn.call( proto, "constructor" ) && proto.constructor; - return typeof Ctor === "function" && fnToString.call( Ctor ) === ObjectFunctionString; - }, - - isEmptyObject: function( obj ) { - var name; - - for ( name in obj ) { - return false; - } - return true; - }, - - // Evaluates a script in a provided context; falls back to the global one - // if not specified. - globalEval: function( code, options, doc ) { - DOMEval( code, { nonce: options && options.nonce }, doc ); - }, - - each: function( obj, callback ) { - var length, i = 0; - - if ( isArrayLike( obj ) ) { - length = obj.length; - for ( ; i < length; i++ ) { - if ( callback.call( obj[ i ], i, obj[ i ] ) === false ) { - break; - } - } - } else { - for ( i in obj ) { - if ( callback.call( obj[ i ], i, obj[ i ] ) === false ) { - break; - } - } - } - - return obj; - }, - - // results is for internal usage only - makeArray: function( arr, results ) { - var ret = results || []; - - if ( arr != null ) { - if ( isArrayLike( Object( arr ) ) ) { - jQuery.merge( ret, - typeof arr === "string" ? - [ arr ] : arr - ); - } else { - push.call( ret, arr ); - } - } - - return ret; - }, - - inArray: function( elem, arr, i ) { - return arr == null ? -1 : indexOf.call( arr, elem, i ); - }, - - // Support: Android <=4.0 only, PhantomJS 1 only - // push.apply(_, arraylike) throws on ancient WebKit - merge: function( first, second ) { - var len = +second.length, - j = 0, - i = first.length; - - for ( ; j < len; j++ ) { - first[ i++ ] = second[ j ]; - } - - first.length = i; - - return first; - }, - - grep: function( elems, callback, invert ) { - var callbackInverse, - matches = [], - i = 0, - length = elems.length, - callbackExpect = !invert; - - // Go through the array, only saving the items - // that pass the validator function - for ( ; i < length; i++ ) { - callbackInverse = !callback( elems[ i ], i ); - if ( callbackInverse !== callbackExpect ) { - matches.push( elems[ i ] ); - } - } - - return matches; - }, - - // arg is for internal usage only - map: function( elems, callback, arg ) { - var length, value, - i = 0, - ret = []; - - // Go through the array, translating each of the items to their new values - if ( isArrayLike( elems ) ) { - length = elems.length; - for ( ; i < length; i++ ) { - value = callback( elems[ i ], i, arg ); - - if ( value != null ) { - ret.push( value ); - } - } - - // Go through every key on the object, - } else { - for ( i in elems ) { - value = callback( elems[ i ], i, arg ); - - if ( value != null ) { - ret.push( value ); - } - } - } - - // Flatten any nested arrays - return flat( ret ); - }, - - // A global GUID counter for objects - guid: 1, - - // jQuery.support is not used in Core but other projects attach their - // properties to it so it needs to exist. - support: support -} ); - -if ( typeof Symbol === "function" ) { - jQuery.fn[ Symbol.iterator ] = arr[ Symbol.iterator ]; -} - -// Populate the class2type map -jQuery.each( "Boolean Number String Function Array Date RegExp Object Error Symbol".split( " " ), - function( _i, name ) { - class2type[ "[object " + name + "]" ] = name.toLowerCase(); - } ); - -function isArrayLike( obj ) { - - // Support: real iOS 8.2 only (not reproducible in simulator) - // `in` check used to prevent JIT error (gh-2145) - // hasOwn isn't used here due to false negatives - // regarding Nodelist length in IE - var length = !!obj && "length" in obj && obj.length, - type = toType( obj ); - - if ( isFunction( obj ) || isWindow( obj ) ) { - return false; - } - - return type === "array" || length === 0 || - typeof length === "number" && length > 0 && ( length - 1 ) in obj; -} -var Sizzle = -/*! - * Sizzle CSS Selector Engine v2.3.9 - * https://sizzlejs.com/ - * - * Copyright JS Foundation and other contributors - * Released under the MIT license - * https://js.foundation/ - * - * Date: 2022-12-19 - */ -( function( window ) { -var i, - support, - Expr, - getText, - isXML, - tokenize, - compile, - select, - outermostContext, - sortInput, - hasDuplicate, - - // Local document vars - setDocument, - document, - docElem, - documentIsHTML, - rbuggyQSA, - rbuggyMatches, - matches, - contains, - - // Instance-specific data - expando = "sizzle" + 1 * new Date(), - preferredDoc = window.document, - dirruns = 0, - done = 0, - classCache = createCache(), - tokenCache = createCache(), - compilerCache = createCache(), - nonnativeSelectorCache = createCache(), - sortOrder = function( a, b ) { - if ( a === b ) { - hasDuplicate = true; - } - return 0; - }, - - // Instance methods - hasOwn = ( {} ).hasOwnProperty, - arr = [], - pop = arr.pop, - pushNative = arr.push, - push = arr.push, - slice = arr.slice, - - // Use a stripped-down indexOf as it's faster than native - // https://jsperf.com/thor-indexof-vs-for/5 - indexOf = function( list, elem ) { - var i = 0, - len = list.length; - for ( ; i < len; i++ ) { - if ( list[ i ] === elem ) { - return i; - } - } - return -1; - }, - - booleans = "checked|selected|async|autofocus|autoplay|controls|defer|disabled|hidden|" + - "ismap|loop|multiple|open|readonly|required|scoped", - - // Regular expressions - - // http://www.w3.org/TR/css3-selectors/#whitespace - whitespace = "[\\x20\\t\\r\\n\\f]", - - // https://www.w3.org/TR/css-syntax-3/#ident-token-diagram - identifier = "(?:\\\\[\\da-fA-F]{1,6}" + whitespace + - "?|\\\\[^\\r\\n\\f]|[\\w-]|[^\0-\\x7f])+", - - // Attribute selectors: http://www.w3.org/TR/selectors/#attribute-selectors - attributes = "\\[" + whitespace + "*(" + identifier + ")(?:" + whitespace + - - // Operator (capture 2) - "*([*^$|!~]?=)" + whitespace + - - // "Attribute values must be CSS identifiers [capture 5] - // or strings [capture 3 or capture 4]" - "*(?:'((?:\\\\.|[^\\\\'])*)'|\"((?:\\\\.|[^\\\\\"])*)\"|(" + identifier + "))|)" + - whitespace + "*\\]", - - pseudos = ":(" + identifier + ")(?:\\((" + - - // To reduce the number of selectors needing tokenize in the preFilter, prefer arguments: - // 1. quoted (capture 3; capture 4 or capture 5) - "('((?:\\\\.|[^\\\\'])*)'|\"((?:\\\\.|[^\\\\\"])*)\")|" + - - // 2. simple (capture 6) - "((?:\\\\.|[^\\\\()[\\]]|" + attributes + ")*)|" + - - // 3. anything else (capture 2) - ".*" + - ")\\)|)", - - // Leading and non-escaped trailing whitespace, capturing some non-whitespace characters preceding the latter - rwhitespace = new RegExp( whitespace + "+", "g" ), - rtrim = new RegExp( "^" + whitespace + "+|((?:^|[^\\\\])(?:\\\\.)*)" + - whitespace + "+$", "g" ), - - rcomma = new RegExp( "^" + whitespace + "*," + whitespace + "*" ), - rcombinators = new RegExp( "^" + whitespace + "*([>+~]|" + whitespace + ")" + whitespace + - "*" ), - rdescend = new RegExp( whitespace + "|>" ), - - rpseudo = new RegExp( pseudos ), - ridentifier = new RegExp( "^" + identifier + "$" ), - - matchExpr = { - "ID": new RegExp( "^#(" + identifier + ")" ), - "CLASS": new RegExp( "^\\.(" + identifier + ")" ), - "TAG": new RegExp( "^(" + identifier + "|[*])" ), - "ATTR": new RegExp( "^" + attributes ), - "PSEUDO": new RegExp( "^" + pseudos ), - "CHILD": new RegExp( "^:(only|first|last|nth|nth-last)-(child|of-type)(?:\\(" + - whitespace + "*(even|odd|(([+-]|)(\\d*)n|)" + whitespace + "*(?:([+-]|)" + - whitespace + "*(\\d+)|))" + whitespace + "*\\)|)", "i" ), - "bool": new RegExp( "^(?:" + booleans + ")$", "i" ), - - // For use in libraries implementing .is() - // We use this for POS matching in `select` - "needsContext": new RegExp( "^" + whitespace + - "*[>+~]|:(even|odd|eq|gt|lt|nth|first|last)(?:\\(" + whitespace + - "*((?:-\\d)?\\d*)" + whitespace + "*\\)|)(?=[^-]|$)", "i" ) - }, - - rhtml = /HTML$/i, - rinputs = /^(?:input|select|textarea|button)$/i, - rheader = /^h\d$/i, - - rnative = /^[^{]+\{\s*\[native \w/, - - // Easily-parseable/retrievable ID or TAG or CLASS selectors - rquickExpr = /^(?:#([\w-]+)|(\w+)|\.([\w-]+))$/, - - rsibling = /[+~]/, - - // CSS escapes - // http://www.w3.org/TR/CSS21/syndata.html#escaped-characters - runescape = new RegExp( "\\\\[\\da-fA-F]{1,6}" + whitespace + "?|\\\\([^\\r\\n\\f])", "g" ), - funescape = function( escape, nonHex ) { - var high = "0x" + escape.slice( 1 ) - 0x10000; - - return nonHex ? - - // Strip the backslash prefix from a non-hex escape sequence - nonHex : - - // Replace a hexadecimal escape sequence with the encoded Unicode code point - // Support: IE <=11+ - // For values outside the Basic Multilingual Plane (BMP), manually construct a - // surrogate pair - high < 0 ? - String.fromCharCode( high + 0x10000 ) : - String.fromCharCode( high >> 10 | 0xD800, high & 0x3FF | 0xDC00 ); - }, - - // CSS string/identifier serialization - // https://drafts.csswg.org/cssom/#common-serializing-idioms - rcssescape = /([\0-\x1f\x7f]|^-?\d)|^-$|[^\0-\x1f\x7f-\uFFFF\w-]/g, - fcssescape = function( ch, asCodePoint ) { - if ( asCodePoint ) { - - // U+0000 NULL becomes U+FFFD REPLACEMENT CHARACTER - if ( ch === "\0" ) { - return "\uFFFD"; - } - - // Control characters and (dependent upon position) numbers get escaped as code points - return ch.slice( 0, -1 ) + "\\" + - ch.charCodeAt( ch.length - 1 ).toString( 16 ) + " "; - } - - // Other potentially-special ASCII characters get backslash-escaped - return "\\" + ch; - }, - - // Used for iframes - // See setDocument() - // Removing the function wrapper causes a "Permission Denied" - // error in IE - unloadHandler = function() { - setDocument(); - }, - - inDisabledFieldset = addCombinator( - function( elem ) { - return elem.disabled === true && elem.nodeName.toLowerCase() === "fieldset"; - }, - { dir: "parentNode", next: "legend" } - ); - -// Optimize for push.apply( _, NodeList ) -try { - push.apply( - ( arr = slice.call( preferredDoc.childNodes ) ), - preferredDoc.childNodes - ); - - // Support: Android<4.0 - // Detect silently failing push.apply - // eslint-disable-next-line no-unused-expressions - arr[ preferredDoc.childNodes.length ].nodeType; -} catch ( e ) { - push = { apply: arr.length ? - - // Leverage slice if possible - function( target, els ) { - pushNative.apply( target, slice.call( els ) ); - } : - - // Support: IE<9 - // Otherwise append directly - function( target, els ) { - var j = target.length, - i = 0; - - // Can't trust NodeList.length - while ( ( target[ j++ ] = els[ i++ ] ) ) {} - target.length = j - 1; - } - }; -} - -function Sizzle( selector, context, results, seed ) { - var m, i, elem, nid, match, groups, newSelector, - newContext = context && context.ownerDocument, - - // nodeType defaults to 9, since context defaults to document - nodeType = context ? context.nodeType : 9; - - results = results || []; - - // Return early from calls with invalid selector or context - if ( typeof selector !== "string" || !selector || - nodeType !== 1 && nodeType !== 9 && nodeType !== 11 ) { - - return results; - } - - // Try to shortcut find operations (as opposed to filters) in HTML documents - if ( !seed ) { - setDocument( context ); - context = context || document; - - if ( documentIsHTML ) { - - // If the selector is sufficiently simple, try using a "get*By*" DOM method - // (excepting DocumentFragment context, where the methods don't exist) - if ( nodeType !== 11 && ( match = rquickExpr.exec( selector ) ) ) { - - // ID selector - if ( ( m = match[ 1 ] ) ) { - - // Document context - if ( nodeType === 9 ) { - if ( ( elem = context.getElementById( m ) ) ) { - - // Support: IE, Opera, Webkit - // TODO: identify versions - // getElementById can match elements by name instead of ID - if ( elem.id === m ) { - results.push( elem ); - return results; - } - } else { - return results; - } - - // Element context - } else { - - // Support: IE, Opera, Webkit - // TODO: identify versions - // getElementById can match elements by name instead of ID - if ( newContext && ( elem = newContext.getElementById( m ) ) && - contains( context, elem ) && - elem.id === m ) { - - results.push( elem ); - return results; - } - } - - // Type selector - } else if ( match[ 2 ] ) { - push.apply( results, context.getElementsByTagName( selector ) ); - return results; - - // Class selector - } else if ( ( m = match[ 3 ] ) && support.getElementsByClassName && - context.getElementsByClassName ) { - - push.apply( results, context.getElementsByClassName( m ) ); - return results; - } - } - - // Take advantage of querySelectorAll - if ( support.qsa && - !nonnativeSelectorCache[ selector + " " ] && - ( !rbuggyQSA || !rbuggyQSA.test( selector ) ) && - - // Support: IE 8 only - // Exclude object elements - ( nodeType !== 1 || context.nodeName.toLowerCase() !== "object" ) ) { - - newSelector = selector; - newContext = context; - - // qSA considers elements outside a scoping root when evaluating child or - // descendant combinators, which is not what we want. - // In such cases, we work around the behavior by prefixing every selector in the - // list with an ID selector referencing the scope context. - // The technique has to be used as well when a leading combinator is used - // as such selectors are not recognized by querySelectorAll. - // Thanks to Andrew Dupont for this technique. - if ( nodeType === 1 && - ( rdescend.test( selector ) || rcombinators.test( selector ) ) ) { - - // Expand context for sibling selectors - newContext = rsibling.test( selector ) && testContext( context.parentNode ) || - context; - - // We can use :scope instead of the ID hack if the browser - // supports it & if we're not changing the context. - if ( newContext !== context || !support.scope ) { - - // Capture the context ID, setting it first if necessary - if ( ( nid = context.getAttribute( "id" ) ) ) { - nid = nid.replace( rcssescape, fcssescape ); - } else { - context.setAttribute( "id", ( nid = expando ) ); - } - } - - // Prefix every selector in the list - groups = tokenize( selector ); - i = groups.length; - while ( i-- ) { - groups[ i ] = ( nid ? "#" + nid : ":scope" ) + " " + - toSelector( groups[ i ] ); - } - newSelector = groups.join( "," ); - } - - try { - - // `qSA` may not throw for unrecognized parts using forgiving parsing: - // https://drafts.csswg.org/selectors/#forgiving-selector - // like the `:has()` pseudo-class: - // https://drafts.csswg.org/selectors/#relational - // `CSS.supports` is still expected to return `false` then: - // https://drafts.csswg.org/css-conditional-4/#typedef-supports-selector-fn - // https://drafts.csswg.org/css-conditional-4/#dfn-support-selector - if ( support.cssSupportsSelector && - - // eslint-disable-next-line no-undef - !CSS.supports( "selector(:is(" + newSelector + "))" ) ) { - - // Support: IE 11+ - // Throw to get to the same code path as an error directly in qSA. - // Note: once we only support browser supporting - // `CSS.supports('selector(...)')`, we can most likely drop - // the `try-catch`. IE doesn't implement the API. - throw new Error(); - } - - push.apply( results, - newContext.querySelectorAll( newSelector ) - ); - return results; - } catch ( qsaError ) { - nonnativeSelectorCache( selector, true ); - } finally { - if ( nid === expando ) { - context.removeAttribute( "id" ); - } - } - } - } - } - - // All others - return select( selector.replace( rtrim, "$1" ), context, results, seed ); -} - -/** - * Create key-value caches of limited size - * @returns {function(string, object)} Returns the Object data after storing it on itself with - * property name the (space-suffixed) string and (if the cache is larger than Expr.cacheLength) - * deleting the oldest entry - */ -function createCache() { - var keys = []; - - function cache( key, value ) { - - // Use (key + " ") to avoid collision with native prototype properties (see Issue #157) - if ( keys.push( key + " " ) > Expr.cacheLength ) { - - // Only keep the most recent entries - delete cache[ keys.shift() ]; - } - return ( cache[ key + " " ] = value ); - } - return cache; -} - -/** - * Mark a function for special use by Sizzle - * @param {Function} fn The function to mark - */ -function markFunction( fn ) { - fn[ expando ] = true; - return fn; -} - -/** - * Support testing using an element - * @param {Function} fn Passed the created element and returns a boolean result - */ -function assert( fn ) { - var el = document.createElement( "fieldset" ); - - try { - return !!fn( el ); - } catch ( e ) { - return false; - } finally { - - // Remove from its parent by default - if ( el.parentNode ) { - el.parentNode.removeChild( el ); - } - - // release memory in IE - el = null; - } -} - -/** - * Adds the same handler for all of the specified attrs - * @param {String} attrs Pipe-separated list of attributes - * @param {Function} handler The method that will be applied - */ -function addHandle( attrs, handler ) { - var arr = attrs.split( "|" ), - i = arr.length; - - while ( i-- ) { - Expr.attrHandle[ arr[ i ] ] = handler; - } -} - -/** - * Checks document order of two siblings - * @param {Element} a - * @param {Element} b - * @returns {Number} Returns less than 0 if a precedes b, greater than 0 if a follows b - */ -function siblingCheck( a, b ) { - var cur = b && a, - diff = cur && a.nodeType === 1 && b.nodeType === 1 && - a.sourceIndex - b.sourceIndex; - - // Use IE sourceIndex if available on both nodes - if ( diff ) { - return diff; - } - - // Check if b follows a - if ( cur ) { - while ( ( cur = cur.nextSibling ) ) { - if ( cur === b ) { - return -1; - } - } - } - - return a ? 1 : -1; -} - -/** - * Returns a function to use in pseudos for input types - * @param {String} type - */ -function createInputPseudo( type ) { - return function( elem ) { - var name = elem.nodeName.toLowerCase(); - return name === "input" && elem.type === type; - }; -} - -/** - * Returns a function to use in pseudos for buttons - * @param {String} type - */ -function createButtonPseudo( type ) { - return function( elem ) { - var name = elem.nodeName.toLowerCase(); - return ( name === "input" || name === "button" ) && elem.type === type; - }; -} - -/** - * Returns a function to use in pseudos for :enabled/:disabled - * @param {Boolean} disabled true for :disabled; false for :enabled - */ -function createDisabledPseudo( disabled ) { - - // Known :disabled false positives: fieldset[disabled] > legend:nth-of-type(n+2) :can-disable - return function( elem ) { - - // Only certain elements can match :enabled or :disabled - // https://html.spec.whatwg.org/multipage/scripting.html#selector-enabled - // https://html.spec.whatwg.org/multipage/scripting.html#selector-disabled - if ( "form" in elem ) { - - // Check for inherited disabledness on relevant non-disabled elements: - // * listed form-associated elements in a disabled fieldset - // https://html.spec.whatwg.org/multipage/forms.html#category-listed - // https://html.spec.whatwg.org/multipage/forms.html#concept-fe-disabled - // * option elements in a disabled optgroup - // https://html.spec.whatwg.org/multipage/forms.html#concept-option-disabled - // All such elements have a "form" property. - if ( elem.parentNode && elem.disabled === false ) { - - // Option elements defer to a parent optgroup if present - if ( "label" in elem ) { - if ( "label" in elem.parentNode ) { - return elem.parentNode.disabled === disabled; - } else { - return elem.disabled === disabled; - } - } - - // Support: IE 6 - 11 - // Use the isDisabled shortcut property to check for disabled fieldset ancestors - return elem.isDisabled === disabled || - - // Where there is no isDisabled, check manually - /* jshint -W018 */ - elem.isDisabled !== !disabled && - inDisabledFieldset( elem ) === disabled; - } - - return elem.disabled === disabled; - - // Try to winnow out elements that can't be disabled before trusting the disabled property. - // Some victims get caught in our net (label, legend, menu, track), but it shouldn't - // even exist on them, let alone have a boolean value. - } else if ( "label" in elem ) { - return elem.disabled === disabled; - } - - // Remaining elements are neither :enabled nor :disabled - return false; - }; -} - -/** - * Returns a function to use in pseudos for positionals - * @param {Function} fn - */ -function createPositionalPseudo( fn ) { - return markFunction( function( argument ) { - argument = +argument; - return markFunction( function( seed, matches ) { - var j, - matchIndexes = fn( [], seed.length, argument ), - i = matchIndexes.length; - - // Match elements found at the specified indexes - while ( i-- ) { - if ( seed[ ( j = matchIndexes[ i ] ) ] ) { - seed[ j ] = !( matches[ j ] = seed[ j ] ); - } - } - } ); - } ); -} - -/** - * Checks a node for validity as a Sizzle context - * @param {Element|Object=} context - * @returns {Element|Object|Boolean} The input node if acceptable, otherwise a falsy value - */ -function testContext( context ) { - return context && typeof context.getElementsByTagName !== "undefined" && context; -} - -// Expose support vars for convenience -support = Sizzle.support = {}; - -/** - * Detects XML nodes - * @param {Element|Object} elem An element or a document - * @returns {Boolean} True iff elem is a non-HTML XML node - */ -isXML = Sizzle.isXML = function( elem ) { - var namespace = elem && elem.namespaceURI, - docElem = elem && ( elem.ownerDocument || elem ).documentElement; - - // Support: IE <=8 - // Assume HTML when documentElement doesn't yet exist, such as inside loading iframes - // https://bugs.jquery.com/ticket/4833 - return !rhtml.test( namespace || docElem && docElem.nodeName || "HTML" ); -}; - -/** - * Sets document-related variables once based on the current document - * @param {Element|Object} [doc] An element or document object to use to set the document - * @returns {Object} Returns the current document - */ -setDocument = Sizzle.setDocument = function( node ) { - var hasCompare, subWindow, - doc = node ? node.ownerDocument || node : preferredDoc; - - // Return early if doc is invalid or already selected - // Support: IE 11+, Edge 17 - 18+ - // IE/Edge sometimes throw a "Permission denied" error when strict-comparing - // two documents; shallow comparisons work. - // eslint-disable-next-line eqeqeq - if ( doc == document || doc.nodeType !== 9 || !doc.documentElement ) { - return document; - } - - // Update global variables - document = doc; - docElem = document.documentElement; - documentIsHTML = !isXML( document ); - - // Support: IE 9 - 11+, Edge 12 - 18+ - // Accessing iframe documents after unload throws "permission denied" errors (jQuery #13936) - // Support: IE 11+, Edge 17 - 18+ - // IE/Edge sometimes throw a "Permission denied" error when strict-comparing - // two documents; shallow comparisons work. - // eslint-disable-next-line eqeqeq - if ( preferredDoc != document && - ( subWindow = document.defaultView ) && subWindow.top !== subWindow ) { - - // Support: IE 11, Edge - if ( subWindow.addEventListener ) { - subWindow.addEventListener( "unload", unloadHandler, false ); - - // Support: IE 9 - 10 only - } else if ( subWindow.attachEvent ) { - subWindow.attachEvent( "onunload", unloadHandler ); - } - } - - // Support: IE 8 - 11+, Edge 12 - 18+, Chrome <=16 - 25 only, Firefox <=3.6 - 31 only, - // Safari 4 - 5 only, Opera <=11.6 - 12.x only - // IE/Edge & older browsers don't support the :scope pseudo-class. - // Support: Safari 6.0 only - // Safari 6.0 supports :scope but it's an alias of :root there. - support.scope = assert( function( el ) { - docElem.appendChild( el ).appendChild( document.createElement( "div" ) ); - return typeof el.querySelectorAll !== "undefined" && - !el.querySelectorAll( ":scope fieldset div" ).length; - } ); - - // Support: Chrome 105+, Firefox 104+, Safari 15.4+ - // Make sure forgiving mode is not used in `CSS.supports( "selector(...)" )`. - // - // `:is()` uses a forgiving selector list as an argument and is widely - // implemented, so it's a good one to test against. - support.cssSupportsSelector = assert( function() { - /* eslint-disable no-undef */ - - return CSS.supports( "selector(*)" ) && - - // Support: Firefox 78-81 only - // In old Firefox, `:is()` didn't use forgiving parsing. In that case, - // fail this test as there's no selector to test against that. - // `CSS.supports` uses unforgiving parsing - document.querySelectorAll( ":is(:jqfake)" ) && - - // `*` is needed as Safari & newer Chrome implemented something in between - // for `:has()` - it throws in `qSA` if it only contains an unsupported - // argument but multiple ones, one of which is supported, are fine. - // We want to play safe in case `:is()` gets the same treatment. - !CSS.supports( "selector(:is(*,:jqfake))" ); - - /* eslint-enable */ - } ); - - /* Attributes - ---------------------------------------------------------------------- */ - - // Support: IE<8 - // Verify that getAttribute really returns attributes and not properties - // (excepting IE8 booleans) - support.attributes = assert( function( el ) { - el.className = "i"; - return !el.getAttribute( "className" ); - } ); - - /* getElement(s)By* - ---------------------------------------------------------------------- */ - - // Check if getElementsByTagName("*") returns only elements - support.getElementsByTagName = assert( function( el ) { - el.appendChild( document.createComment( "" ) ); - return !el.getElementsByTagName( "*" ).length; - } ); - - // Support: IE<9 - support.getElementsByClassName = rnative.test( document.getElementsByClassName ); - - // Support: IE<10 - // Check if getElementById returns elements by name - // The broken getElementById methods don't pick up programmatically-set names, - // so use a roundabout getElementsByName test - support.getById = assert( function( el ) { - docElem.appendChild( el ).id = expando; - return !document.getElementsByName || !document.getElementsByName( expando ).length; - } ); - - // ID filter and find - if ( support.getById ) { - Expr.filter[ "ID" ] = function( id ) { - var attrId = id.replace( runescape, funescape ); - return function( elem ) { - return elem.getAttribute( "id" ) === attrId; - }; - }; - Expr.find[ "ID" ] = function( id, context ) { - if ( typeof context.getElementById !== "undefined" && documentIsHTML ) { - var elem = context.getElementById( id ); - return elem ? [ elem ] : []; - } - }; - } else { - Expr.filter[ "ID" ] = function( id ) { - var attrId = id.replace( runescape, funescape ); - return function( elem ) { - var node = typeof elem.getAttributeNode !== "undefined" && - elem.getAttributeNode( "id" ); - return node && node.value === attrId; - }; - }; - - // Support: IE 6 - 7 only - // getElementById is not reliable as a find shortcut - Expr.find[ "ID" ] = function( id, context ) { - if ( typeof context.getElementById !== "undefined" && documentIsHTML ) { - var node, i, elems, - elem = context.getElementById( id ); - - if ( elem ) { - - // Verify the id attribute - node = elem.getAttributeNode( "id" ); - if ( node && node.value === id ) { - return [ elem ]; - } - - // Fall back on getElementsByName - elems = context.getElementsByName( id ); - i = 0; - while ( ( elem = elems[ i++ ] ) ) { - node = elem.getAttributeNode( "id" ); - if ( node && node.value === id ) { - return [ elem ]; - } - } - } - - return []; - } - }; - } - - // Tag - Expr.find[ "TAG" ] = support.getElementsByTagName ? - function( tag, context ) { - if ( typeof context.getElementsByTagName !== "undefined" ) { - return context.getElementsByTagName( tag ); - - // DocumentFragment nodes don't have gEBTN - } else if ( support.qsa ) { - return context.querySelectorAll( tag ); - } - } : - - function( tag, context ) { - var elem, - tmp = [], - i = 0, - - // By happy coincidence, a (broken) gEBTN appears on DocumentFragment nodes too - results = context.getElementsByTagName( tag ); - - // Filter out possible comments - if ( tag === "*" ) { - while ( ( elem = results[ i++ ] ) ) { - if ( elem.nodeType === 1 ) { - tmp.push( elem ); - } - } - - return tmp; - } - return results; - }; - - // Class - Expr.find[ "CLASS" ] = support.getElementsByClassName && function( className, context ) { - if ( typeof context.getElementsByClassName !== "undefined" && documentIsHTML ) { - return context.getElementsByClassName( className ); - } - }; - - /* QSA/matchesSelector - ---------------------------------------------------------------------- */ - - // QSA and matchesSelector support - - // matchesSelector(:active) reports false when true (IE9/Opera 11.5) - rbuggyMatches = []; - - // qSa(:focus) reports false when true (Chrome 21) - // We allow this because of a bug in IE8/9 that throws an error - // whenever `document.activeElement` is accessed on an iframe - // So, we allow :focus to pass through QSA all the time to avoid the IE error - // See https://bugs.jquery.com/ticket/13378 - rbuggyQSA = []; - - if ( ( support.qsa = rnative.test( document.querySelectorAll ) ) ) { - - // Build QSA regex - // Regex strategy adopted from Diego Perini - assert( function( el ) { - - var input; - - // Select is set to empty string on purpose - // This is to test IE's treatment of not explicitly - // setting a boolean content attribute, - // since its presence should be enough - // https://bugs.jquery.com/ticket/12359 - docElem.appendChild( el ).innerHTML = "" + - ""; - - // Support: IE8, Opera 11-12.16 - // Nothing should be selected when empty strings follow ^= or $= or *= - // The test attribute must be unknown in Opera but "safe" for WinRT - // https://msdn.microsoft.com/en-us/library/ie/hh465388.aspx#attribute_section - if ( el.querySelectorAll( "[msallowcapture^='']" ).length ) { - rbuggyQSA.push( "[*^$]=" + whitespace + "*(?:''|\"\")" ); - } - - // Support: IE8 - // Boolean attributes and "value" are not treated correctly - if ( !el.querySelectorAll( "[selected]" ).length ) { - rbuggyQSA.push( "\\[" + whitespace + "*(?:value|" + booleans + ")" ); - } - - // Support: Chrome<29, Android<4.4, Safari<7.0+, iOS<7.0+, PhantomJS<1.9.8+ - if ( !el.querySelectorAll( "[id~=" + expando + "-]" ).length ) { - rbuggyQSA.push( "~=" ); - } - - // Support: IE 11+, Edge 15 - 18+ - // IE 11/Edge don't find elements on a `[name='']` query in some cases. - // Adding a temporary attribute to the document before the selection works - // around the issue. - // Interestingly, IE 10 & older don't seem to have the issue. - input = document.createElement( "input" ); - input.setAttribute( "name", "" ); - el.appendChild( input ); - if ( !el.querySelectorAll( "[name='']" ).length ) { - rbuggyQSA.push( "\\[" + whitespace + "*name" + whitespace + "*=" + - whitespace + "*(?:''|\"\")" ); - } - - // Webkit/Opera - :checked should return selected option elements - // http://www.w3.org/TR/2011/REC-css3-selectors-20110929/#checked - // IE8 throws error here and will not see later tests - if ( !el.querySelectorAll( ":checked" ).length ) { - rbuggyQSA.push( ":checked" ); - } - - // Support: Safari 8+, iOS 8+ - // https://bugs.webkit.org/show_bug.cgi?id=136851 - // In-page `selector#id sibling-combinator selector` fails - if ( !el.querySelectorAll( "a#" + expando + "+*" ).length ) { - rbuggyQSA.push( ".#.+[+~]" ); - } - - // Support: Firefox <=3.6 - 5 only - // Old Firefox doesn't throw on a badly-escaped identifier. - el.querySelectorAll( "\\\f" ); - rbuggyQSA.push( "[\\r\\n\\f]" ); - } ); - - assert( function( el ) { - el.innerHTML = "" + - ""; - - // Support: Windows 8 Native Apps - // The type and name attributes are restricted during .innerHTML assignment - var input = document.createElement( "input" ); - input.setAttribute( "type", "hidden" ); - el.appendChild( input ).setAttribute( "name", "D" ); - - // Support: IE8 - // Enforce case-sensitivity of name attribute - if ( el.querySelectorAll( "[name=d]" ).length ) { - rbuggyQSA.push( "name" + whitespace + "*[*^$|!~]?=" ); - } - - // FF 3.5 - :enabled/:disabled and hidden elements (hidden elements are still enabled) - // IE8 throws error here and will not see later tests - if ( el.querySelectorAll( ":enabled" ).length !== 2 ) { - rbuggyQSA.push( ":enabled", ":disabled" ); - } - - // Support: IE9-11+ - // IE's :disabled selector does not pick up the children of disabled fieldsets - docElem.appendChild( el ).disabled = true; - if ( el.querySelectorAll( ":disabled" ).length !== 2 ) { - rbuggyQSA.push( ":enabled", ":disabled" ); - } - - // Support: Opera 10 - 11 only - // Opera 10-11 does not throw on post-comma invalid pseudos - el.querySelectorAll( "*,:x" ); - rbuggyQSA.push( ",.*:" ); - } ); - } - - if ( ( support.matchesSelector = rnative.test( ( matches = docElem.matches || - docElem.webkitMatchesSelector || - docElem.mozMatchesSelector || - docElem.oMatchesSelector || - docElem.msMatchesSelector ) ) ) ) { - - assert( function( el ) { - - // Check to see if it's possible to do matchesSelector - // on a disconnected node (IE 9) - support.disconnectedMatch = matches.call( el, "*" ); - - // This should fail with an exception - // Gecko does not error, returns false instead - matches.call( el, "[s!='']:x" ); - rbuggyMatches.push( "!=", pseudos ); - } ); - } - - if ( !support.cssSupportsSelector ) { - - // Support: Chrome 105+, Safari 15.4+ - // `:has()` uses a forgiving selector list as an argument so our regular - // `try-catch` mechanism fails to catch `:has()` with arguments not supported - // natively like `:has(:contains("Foo"))`. Where supported & spec-compliant, - // we now use `CSS.supports("selector(:is(SELECTOR_TO_BE_TESTED))")`, but - // outside that we mark `:has` as buggy. - rbuggyQSA.push( ":has" ); - } - - rbuggyQSA = rbuggyQSA.length && new RegExp( rbuggyQSA.join( "|" ) ); - rbuggyMatches = rbuggyMatches.length && new RegExp( rbuggyMatches.join( "|" ) ); - - /* Contains - ---------------------------------------------------------------------- */ - hasCompare = rnative.test( docElem.compareDocumentPosition ); - - // Element contains another - // Purposefully self-exclusive - // As in, an element does not contain itself - contains = hasCompare || rnative.test( docElem.contains ) ? - function( a, b ) { - - // Support: IE <9 only - // IE doesn't have `contains` on `document` so we need to check for - // `documentElement` presence. - // We need to fall back to `a` when `documentElement` is missing - // as `ownerDocument` of elements within `