w1cdn/spothole
1
0
Fork 0
mirror of https://git.ianrenton.com/ian/spothole.git synced 2026-08-08 19:41:41 +00:00
Code Issues 13 Packages Projects Releases Wiki Activity

Prepare for misc.ianrenton.com going offline

Browse Source
This commit is contained in:
Ian Renton
2026-08-07 18:46:32 +01:00
parent 7dfcdb00e8
commit 602551330e
5 changed files with 1112 additions and 3 deletions
Download Patch File Download Diff File Expand all files Collapse all files
static/vendor/js/leaflet-workedallbritainireland.js
Vendored
+3 -3
View File
@@ -16,7 +16,7 @@ L.WorkedAllBritainIreland = L.LayerGroup.extend({
// Workaround to load the geodesy modules in non-modular code. Once we have loaded all three modules, trigger a
// first draw.
import("https://misc.ianrenton.com/Leaflet.WorkedAllBritainIreland/modules/geodesy/osgridref.js")
import(new URL('./modules/geodesy/osgridref.js', import.meta.url).href)
.then(module => {
this._osGridLibrary = module;
if (this._ieGridLibrary && this._utmLibrary) {
@@ -27,7 +27,7 @@ L.WorkedAllBritainIreland = L.LayerGroup.extend({
console.log("Error loading OS Grid Ref library, GB WAB squares may not be available.");
console.log(error);
});
import("https://misc.ianrenton.com/Leaflet.WorkedAllBritainIreland/modules/geodesy/iegridref.js")
import(new URL('./modules/geodesy/iegridref.js', import.meta.url).href)
.then(module => {
this._ieGridLibrary = module;
if (this._osGridLibrary && this._utmLibrary) {
@@ -38,7 +38,7 @@ L.WorkedAllBritainIreland = L.LayerGroup.extend({
console.log("Error loading IE Grid Ref library, NI WAB squares may not be available.");
console.log(error);
});
import("https://misc.ianrenton.com/Leaflet.WorkedAllBritainIreland/modules/geodesy/utm_ci.js")
import(new URL('./modules/geodesy/utm_ci.js', import.meta.url).href)
.then(module => {
this._utmLibrary = module;
if (this._osGridLibrary && this._ieGridLibrary) {
static/vendor/js/modules/geodesy/LICENSE
Vendored
+22
View File
@@ -0,0 +1,22 @@
The MIT License (MIT)
Copyright (c) 2014 Chris Veness
With some additional code & modifications by Ian Renton, 2025
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
static/vendor/js/modules/geodesy/iegridref.js
Vendored
+326
View File
@@ -0,0 +1,326 @@
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* Ordnance Survey of Ireland Grid Reference funcs (c) Chris Veness 2005-2021 & Ian Renton 2025 */
/* MIT Licence */
/* www.movable-type.co.uk/scripts/latlong-gridref.html */
/* www.movable-type.co.uk/scripts/geodesy-library.html#IeGridRef */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
import LatLonEllipsoidal, { Dms } from 'https://cdn.jsdelivr.net/npm/geodesy@2/latlon-ellipsoidal-datum.js';
/**
* Ordnance Survey of Ireland & Northern Ireland grid reference calculations, based on the
* IeGridRef class in the geodesy library at https://github.com/chrisveness/geodesy
*/
/* IeGridRef - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
const nationalGrid = {
trueOrigin: { lat: 53.5, lon: -8 }, // true origin of Irish grid 53°30′N, 8°W
falseOrigin: { easting: -200e3, northing: -250e3 }, // easting & northing of false origin, metres from true origin
scaleFactor: 1.000035, // scale factor on central meridian
ellipsoid: LatLonEllipsoidal.ellipsoids.Airy1830,
};
/**
* Irish Grid References with methods to parse and convert them to latitude/longitude points.
*/
class IeGridRef {
/**
* Creates an IeGridRef 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 IeGridRef from '/js/geodesy/IeGridRef.js';
* const gridref = new IeGridRef(651409, 313177);
*/
constructor(easting, northing) {
this.easting = Number(easting);
this.northing = Number(northing);
if (isNaN(easting) || this.easting<0 || this.easting>7000e3) throw new RangeError(`invalid easting ‘${easting}’`);
if (isNaN(northing) || this.northing<0 || this.northing>13000e3) throw new RangeError(`invalid northing ‘${northing}’`);
}
/**
* Converts ‘this’ Irish 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 IeGridRef(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
const λ0 = nationalGrid.trueOrigin.lon.toRadians(); // longitude of true origin
const E0 = -nationalGrid.falseOrigin.easting; // easting of true origin
const N0 = -nationalGrid.falseOrigin.northing; // northing of true origin
const F0 = nationalGrid.scaleFactor; // scale factor
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*n2 + (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_IeGridRef(φ.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_IeGridRef which includes toOsGrid()
point = new LatLon_IeGridRef(point.lat, point.lon, point.height, point.datum);
}
return point;
}
/**
* Parses grid reference to IeGridRef object.
*
* Accepts standard grid references (eg 'G 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 {IeGridRef} Numeric version of grid reference in metres from false origin (SW corner of
* supplied grid square).
* @throws {Error} Invalid grid reference.
*
* @example
* const grid = IeGridRef.parse('G 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 IeGridRef(match[1], match[2]);
// validate format
match = gridref.match(/^[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); // 100km square
// shuffle down letters after 'I' since 'I' is not used in grid:
if (l1 > 7) l1--;
// convert grid letters into 100km-square indexes from false origin (grid square SV):
const e100km = l1 % 5;
const n100km = 4 - Math.floor(l1 / 5);
// skip grid letters to get numeric (easting/northing) part of ref
let en = gridref.slice(1).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 IeGridRef(e, n);
}
/**
* Converts ‘this’ numeric grid reference to standard OS of Ireland 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 IeGridRef(651409, 313177).toString(8); // 'TG 5140 1317'
* const gridref = new IeGridRef(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 the numeric equivalent of the grid letters
let l1 = (n100km) * 5 % 25 + e100km % 5;
return null; // haven't done this maths yet
// compensate for skipped 'I' and calculate grid letter
if (l1 > 7) l1++;
const letter = String.fromCharCode(l1 + '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 `${letter} ${e} ${n}`;
}
}
/* LatLon_IeGridRef - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/**
* Extends LatLon class with method to convert LatLon point to Irish Grid Reference.
*
* @extends LatLonEllipsoidal
*/
class LatLon_IeGridRef extends LatLonEllipsoidal {
/**
* Converts latitude/longitude to Ordnance Survey of Ireland grid reference easting/northing coordinate.
*
* @returns {IeGridRef} Irish 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
const λ0 = nationalGrid.trueOrigin.lon.toRadians(); // longitude of true origin
const E0 = -nationalGrid.falseOrigin.easting; // easting of true origin
const N0 = -nationalGrid.falseOrigin.northing; // northing of true origin
const F0 = nationalGrid.scaleFactor; // scale factor
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*n2 + (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 IeGridRef(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_IeGridRef.
*/
convertDatum(toDatum) {
const osieED = super.convertDatum(toDatum); // returns LatLonEllipsoidal_Datum
const osieOSGR = new LatLon_IeGridRef(osieED.lat, osieED.lon, osieED.height, osieED.datum);
return osieOSGR;
}
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
export { IeGridRef as default, LatLon_IeGridRef as LatLon, Dms };
static/vendor/js/modules/geodesy/osgridref.js
Vendored
+348
View File
@@ -0,0 +1,348 @@
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* 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 'https://cdn.jsdelivr.net/npm/geodesy@2/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*n2 + (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(/^[HNOST][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*n2 + (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 };
static/vendor/js/modules/geodesy/utm_ci.js
Vendored
+413
View File
@@ -0,0 +1,413 @@
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* 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 };