/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* Ordnance Survey Grid Reference functions (c) Chris Veness 2005-2016 */
/* MIT Licence */
/* www.movable-type.co.uk/scripts/latlong-gridref.html */
/* www.movable-type.co.uk/scripts/geodesy/docs/module-osgridref.html */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
import LatLon from './latlon-ellipsoidal.js'
/**
* Ordnance Survey grid references provide geocoordinate references for UK mapping purposes.
*
* www.ordnancesurvey.co.uk/docs/support/guide-coordinate-systems-great-britain.pdf.
*
* @module osgridref
* @requires latlon-ellipsoidal
*/
/*
* 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 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/**
* 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 false origin.
* @param {number} northing - Northing in metres from OS false origin.
*
* @example
* import OsGridRef from 'osgridref';
* var gridref = new OsGridRef(651409, 313177);
*/
constructor(easting, northing) {
this.easting = Number(easting);
this.northing = Number(northing);
}
/**
* Converts ‘this’ Ordnance Survey grid reference easting/northing coordinate to latitude/longitude
* (SW corner of grid square).
*
* 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
* var gridref = new OsGridRef(651409.903, 313177.270);
* var pWgs84 = OsGridRef.osGridToLatLon(gridref); // 52°39′28.723″N, 001°42′57.787″E
* // to obtain (historical) OSGB36 latitude/longitude point:
* var pOsgb = OsGridRef.osGridToLatLon(gridref, LatLon.datum.OSGB36); // 52°39′27.253″N, 001°43′04.518″E
*/
toLatLon(datum=LatLon.datums.WGS84) {
var E = this.easting;
var N = this.northing;
var a = 6377563.396, b = 6356256.909; // Airy 1830 major & minor semi-axes
var F0 = 0.9996012717; // NatGrid scale factor on central meridian
var φ0 = (49).toRadians(), λ0 = (-2).toRadians(); // NatGrid true origin is 49°N,2°W
var N0 = -100000, E0 = 400000; // northing & easting of true origin, metres
var e2 = 1 - (b*b)/(a*a); // eccentricity squared
var n = (a-b)/(a+b), n2 = n*n, n3 = n*n*n; // n, n², n³
var φ=φ0, M=0;
do {
φ = (N-N0-M)/(a*F0) + φ;
var Ma = (1 + n + (5/4)*n2 + (5/4)*n3) * (φ-φ0);
var Mb = (3*n + 3*n*n + (21/8)*n3) * Math.sin(φ-φ0) * Math.cos(φ+φ0);
var Mc = ((15/8)*n2 + (15/8)*n3) * Math.sin(2*(φ-φ0)) * Math.cos(2*(φ+φ0));
var Md = (35/24)*n3 * Math.sin(3*(φ-φ0)) * Math.cos(3*(φ+φ0));
M = b * F0 * (Ma - Mb + Mc - Md); // meridional arc
} while (N-N0-M >= 0.00001); // ie until < 0.01mm
var cosφ = Math.cos(φ), sinφ = Math.sin(φ);
var ν = a*F0/Math.sqrt(1-e2*sinφ*sinφ); // nu = transverse radius of curvature
var ρ = a*F0*(1-e2)/Math.pow(1-e2*sinφ*sinφ, 1.5); // rho = meridional radius of curvature
var η2 = ν/ρ-1; // eta = ?
var tanφ = Math.tan(φ);
var tan2φ = tanφ*tanφ, tan4φ = tan2φ*tan2φ, tan6φ = tan4φ*tan2φ;
var secφ = 1/cosφ;
var ν3 = ν*ν*ν, ν5 = ν3*ν*ν, ν7 = ν5*ν*ν;
var VII = tanφ/(2*ρ*ν);
var VIII = tanφ/(24*ρ*ν3)*(5+3*tan2φ+η2-9*tan2φ*η2);
var IX = tanφ/(720*ρ*ν5)*(61+90*tan2φ+45*tan4φ);
var X = secφ/ν;
var XI = secφ/(6*ν3)*(ν/ρ+2*tan2φ);
var XII = secφ/(120*ν5)*(5+28*tan2φ+24*tan4φ);
var XIIA = secφ/(5040*ν7)*(61+662*tan2φ+1320*tan4φ+720*tan6φ);
var 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;
var λ = λ0 + X*dE - XI*dE3 + XII*dE5 - XIIA*dE7;
var point = new LatLon_OsGridRef(φ.toDegrees(), λ.toDegrees(), 0, LatLon.datums.OSGB36);
if (datum != LatLon.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 on Invalid grid reference.
*
* @example
* var 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
var match = gridref.match(/^(\d+),\s*(\d+)$/);
if (match) return new OsGridRef(match[1], match[2]);
// validate format
match = gridref.match(/^[A-Z]{2}\s*[0-9]+\s*[0-9]+$/i);
if (!match) throw new Error('Invalid grid reference');
// get numeric values of letter references, mapping A->0, B->1, C->2, etc:
var l1 = gridref.toUpperCase().charCodeAt(0) - 'A'.charCodeAt(0);
var l2 = gridref.toUpperCase().charCodeAt(1) - 'A'.charCodeAt(0);
// 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):
var e100km = ((l1 - 2) % 5) * 5 + (l2 % 5);
var n100km = (19 - Math.floor(l1 / 5) * 5) - Math.floor(l2 / 5);
// skip grid letters to get numeric (easting/northing) part of ref
var 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 (e100km < 0 || e100km > 6 || n100km < 0 || n100km > 12) throw new Error('Invalid grid reference');
if (en.length != 2) throw new Error('Invalid grid reference');
if (en[0].length != en[1].length) throw new Error('Invalid grid reference');
// standardise to 10-digit refs (metres)
en[0] = (en[0] + '00000').slice(0, 5);
en[1] = (en[1] + '00000').slice(0, 5);
var e = e100km + en[0];
var 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).
* @returns {string} This grid reference in standard format.
*/
toString(digits) {
digits = (digits === undefined) ? 10 : Number(digits);
if (isNaN(digits)) throw new Error('Invalid precision');
var e = this.easting;
var n = this.northing;
if (isNaN(e) || isNaN(n)) throw new Error('Invalid grid reference');
// use digits = 0 to return numeric format (in metres)
if (digits == 0) return e.pad(6) + ',' + n.pad(6);
// get the 100km-grid indices
var e100k = Math.floor(e / 100000), n100k = Math.floor(n / 100000);
if (e100k < 0 || e100k > 6 || n100k < 0 || n100k > 12) return '';
// translate those into numeric equivalents of the grid letters
var l1 = (19 - n100k) - (19 - n100k) % 5 + Math.floor((e100k + 10) / 5);
var l2 = (19 - n100k) * 5 % 25 + e100k % 5;
// compensate for skipped 'I' and build grid letter-pairs
if (l1 > 7) l1++;
if (l2 > 7) l2++;
var letPair = 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));
var gridRef = letPair + ' ' + e.pad(digits / 2) + ' ' + n.pad(digits / 2);
return gridRef;
}
}
/* LatLon_OsGridRef - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/**
* Extends LatLon class with method to convert LatLon point to OS grid reference.
*
* @extends LatLon
*/
class LatLon_OsGridRef extends LatLon {
/**
* Converts latitude/longitude to Ordnance Survey grid reference easting/northing coordinate.
*
* @returns {OsGridRef} OS Grid Reference easting/northing.
*
* @example
* var p = new LatLon(52.65798, 1.71605);
* var grid = OsGridRef.latLonToOsGrid(p); // TG 51409 13177
* // for conversion of (historical) OSGB36 latitude/longitude point:
* var p = new LatLon(52.65757, 1.71791, LatLon.datum.OSGB36);
*/
toOsGrid() {
// if necessary convert to OSGB36 first
var point = this.datum == LatLon.datums.OSGB36 ? this : this.convertDatum(LatLon.datums.OSGB36);
var φ = point.lat.toRadians();
var λ = point.lon.toRadians();
var a = 6377563.396, b = 6356256.909; // Airy 1830 major & minor semi-axes
var F0 = 0.9996012717; // NatGrid scale factor on central meridian
var φ0 = (49).toRadians(), λ0 = (-2).toRadians(); // NatGrid true origin is 49°N,2°W
var N0 = -100000, E0 = 400000; // northing & easting of true origin, metres
var e2 = 1 - (b*b)/(a*a); // eccentricity squared
var n = (a-b)/(a+b), n2 = n*n, n3 = n*n*n; // n, n², n³
var cosφ = Math.cos(φ), sinφ = Math.sin(φ);
var ν = a*F0/Math.sqrt(1-e2*sinφ*sinφ); // nu = transverse radius of curvature
var ρ = a*F0*(1-e2)/Math.pow(1-e2*sinφ*sinφ, 1.5); // rho = meridional radius of curvature
var η2 = ν/ρ-1; // eta = ?
var Ma = (1 + n + (5/4)*n2 + (5/4)*n3) * (φ-φ0);
var Mb = (3*n + 3*n*n + (21/8)*n3) * Math.sin(φ-φ0) * Math.cos(φ+φ0);
var Mc = ((15/8)*n2 + (15/8)*n3) * Math.sin(2*(φ-φ0)) * Math.cos(2*(φ+φ0));
var Md = (35/24)*n3 * Math.sin(3*(φ-φ0)) * Math.cos(3*(φ+φ0));
var M = b * F0 * (Ma - Mb + Mc - Md); // meridional arc
var cos3φ = cosφ*cosφ*cosφ;
var cos5φ = cos3φ*cosφ*cosφ;
var tan2φ = Math.tan(φ)*Math.tan(φ);
var tan4φ = tan2φ*tan2φ;
var I = M + N0;
var II = (ν/2)*sinφ*cosφ;
var III = (ν/24)*sinφ*cos3φ*(5-tan2φ+9*η2);
var IIIA = (ν/720)*sinφ*cos5φ*(61-58*tan2φ+tan4φ);
var IV = ν*cosφ;
var V = (ν/6)*cos3φ*(ν/ρ-tan2φ);
var VI = (ν/120) * cos5φ * (5 - 18*tan2φ + tan4φ + 14*η2 - 58*tan2φ*η2);
var Δλ = λ-λ0;
var Δλ2 = Δλ*Δλ, Δλ3 = Δλ2*Δλ, Δλ4 = Δλ3*Δλ, Δλ5 = Δλ4*Δλ, Δλ6 = Δλ5*Δλ;
var N = I + II*Δλ2 + III*Δλ4 + IIIA*Δλ6;
var E = E0 + IV*Δλ + V*Δλ3 + VI*Δλ5;
N = Number(N.toFixed(3)); // round to mm precision
E = Number(E.toFixed(3));
return new OsGridRef(E, N); // gets truncated to SW corner of 1m grid square
}
}
/** Extend Number object with method to pad with leading zeros to make it w chars wide
* (q.v. stackoverflow.com/questions/2998784 */
Number.prototype.pad = function(w) {
var n = this.toString();
while (n.length < w) n = '0' + n;
return n;
};
export { OsGridRef as default, LatLon_OsGridRef as LatLon };
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */