geodesy - npm Package Compare versions

Comparing version 2.1.0 to 2.2.0

CHANGELOG.md

 # Changelog
 
+## [2.2.0] - 2019-07-08
+
+### Fixed
+
+- Fix vincenty inverse calculation for antipodal points
+- Provide convertDatum() method on a LatLon obtained from Utm.toLatLon()
+
+### Added
+
+- Option to override UTM zone in LatLon.toUtm(), option to suppress UTM easting/northing checks
+- ETRS89 datum (≡ WGS84 @ 1m level)
+
 ## [2.1.0] - 2019-06-03
@@ -4,0 +16,0 @@

latlon-ellipsoidal-datum.js

@@ -8,3 +8,2 @@ /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */
 
-
 import LatLonEllipsoidal, { Cartesian, Dms } from './latlon-ellipsoidal.js';
@@ -59,3 +58,3 @@
     ED50:       { ellipsoid: ellipsoids.Intl1924,      transform: [   89.5,    93.8,    123.1,    -1.2,     0.0,      0.0,      0.156    ] }, // epsg.io/1311
-    // en.wikipedia.org/wiki/European_Terrestrial_Reference_System_1989
+    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
@@ -236,6 +235,9 @@ NAD27:      { ellipsoid: ellipsoids.Clarke1866,    transform: [    8,    -160,     -176,       0,       0,        0,        0        ] },
      * Converts ‘this’ point from (geodetic) latitude/longitude coordinates to (geocentric) cartesian
-     * (x/y/z) coordinates.
+     * (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.
+     *   earth centre, augmented with reference frame conversion methods and properties.
      */
@@ -242,0 +244,0 @@ toCartesian() {

latlon-ellipsoidal-vincenty.js

@@ -11,2 +11,3 @@ /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */
 const π = Math.PI;
+const ε = Number.EPSILON;
 
@@ -179,4 +180,4 @@
         const tanU1 = (1-f) * Math.tan(φ1), cosU1 = 1 / Math.sqrt((1 + tanU1*tanU1)), sinU1 = tanU1 * cosU1;
-        const σ1 = Math.atan2(tanU1, cosα1);
-        const sinα = cosU1 * sinα1;
+        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α;
@@ -187,11 +188,12 @@ const uSq = cosSqα * (a*a - b*b) / (b*b);
 
-        let cos2σM = null, sinσ = null, cosσ = null, Δσ = null;
+        let σ = s / (b*A), sinσ = null, cosσ = null, Δσ = 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 σ = s / (b*A), σʹ, iterations = 0;
+        let σʹ = null, iterations = 0;
         do {
-            cos2σM = Math.cos(2*σ1 + σ);
+            cos2σₘ = Math.cos(2*σ1 + σ);
             sinσ = Math.sin(σ);
             cosσ = Math.cos(σ);
-            Δσ = B*sinσ*(cos2σM+B/4*(cosσ*(-1+2*cos2σM*cos2σM)-
-                B/6*cos2σM*(-3+4*sinσ*sinσ)*(-3+4*cos2σM*cos2σM)));
+            Δσ = B*sinσ*(cos2σₘ+B/4*(cosσ*(-1+2*cos2σₘ*cos2σₘ)-
+                B/6*cos2σₘ*(-3+4*sinσ*sinσ)*(-3+4*cos2σₘ*cos2σₘ)));
             σʹ = σ;
@@ -206,3 +208,3 @@ σ = s / (b*A) + Δσ;
         const C = f/16*cosSqα*(4+f*(4-3*cosSqα));
-        const L = λ - (1-C) * f * sinα * (σ + C*sinσ*(cos2σM+C*cosσ*(-1+2*cos2σM*cos2σM)));
+        const L = λ - (1-C) * f * sinα * (σ + C*sinσ*(cos2σₘ+C*cosσ*(-1+2*cos2σₘ*cos2σₘ)));
         const λ2 = λ1 + L;
@@ -246,11 +248,15 @@
 
-        const L = λ2 - λ1;
+        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;
 
-        let sinλ = null, cosλ = null, sinσ = 0, cosσ = 0, sinα = null;
-        let sinSqσ = null, cosSqα = 0, cos2σM = 0, σ = null, C = null;
+        const antipodal = Math.abs(L) > π/2 || Math.abs(φ2-φ1) > π/2;
 
-        let λ = L, λʹ, iterations = 0;
-        const antimeridian = Math.abs(L) > π;
+        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;                      // σM = angular distance on the sphere from the equator to the midpoint of the line
+        let sinα = null, cosSqα = 1;         // α = azimuth of the geodesic at the equator
+        let C = null;
+
+        let λʹ = null, iterations = 0;
         do {
@@ -260,3 +266,3 @@ sinλ = Math.sin(λ);
             sinSqσ = (cosU2*sinλ) * (cosU2*sinλ) + (cosU1*sinU2-sinU1*cosU2*cosλ) * (cosU1*sinU2-sinU1*cosU2*cosλ);
-            if (Math.abs(sinSqσ) < Number.EPSILON) break;  // co-incident points
+            if (Math.abs(sinSqσ) < ε) break;  // co-incident/antipodal points (falls back on λ/σ = L)
             sinσ = Math.sqrt(sinSqσ);
@@ -267,7 +273,7 @@ cosσ = sinU1*sinU2 + cosU1*cosU2*cosλ;
             cosSqα = 1 - sinα*sinα;
-            cos2σM = (cosSqα != 0) ? (cosσ - 2*sinU1*sinU2/cosSqα) : 0; // on equatorial line cos²α = 0 (§6)
+            cos2σₘ = (cosSqα != 0) ? (cosσ - 2*sinU1*sinU2/cosSqα) : 0; // on equatorial line cos²α = 0 (§6)
             C = f/16*cosSqα*(4+f*(4-3*cosSqα));
             λʹ = λ;
-            λ = L + (1-C) * f * sinα * (σ + C*sinσ*(cos2σM+C*cosσ*(-1+2*cos2σM*cos2σM)));
-            const iterationCheck = antimeridian ? Math.abs(λ)-π : Math.abs(λ);
+            λ = 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('λ > π');
@@ -280,14 +286,18 @@ } while (Math.abs(λ-λʹ) > 1e-12 && ++iterations<1000);
         const B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)));
-        const Δσ = B*sinσ*(cos2σM+B/4*(cosσ*(-1+2*cos2σM*cos2σM)-
-            B/6*cos2σM*(-3+4*sinσ*sinσ)*(-3+4*cos2σM*cos2σM)));
+        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*(σ-Δσ);
+        const s = b*A*(σ-Δσ); // s = length of the geodesic
 
-        const α1 = Math.atan2(cosU2*sinλ,  cosU1*sinU2-sinU1*cosU2*cosλ);
-        const α2 = Math.atan2(cosU1*sinλ, -sinU1*cosU2+cosU1*sinU2*cosλ);
+        // 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) < Number.EPSILON ? NaN : Dms.wrap360(α1.toDegrees()),
-            finalBearing:   Math.abs(s) < Number.EPSILON ? NaN : Dms.wrap360(α2.toDegrees()),
+            initialBearing: Math.abs(s) < ε ? NaN : Dms.wrap360(α1.toDegrees()),
+            finalBearing:   Math.abs(s) < ε ? NaN : Dms.wrap360(α2.toDegrees()),
             iterations:     iterations,
@@ -294,0 +304,0 @@ };

mgrs.js

@@ -71,3 +71,4 @@ /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */
     constructor(zone, band, e100k, n100k, easting, northing, datum=LatLonEllipsoidal.datums.WGS84) {
-        if (!(1<=Number(zone) && Number(zone)<=60)) throw new RangeError(`MGRS zone ‘${zone}’ out of range`);
+        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
@@ -290,2 +291,6 @@ if (band.length!=1 || latBands.indexOf(band) == -1) errors.push(`invalid MGRS band ‘${band}’`);
      *
+     * @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.
@@ -298,4 +303,4 @@ * @throws  {Error} If point not valid, if point outside latitude range.
      */
-    toUtm() {
-        const utm = super.toUtm();
+    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);
@@ -302,0 +307,0 @@ }

package.json

@@ -11,4 +11,5 @@ {
   "author": "Chris Veness",
-  "version": "2.1.0",
+  "version": "2.2.0",
   "license": "MIT",
+  "type": "module",
   "engines": {
@@ -28,3 +29,3 @@ "node": ">=8.0.0"
     "coveralls": "^3.0.0",
-    "eslint": "^5.0.0",
+    "eslint": "^6.0.0",
     "esm": "^3.0.0",
@@ -31,0 +32,0 @@ "jsdoc": "^3.0.0",

README.md

@@ -99,3 +99,3 @@ Geodesy functions
 <script type="module">
-    import LatLon from 'https://cdn.jsdelivr.net/npm/geodesy@2.1.0/latlon-spherical.min.js';
+    import LatLon from 'https://cdn.jsdelivr.net/npm/geodesy@2.2.0/latlon-spherical.min.js';
 
@@ -102,0 +102,0 @@ const p1 = new LatLon(50.06632, -5.71475);

test/latlon-ellipsoidal-vincenty-tests.js

@@ -19,2 +19,5 @@ /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */
 
+    const circEquatorial = 40075016.686; // eslint-disable-line no-unused-vars
+    const circMeridional = 40007862.918;
+
     describe('UK', function() {
@@ -45,6 +48,10 @@ const le = new LatLon(50.06632, -5.71475), jog = new LatLon(58.64402, -3.07009);
     describe('antipodal', function() {
-        test('antipodal distance',                   () => new LatLon(0, 0).distanceTo(new LatLon(0.5, 179.5)).should.equal(19936288.579));
+        test('near-antipodal distance',              () => new LatLon(0, 0).distanceTo(new LatLon(0.5, 179.5)).should.equal(19936288.579));
         test('antipodal convergence failure dist',   () => new LatLon(0, 0).distanceTo(new LatLon(0.5, 179.7)).should.be.NaN);
         test('antipodal convergence failure brng i', () => new LatLon(0, 0).initialBearingTo(new LatLon(0.5, 179.7)).should.be.NaN);
         test('antipodal convergence failure brng f', () => new LatLon(0, 0).finalBearingTo(new LatLon(0.5, 179.7)).should.be.NaN);
+        test('antipodal distance equatorial',        () => new LatLon(0, 0).distanceTo(new LatLon(0, 180)).should.equal(circMeridional/2));
+        test('antipodal brng equatorial',            () => new LatLon(0, 0).initialBearingTo(new LatLon(0, 180)).should.equal(0));
+        test('antipodal distance meridional',        () => new LatLon(90, 0).distanceTo(new LatLon(-90, 0)).should.equal(circMeridional/2));
+        test('antipodal brng meridional',            () => new LatLon(90, 0).initialBearingTo(new LatLon(-90, 0)).should.equal(0));
     });
@@ -51,0 +58,0 @@

test/utm-mgrs-tests.js

@@ -43,2 +43,3 @@ /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */
         test('zone fail',       () => should.Throw(function() { new Utm(61, 'N', 0, 0); }, RangeError, 'invalid UTM zone ‘61’'));
+        test('zone fail',       () => should.Throw(function() { new Utm(1.5, 'N', 0, 0); }, RangeError, 'invalid UTM zone ‘1.5’'));
         test('hemisphere fail', () => should.Throw(function() { new Utm(1, 'E', 0, 0); }, RangeError, 'invalid UTM hemisphere ‘E’'));
@@ -51,3 +52,4 @@ test('easting fail',    () => should.Throw(function() { new Utm(1, 'N', 1001e3, 0); }, RangeError, 'invalid UTM easting ‘1001000’'));
     describe('MGRS constructor fail', function() {
-        test('bad zone',             () => should.Throw(function() { new Mgrs(0, 'C', 'A', 'A', 0, 0); }, RangeError, 'MGRS zone ‘0’ out of range'));
+        test('bad zone',             () => should.Throw(function() { new Mgrs(0, 'C', 'A', 'A', 0, 0); }, RangeError, 'invalid MGRS zone ‘0’'));
+        test('bad zone',             () => should.Throw(function() { new Mgrs(1.5, 'C', 'A', 'A', 0, 0); }, RangeError, 'invalid MGRS zone ‘1.5’'));
         test('bad band',             () => should.Throw(function() { new Mgrs(1, 'A', 'A', 'A', 0, 0); }, RangeError, 'invalid MGRS band ‘A’'));
@@ -81,2 +83,3 @@ test('bad grid sq easting',  () => should.Throw(function() { new Mgrs(1, 'C', 'I', 'A', 0, 0); }, RangeError, 'invalid MGRS 100km grid square column ‘I’ for zone 1'));
         test('-1,-1',              () => new LatLon(-1, -1).toUtm().toString(5).should.equal('30 S 722561.73648 9889402.02748'));
+        test('1,1 Z31',            () => new LatLon( 1,  1).toUtm(30).toString(5).should.equal('30 N 945396.68398 110801.83255'));
         test('eiffel tower',       () => new LatLon( 48.8583,   2.2945).toUtm().toString(3).should.equal('31 N 448251.898 5411943.794'));
@@ -129,2 +132,11 @@ test('sidney o/h',         () => new LatLon(-33.857,  151.215 ).toUtm().toString(3).should.equal('56 S 334873.199 6252266.092'));
 
+    describe('ED50 conversion', function() {
+        const helmertturm = new Utm(33, 'N', 368381.402, 5805291.614, LatLon.datums.ED50); // epsg.io/23033
+        const llED50 = helmertturm.toLatLon();
+        const llWGS84 = llED50.convertDatum(LatLon.datums.WGS84);
+        // TODO: no llWGS84.toUtm()!
+        test('helmertturm ED50', () => llED50.toString('dms', 3).should.equal('52°22′51.446″N, 013°03′58.741″E')); // earth-info.nga.mil/GandG/coordsys/datums/datumorigins.html
+        test('helmertturm WGS84', () => llWGS84.toString('dms', 3).should.equal('52°22′48.931″N, 013°03′54.824″E'));
+    });
+
     describe('IBM coordconvert', function() {
@@ -131,0 +143,0 @@ // https://www.ibm.com/developerworks/library/j-coordconvert/#listing7 (note UTM/MGRS confusion; UTM is rounded, MGRS is truncated; UPS not included)

utm.js

@@ -10,3 +10,3 @@ /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */
 
-import LatLonEllipsoidal, { Dms } from './latlon-ellipsoidal.js';
+import LatLonEllipsoidal, { Dms } from './latlon-ellipsoidal-datum.js';
 
@@ -38,3 +38,4 @@
     /**
-     * Creates a Utm coordinate object.
+     * Creates a Utm coordinate object comprising zone, hemisphere, easting, northing on a given
+     * datum (normally WGS84).
      *
@@ -49,2 +50,5 @@ * @param  {number}        zone - UTM 6° longitudinal zone (1..60 covering 180°W..180°E).
      * @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.
@@ -56,8 +60,11 @@ *
      */
-    constructor(zone, hemisphere, easting, northing, datum=LatLonEllipsoidal.datums.WGS84, convergence=null, scale=null) {
+    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 (!(0<=easting && easting<=1000e3)) throw new RangeError(`invalid UTM easting ‘${easting}’`);
-        if (hemisphere.toUpperCase()=='N' && !(0<=northing && northing<9328094)) throw new RangeError(`invalid UTM northing ‘${northing}’`);
-        if (hemisphere.toUpperCase()=='S' && !(1118414<northing && northing<=10000e3)) throw new RangeError(`invalid UTM northing ‘${northing}’`);
+        if (verifyEN) { // range-check E/N values
+            if (!(0<=easting && easting<=1000e3)) throw new RangeError(`invalid UTM easting ‘${easting}’`);
+            if (hemisphere.toUpperCase()=='N' && !(0<=northing && northing<9328094)) throw new RangeError(`invalid UTM northing ‘${northing}’`);
+            if (hemisphere.toUpperCase()=='S' && !(1118414<northing && northing<=10000e3)) throw new RangeError(`invalid UTM northing ‘${northing}’`);
+        }
         if (!datum || datum.ellipsoid==undefined) throw new TypeError(`unrecognised datum ‘${datum}’`);
@@ -78,2 +85,5 @@
      *
+     * 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   {Utm} utmCoord - UTM coordinate to be converted to latitude/longitude.
@@ -254,5 +264,8 @@ * @returns {LatLon} Latitude/longitude of supplied grid reference.
      *
-     * Implements Karney’s method, using Krüger series to order n^6, giving results accurate to 5nm
+     * 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.
@@ -265,3 +278,3 @@ * @throws  {TypeError} Latitude outside UTM limits.
      */
-    toUtm() {
+    toUtm(zoneOverride=undefined) {
         if (!(-80<=this.lat && this.lat<=84)) throw new RangeError(`latitude ‘${this.lat}’ outside UTM limits`);
@@ -271,3 +284,3 @@
 
-        let zone = Math.floor((this.lon+180)/6) + 1; // longitudinal zone
+        let zone = zoneOverride || Math.floor((this.lon+180)/6) + 1; // longitudinal zone
         let λ0 = ((zone-1)*6 - 180 + 3).toRadians(); // longitude of central meridian
@@ -367,3 +380,3 @@
 
-        return new Utm(zone, h, x, y, this.datum, convergence, scale);
+        return new Utm(zone, h, x, y, this.datum, convergence, scale, !!zoneOverride);
     }
@@ -370,0 +383,0 @@ }

New alerts

License Policy Violation

License

This package is not allowed per your license policy. Review the package's license to ensure compliance.

Found 1 instance in 1 package

Fixed alerts

License Policy Violation

License

This package is not allowed per your license policy. Review the package's license to ensure compliance.

Found 1 instance in 1 package

Improved metrics

Total package byte prevSize

426892

increased by1.21%

Lines of code

6358

increased by0.74%

Is type module

Yes

No dependency changes