/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* UTM / WGS-84 Conversion Functions (c) Chris Veness 2014-2016 */
/* MIT Licence */
/* www.movable-type.co.uk/scripts/latlong-utm-mgrs.html */
/* www.movable-type.co.uk/scripts/geodesy/docs/module-utm.html */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
import LatLon from './latlon-ellipsoidal.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.
*
* 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
* @requires latlon-ellipsoidal
*/
/* Utm - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/**
* UTM coordinates, with functions to parse them and convert them to LatLon points.
*/
class Utm { // note prototype-based class not inheritance-based class
/**
* Creates a Utm coordinate object.
*
* @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] - Meridian convergence (bearing of grid north clockwise from true
* north), in degrees
* @param {number} [scale] - Grid scale factor
* @throws {Error} Invalid UTM coordinate
*
* @example
* import Utm from 'utm';
* var utmCoord = new Utm(31, 'N', 448251, 5411932);
*/
constructor(zone, hemisphere, easting, northing, datum=LatLon.datums.WGS84, convergence=null, scale=null) {
if (!(1<=zone && zone<=60)) throw new Error('Invalid UTM zone '+zone);
if (!hemisphere.match(/[NS]/i)) throw new Error('Invalid UTM hemisphere '+hemisphere);
// range-check easting/northing (with 40km overlap between zones) - is this worthwhile?
//if (!(120e3<=easting && easting<=880e3)) throw new Error('Invalid UTM easting '+ easting);
//if (!(0<=northing && northing<=10000e3)) throw new Error('Invalid UTM northing '+ northing);
this.zone = Number(zone);
this.hemisphere = hemisphere.toUpperCase();
this.easting = Number(easting);
this.northing = Number(northing);
this.datum = datum;
this.convergence = convergence===null ? null : Number(convergence);
this.scale = scale===null ? null : Number(scale);
}
/**
* Converts UTM zone/easting/northing coordinate to latitude/longitude
*
* @param {Utm} utmCoord - UTM coordinate to be converted to latitude/longitude.
* @returns {LatLon} Latitude/longitude of supplied grid reference.
*
* @example
* var grid = new Utm(31, 'N', 448251.795, 5411932.678);
* var latlong = grid.toLatLon(); // 48°51′29.52″N, 002°17′40.20″E
*/
toLatLon() {
var z = this.zone;
var h = this.hemisphere;
var x = this.easting;
var y = this.northing;
if (isNaN(z) || isNaN(x) || isNaN(y)) throw new Error('Invalid coordinate');
var falseEasting = 500e3, falseNorthing = 10000e3;
var a = this.datum.ellipsoid.a, f = this.datum.ellipsoid.f;
// WGS 84: a = 6378137, b = 6356752.314245, f = 1/298.257223563;
var k0 = 0.9996; // UTM scale on the central meridian
x = x - falseEasting; // make x ± relative to central meridian
y = h=='S' ? y - falseNorthing : y; // make y ± relative to equator
// ---- from Karney 2011 Eq 15-22, 36:
var e = Math.sqrt(f*(2-f)); // eccentricity
var n = f / (2 - f); // 3rd flattening
var n2 = n*n, n3 = n*n2, n4 = n*n3, n5 = n*n4, n6 = n*n5;
var A = a/(1+n) * (1 + 1/4*n2 + 1/64*n4 + 1/256*n6); // 2πA is the circumference of a meridian
var η = x / (k0*A);
var ξ = y / (k0*A);
var β = [ null, // note β is one-based array (6th order Krüger expressions)
1/2*n - 2/3*n2 + 37/96*n3 - 1/360*n4 - 81/512*n5 + 96199/604800*n6,
1/48*n2 + 1/15*n3 - 437/1440*n4 + 46/105*n5 - 1118711/3870720*n6,
17/480*n3 - 37/840*n4 - 209/4480*n5 + 5569/90720*n6,
4397/161280*n4 - 11/504*n5 - 830251/7257600*n6,
4583/161280*n5 - 108847/3991680*n6,
20648693/638668800*n6 ];
var ξʹ = ξ;
for (var j=1; j<=6; j++) ξʹ -= β[j] * Math.sin(2*j*ξ) * Math.cosh(2*j*η);
var ηʹ = η;
for (var j=1; j<=6; j++) ηʹ -= β[j] * Math.cos(2*j*ξ) * Math.sinh(2*j*η);
var sinhηʹ = Math.sinh(ηʹ);
var sinξʹ = Math.sin(ξʹ), cosξʹ = Math.cos(ξʹ);
var τʹ = sinξʹ / Math.sqrt(sinhηʹ*sinhηʹ + cosξʹ*cosξʹ);
var τi = τʹ;
do {
var σi = Math.sinh(e*Math.atanh(e*τi/Math.sqrt(1+τi*τi)));
var τiʹ = τi * Math.sqrt(1+σi*σi) - σi * Math.sqrt(1+τi*τi);
var δτi = (τʹ - τiʹ)/Math.sqrt(1+τiʹ*τiʹ)
* (1 + (1-e*e)*τi*τi) / ((1-e*e)*Math.sqrt(1+τi*τi));
τi += δτi;
} while (Math.abs(δτi) > 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
var τ = τi;
var φ = Math.atan(τ);
var λ = Math.atan2(sinhηʹ, cosξʹ);
// ---- convergence: Karney 2011 Eq 26, 27
var p = 1;
for (var j=1; j<=6; j++) p -= 2*j*β[j] * Math.cos(2*j*ξ) * Math.cosh(2*j*η);
var q = 0;
for (var j=1; j<=6; j++) q += 2*j*β[j] * Math.sin(2*j*ξ) * Math.sinh(2*j*η);
var γʹ = Math.atan(Math.tan(ξʹ) * Math.tanh(ηʹ));
var γʺ = Math.atan2(q, p);
var γ = γʹ + γʺ;
// ---- scale: Karney 2011 Eq 28
var sinφ = Math.sin(φ);
var kʹ = Math.sqrt(1 - e*e*sinφ*sinφ) * Math.sqrt(1 + τ*τ) * Math.sqrt(sinhηʹ*sinhηʹ + cosξʹ*cosξʹ);
var kʺ = A / a / Math.sqrt(p*p + q*q);
var k = k0 * kʹ * kʺ;
// ------------
var λ0 = ((z-1)*6 - 180 + 3).toRadians(); // longitude of central meridian
λ += λ0; // move λ from zonal to global coordinates
// round to reasonable precision
var lat = Number(φ.toDegrees().toFixed(11)); // nm precision (1nm = 10^-11°)
var lon = Number(λ.toDegrees().toFixed(11)); // (strictly lat rounding should be φ⋅cosφ!)
var convergence = Number(γ.toDegrees().toFixed(9));
var scale = Number(k.toFixed(12));
var latLong = new LatLon(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}
* @throws Error Invalid UTM coordinate
*
* @example
* var utmCoord = Utm.parse('31 N 448251 5411932');
* // utmCoord: {zone: 31, hemisphere: 'N', easting: 448251, northing: 5411932 }
*/
static parse(utmCoord, datum=LatLon.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');
var zone = utmCoord[0], hemisphere = utmCoord[1], easting = utmCoord[2], northing = utmCoord[3];
return new this(zone, hemisphere, easting, northing, datum);
}
/**
* 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.
*
* @param {number} [digits=0] - Number of digits to appear after the decimal point (3 ≡ mm).
* @returns {string} A string representation of the coordinate.
*/
toString(digits) {
digits = Number(digits||0); // default 0 if not supplied
var z = this.zone;
var h = this.hemisphere;
var e = this.easting;
var n = this.northing;
if (isNaN(z) || !h.match(/[NS]/) || isNaN(e) || isNaN(n)) return '';
return z+' '+h+' '+e.toFixed(digits)+' '+n.toFixed(digits);
}
}
/* LatLon_Utm - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/**
* Extends LatLon with method to convert LatLon points to UTM coordinates.
*
* @extends LatLon
*/
class LatLon_Utm extends LatLon {
/**
* Converts latitude/longitude to UTM coordinate.
*
* Implements Karney’s method, using Krüger series to order n^6, giving results accurate to 5nm
* for distances up to 3900km from the central meridian.
*
* @extends LatLon
* @returns {Utm} UTM coordinate.
* @throws {Error} If point not valid, if point outside latitude range.
*
* @example
* var latlong = new LatLon(48.8582, 2.2945);
* var utmCoord = latlong.toUtm(); // 31 N 448252 5411933
*/
toUtm() {
if (isNaN(this.lat) || isNaN(this.lon)) throw new Error('Invalid point');
if (!(-80<=this.lat && this.lat<=84)) throw new Error('Outside UTM limits');
var falseEasting = 500e3, falseNorthing = 10000e3;
var zone = Math.floor((this.lon+180)/6) + 1; // longitudinal zone
var λ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
var mgrsLatBands = 'CDEFGHJKLMNPQRSTUVWXX'; // X is repeated for 80-84°N
var 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(); }
var φ = this.lat.toRadians(); // latitude ± from equator
var λ = this.lon.toRadians() - λ0; // longitude ± from central meridian
var a = this.datum.ellipsoid.a, f = this.datum.ellipsoid.f;
// WGS 84: a = 6378137, b = 6356752.314245, f = 1/298.257223563;
var k0 = 0.9996; // UTM scale on the central meridian
// ---- easting, northing: Karney 2011 Eq 7-14, 29, 35:
var e = Math.sqrt(f*(2-f)); // eccentricity
var n = f / (2 - f); // 3rd flattening
var n2 = n*n, n3 = n*n2, n4 = n*n3, n5 = n*n4, n6 = n*n5; // TODO: compare Horner-form accuracy?
var cosλ = Math.cos(λ), sinλ = Math.sin(λ), tanλ = Math.tan(λ);
var τ = Math.tan(φ); // τ ≡ tanφ, τʹ ≡ tanφʹ; prime (ʹ) indicates angles on the conformal sphere
var σ = Math.sinh(e*Math.atanh(e*τ/Math.sqrt(1+τ*τ)));
var τʹ = τ*Math.sqrt(1+σ*σ) - σ*Math.sqrt(1+τ*τ);
var ξʹ = Math.atan2(τʹ, cosλ);
var ηʹ = Math.asinh(sinλ / Math.sqrt(τʹ*τʹ + cosλ*cosλ));
var A = a/(1+n) * (1 + 1/4*n2 + 1/64*n4 + 1/256*n6); // 2πA is the circumference of a meridian
var α = [ 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 ];
var ξ = ξʹ;
for (var j=1; j<=6; j++) ξ += α[j] * Math.sin(2*j*ξʹ) * Math.cosh(2*j*ηʹ);
var η = ηʹ;
for (var j=1; j<=6; j++) η += α[j] * Math.cos(2*j*ξʹ) * Math.sinh(2*j*ηʹ);
var x = k0 * A * η;
var y = k0 * A * ξ;
// ---- convergence: Karney 2011 Eq 23, 24
var pʹ = 1;
for (var j=1; j<=6; j++) pʹ += 2*j*α[j] * Math.cos(2*j*ξʹ) * Math.cosh(2*j*ηʹ);
var qʹ = 0;
for (var j=1; j<=6; j++) qʹ += 2*j*α[j] * Math.sin(2*j*ξʹ) * Math.sinh(2*j*ηʹ);
var γʹ = Math.atan(τʹ / Math.sqrt(1+τʹ*τʹ)*tanλ);
var γʺ = Math.atan2(qʹ, pʹ);
var γ = γʹ + γʺ;
// ---- scale: Karney 2011 Eq 25
var sinφ = Math.sin(φ);
var kʹ = Math.sqrt(1 - e*e*sinφ*sinφ) * Math.sqrt(1 + τ*τ) / Math.sqrt(τʹ*τʹ + cosλ*cosλ);
var kʺ = A / a * Math.sqrt(pʹ*pʹ + qʹ*qʹ);
var 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(6)); // nm precision
y = Number(y.toFixed(6)); // nm precision
var convergence = Number(γ.toDegrees().toFixed(9));
var scale = Number(k.toFixed(12));
var h = this.lat>=0 ? 'N' : 'S'; // hemisphere
return new Utm(zone, h, x, y, this.datum, convergence, scale);
}
}
/** Extend Number object with method to pad with leading zeros to make it w chars wide
* (q.v. stackoverflow.com/questions/10073699 */
Number.prototype.pad = function(w, z='0') {
const n = this+'';
return n.length >= w ? n : new Array(w - n.length + 1).join(z) + n;
};
export { Utm as default, LatLon_Utm as LatLon };
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */