/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *//* Geodesy tools for an ellipsoidal earth model (c) Chris Veness 2005-2016 *//* MIT Licence *//* www.movable-type.co.uk/scripts/latlong-convert-coords.html *//* www.movable-type.co.uk/scripts/geodesy/docs/module-latlon-ellipsoidal.html *//* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */import Dms from './dms.js';/** * A latitude/longitude point defines a geographic location on or relative to the earth’s surface, * measured in degrees from the equator and the Greenwich meridian and metres above the ellipsoid, * and based on a given datum. The datum defines the reference ellipsoid, and the transformation * parameters between different datums. * * The module includes ellipsoid parameters and datums for different coordinate systems, and methods * for converting between them and to/from cartesian coordinates. * * 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. * * This module is used for both trigonometric methods (eg latlon-vincenty) and n-vector methods (eg * latlon-nvector-ellipsoidal). * * @module latlon-ellipsoidal * @requires dms *//* * Ellipsoid parameters; semi-major axis (a), semi-minor axis (b), and flattening (f) for each ellipsoid. * Flattening f = (a−b)/a; at least one of these parameters is derived from defining constants. */const ellipsoids = { WGS84: { a: 6378137, b: 6356752.314245, f: 1/298.257223563 }, GRS80: { a: 6378137, b: 6356752.314140, f: 1/298.257222101 }, 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.962818, f: 1/299.1528128 }, Clarke1866: { a: 6378206.4, b: 6356583.8, f: 1/294.978698214 }, Intl1924: { a: 6378388, b: 6356911.946, f: 1/297 }, // aka Hayford WGS72: { a: 6378135, b: 6356750.5, f: 1/298.26 },};/* * Datums; with associated ellipsoid, and Helmert transform parameters to convert from WGS-84 into * given datum. * * More are available from earth-info.nga.mil/GandG/coordsys/datums/NATO_DT.pdf, * www.fieldenmaps.info/cconv/web/cconv_params.js */const datums = { /* eslint key-spacing: 0, comma-dangle: 0 */ WGS84: { ellipsoid: ellipsoids.WGS84, transform: { tx: 0.0, ty: 0.0, tz: 0.0, // m rx: 0.0, ry: 0.0, rz: 0.0, // sec s: 0.0 } // ppm }, NAD83: { // (2009); functionally ≡ WGS84 - www.uvm.edu/giv/resources/WGS84_NAD83.pdf ellipsoid: ellipsoids.GRS80, transform: { tx: 1.004, ty: -1.910, tz: -0.515, // m rx: 0.0267, ry: 0.00034, rz: 0.011, // sec s: -0.0015 } // ppm }, // note: if you *really* need to convert WGS84<->NAD83, you need more knowledge than this! OSGB36: { // www.ordnancesurvey.co.uk/docs/support/guide-coordinate-systems-great-britain.pdf ellipsoid: ellipsoids.Airy1830, transform: { tx: -446.448, ty: 125.157, tz: -542.060, // m rx: -0.1502, ry: -0.2470, rz: -0.8421, // sec s: 20.4894 } // ppm }, ED50: { // og.decc.gov.uk/en/olgs/cms/pons_and_cop/pons/pon4/pon4.aspx ellipsoid: ellipsoids.Intl1924, transform: { tx: 89.5, ty: 93.8, tz: 123.1, // m rx: 0.0, ry: 0.0, rz: 0.156, // sec s: -1.2 } // ppm }, Irl1975: { // osi.ie/OSI/media/OSI/Content/Publications/transformations_booklet.pdf ellipsoid: ellipsoids.AiryModified, transform: { tx: -482.530, ty: 130.596, tz: -564.557, // m rx: -1.042, ry: -0.214, rz: -0.631, // sec s: -8.150 } // ppm }, // TODO: many sources have opposite sign to rotations - to be checked! TokyoJapan: { // www.geocachingtoolbox.com?page=datumEllipsoidDetails ellipsoid: ellipsoids.Bessel1841, transform: { tx: 148, ty: -507, tz: -685, // m rx: 0, ry: 0, rz: 0, // sec s: 0 } // ppm }, NAD27: { // en.wikipedia.org/wiki/Helmert_transformation ellipsoid: ellipsoids.Clarke1866, transform: { tx: 8, ty: -160, tz: -176, // m rx: 0, ry: 0, rz: 0, // sec s: 0 } // ppm }, WGS72: { // www.icao.int/safety/pbn/documentation/eurocontrol/eurocontrol wgs 84 implementation manual.pdf ellipsoid: ellipsoids.WGS72, transform: { tx: 0, ty: 0, tz: -4.5, // m rx: 0, ry: 0, rz: 0.554, // sec s: -0.22 } // ppm },};/* LatLonEllipsoidal - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *//** * Latitude/longitude points on an ellipsoidal model earth, with ellipsoid parameters and methods * for converting between datums and to cartesian (ECEF) coordinates. */class LatLonEllipsoidal { // note prototype-based class not inheritance-based class /** * Creates lat/lon (polar) point with latitude & longitude values, on a specified datum. * * @param {number} lat - Geodetic latitude in degrees. * @param {number} lon - Longitude in degrees. * @param {number} [height=0] - Height above ellipsoid in metres. * @param {LatLon.datums} [datum=WGS84] - Datum this point is defined within. * * @example * import LatLon from 'latlon-ellipsoidal'; * var p1 = new LatLon(51.4778, -0.0016, 0, LatLon.datums.WGS84); */ constructor(lat, lon, height=0, datum=datums.WGS84) { this.lat = Number(lat); this.lon = Number(lon); this.height = Number(height); this.datum = datum; } /** * Ellipsoid parameters; major axis (a), minor axis (b), and flattening (f) for each ellipsoid. * * @example * var 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. * * @example * var a = LatLon.datums.OSGB36.ellipsoid.a; // 6377563.396 * var tx = LatLon.datums.OSGB36.transform; // to WGS-84 */ static get datums() { return datums; } /** * 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. * * @example * var pWGS84 = new LatLon(51.4778, -0.0016, 0, LatLon.datums.WGS84); * var pOSGB = pWGS84.convertDatum(LatLon.datums.OSGB36); // 51.4773°N, 000.0000°E */ convertDatum(toDatum) { var oldLatLon = this; var transform; if (oldLatLon.datum == datums.WGS84) { // converting from WGS 84 transform = toDatum.transform; } if (toDatum == datums.WGS84) { // converting to WGS 84; use inverse transform (don't overwrite original!) transform = {}; for (var param in oldLatLon.datum.transform) { if (oldLatLon.datum.transform.hasOwnProperty(param)) { transform[param] = -oldLatLon.datum.transform[param]; } } } if (transform === undefined) { // neither this.datum nor toDatum are WGS84: convert this to WGS84 first oldLatLon = this.convertDatum(datums.WGS84); transform = toDatum.transform; } var oldCartesian = oldLatLon.toCartesian(); // convert polar to cartesian... var newCartesian = oldCartesian.applyTransform(transform); // ...apply transform... var newLatLon = newCartesian.toLatLon(toDatum); // ...and convert cartesian to polar return newLatLon; } /** * 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() { var φ = this.lat.toRadians(), λ = this.lon.toRadians(); var h = this.height; // height above ellipsoid var a = this.datum.ellipsoid.a, f = this.datum.ellipsoid.f; var sinφ = Math.sin(φ), cosφ = Math.cos(φ); var sinλ = Math.sin(λ), cosλ = Math.cos(λ); var eSq = 2*f - f*f; // 1st eccentricity squared ≡ (a²-b²)/a² var ν = a / Math.sqrt(1 - eSq*sinφ*sinφ); // radius of curvature in prime vertical var x = (ν+h) * cosφ * cosλ; var y = (ν+h) * cosφ * sinλ; var z = (ν*(1-eSq)+h) * sinφ; var p = new Cartesian(x, y, z); return p; } /** * Returns a string representation of ‘this’ point, formatted as degrees, degrees+minutes, or * degrees+minutes+seconds. * * @param {string} [format=dms] - Format point as 'd', 'dm', 'dms'. * @param {number} [dp=0|2|4] - Number of decimal places to use: default 0 for dms, 2 for dm, 4 for d. * @param {number} [heightDp=null] - Number of decimal places to use for height; default is no height display. * @returns {string} Comma-separated formatted latitude/longitude. */ toString(format='dms', dp=undefined, heightDp=null) { var height = heightDp==null ? '' : ' ' + (this.height>0 ? '+' : '') + this.height.toFixed(Number(heightDp)) + 'm'; return Dms.toLat(this.lat, format, dp) + ', ' + Dms.toLon(this.lon, format, dp) + height; }}/* Cartesian - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *//** * Converts ECEF (earth-centered earth-fixed) cartesian coordinates to LatLon points, applies * Helmert transformations. */class Cartesian { // note prototype-based class not inheritance-based class /** * 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 'latlon-ellipsoidal'; * var coord = new Cartesian(3980581.210, -111.159, 4966824.522); */ constructor(x, y, z) { this.x = Number(x); this.y = Number(y); this.z = Number(z); } /** * Converts ‘this’ (geocentric) cartesian (x/y/z) coordinate to (ellipsoidal geodetic) * latitude/longitude point on specified datum. * * Uses Bowring’s (1985) formulation for μm precision in concise form. * * @param {LatLon.datums} [datum=WGS84] - Datum to use when converting point. */ toLatLon(datum=LatLonEllipsoidal.datums.WGS84) { var x = this.x, y = this.y, z = this.z; var a = datum.ellipsoid.a, b = datum.ellipsoid.b, f = datum.ellipsoid.f; var e2 = 2*f - f*f; // 1st eccentricity squared ≡ (a²-b²)/a² var ε2 = e2 / (1-e2); // 2nd eccentricity squared ≡ (a²-b²)/b² var p = Math.sqrt(x*x + y*y); // distance from minor axis var R = Math.sqrt(p*p + z*z); // polar radius // parametric latitude (Bowring eqn 17, replacing tanβ = z·a / p·b) var tanβ = (b*z)/(a*p) * (1+ε2*b/R); var sinβ = tanβ / Math.sqrt(1+tanβ*tanβ); var cosβ = sinβ / tanβ; // geodetic latitude (Bowring eqn 18: tanφ = z+ε²bsin³β / p−e²cos³β) var φ = isNaN(cosβ) ? 0 : Math.atan2(z + ε2*b*sinβ*sinβ*sinβ, p - e2*a*cosβ*cosβ*cosβ); // longitude var λ = Math.atan2(y, x); // height above ellipsoid (Bowring eqn 7) var sinφ = Math.sin(φ), cosφ = Math.cos(φ); var ν = a / Math.sqrt(1-e2*sinφ*sinφ); // length of the normal terminated by the minor axis var h = p*cosφ + z*sinφ - (a*a/ν); var point = new LatLonEllipsoidal(φ.toDegrees(), λ.toDegrees(), h, datum); return point; } /** * Applies Helmert (seven-parameter) transformation to ‘this’ coordinate using transform * parameters t. * * @param {LatLon.datums.transform} t - Transformation to apply to this coordinate. */ applyTransform(t) { var x1 = this.x, y1 = this.y, z1 = this.z; var tx = t.tx, ty = t.ty, tz = t.tz; var rx = (t.rx/3600).toRadians(); // normalise seconds to radians var ry = (t.ry/3600).toRadians(); // normalise seconds to radians var rz = (t.rz/3600).toRadians(); // normalise seconds to radians var s1 = t.s/1e6 + 1; // normalise ppm to (s+1) // apply transform var x2 = tx + x1*s1 - y1*rz + z1*ry; var y2 = ty + x1*rz + y1*s1 - z1*rx; var z2 = tz - x1*ry + y1*rx + z1*s1; var point = new Cartesian(x2, y2, z2); 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) { dp = Number(dp); return '['+this.x.toFixed(dp)+','+this.y.toFixed(dp)+','+this.z.toFixed(dp)+']'; }}// Extend Number object with methods to convert between degrees & radiansNumber.prototype.toRadians = function() { return this * Math.PI / 180; };Number.prototype.toDegrees = function() { return this * 180 / Math.PI; };export { LatLonEllipsoidal as default, Cartesian };/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */