mirror of
https://git.ianrenton.com/ian/spothole.git
synced 2026-08-08 19:41:41 +00:00
414 lines
18 KiB
JavaScript
414 lines
18 KiB
JavaScript
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
|
||
/* UTM / WGS-84 Conversion Functions (c) Chris Veness 2014-2022 & Ian Renton 2025 */
|
||
/* 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 'https://cdn.jsdelivr.net/npm/geodesy@2/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<northing && northing<=10000e3)) throw new RangeError(`invalid UTM northing ‘${northing}’`);
|
||
}
|
||
if (!datum || datum.ellipsoid==undefined) throw new TypeError(`unrecognised datum ‘${datum}’`);
|
||
|
||
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.
|
||
*
|
||
* 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.
|
||
* @returns {LatLon} Latitude/longitude of supplied grid reference.
|
||
*
|
||
* @example
|
||
* const grid = new Utm(31, 'N', 448251.795, 5411932.678);
|
||
* const latlong = grid.toLatLon(); // 48°51′29.52″N, 002°17′40.20″E
|
||
*/
|
||
toLatLon() {
|
||
const { zone: z, hemisphere: h } = this;
|
||
|
||
const falseEasting = 500e3, falseNorthing = 10000e3;
|
||
|
||
const { a, f } = this.datum.ellipsoid; // WGS-84: a = 6378137, f = 1/298.257223563;
|
||
|
||
const k0 = 0.9996; // UTM scale on the central meridian
|
||
|
||
const x = this.easting - falseEasting; // make x ± relative to central meridian
|
||
const y = h=='S' ? this.northing - falseNorthing : this.northing; // make y ± relative to equator
|
||
|
||
// ---- from Karney 2011 Eq 15-22, 36:
|
||
|
||
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 A = a/(1+n) * (1 + 1/4*n2 + 1/64*n4 + 1/256*n6); // 2πA is the circumference of a meridian
|
||
|
||
const η = x / (k0*A);
|
||
const ξ = y / (k0*A);
|
||
|
||
const β = [ 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 ];
|
||
|
||
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*η);
|
||
|
||
const sinhηʹ = Math.sinh(ηʹ);
|
||
const sinξʹ = Math.sin(ξʹ), cosξʹ = Math.cos(ξʹ);
|
||
|
||
const τʹ = sinξʹ / Math.sqrt(sinhηʹ*sinhηʹ + cosξʹ*cosξʹ);
|
||
|
||
let δτi = null;
|
||
let τi = τʹ;
|
||
do {
|
||
const σi = Math.sinh(e*Math.atanh(e*τi/Math.sqrt(1+τi*τi)));
|
||
const τiʹ = τi * Math.sqrt(1+σi*σi) - σi * Math.sqrt(1+τi*τi);
|
||
δτ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
|
||
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 a Channel Islands (WA/WV) grid reference.
|
||
*/
|
||
static parseChannelIslandGrid(gridref) {
|
||
// validate format
|
||
let match = gridref.match(/^W[AV]\s*[0-9]+\s*[0-9]+$/i);
|
||
if (!match) throw new Error(`invalid grid reference ‘${gridref}’`);
|
||
|
||
// 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');
|
||
|
||
let utmCoord = "30 N ";
|
||
const e = 5 + en[0];
|
||
utmCoord += e + " ";
|
||
if (gridref.substring(0, 2) === "WA") {
|
||
const n = 55 + en[1];
|
||
utmCoord += n;
|
||
} else if (gridref.substring(0, 2) === "WV") {
|
||
const n = 54 + en[1];
|
||
utmCoord += n;
|
||
}
|
||
|
||
return Utm.parse(utmCoord);
|
||
}
|
||
|
||
|
||
/**
|
||
* 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 };
|