forked from chrisveness/geodesy
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathosgridref.js
256 lines (208 loc) · 10.5 KB
/
osgridref.js
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* Ordnance Survey Grid Reference functions (c) Chris Veness 2005-2014 */
/* - www.movable-type.co.uk/scripts/latlon-gridref.html */
/* - www.ordnancesurvey.co.uk/docs/support/guide-coordinate-systems-great-britain.pdf */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* jshint node:true *//* global define */
'use strict';
if (typeof module!='undefined' && module.exports) var LatLonE = require('./latlon-ellipsoid.js'); // CommonJS (Node.js)
/**
* Creates an OsGridRef object.
*
* @classdesc Convert OS grid references to/from OSGB latitude/longitude points.
* @requires LatLonE
*
* @constructor
* @param {number} easting - Easting in metres from OS false origin.
* @param {number} northing - Northing in metres from OS false origin.
*
* @example
* var grid = new OsGridRef(651409, 313177);
*/
function OsGridRef(easting, northing) {
// allow instantiation without 'new'
if (!(this instanceof OsGridRef)) return new OsGridRef(easting, northing);
this.easting = Math.floor(Number(easting));
this.northing = Math.floor(Number(northing));
}
/**
* Converts (OSGB36) latitude/longitude to Ordnance Survey grid reference easting/northing coordinate.
*
* @param {LatLonE} point - OSGB36 latitude/longitude.
* @returns {OsGridRef} OS Grid Reference easting/northing.
* @throws {Error} If datum of point is not OSGB36.
*
* @example
* var p = new LatLonE(52.65757, 1.71791, LatLonE.datum.OSGB36);
* var grid = OsGridRef.latLonToOsGrid(p); // grid.toString(): TG 51409 13177
*/
OsGridRef.latLonToOsGrid = function(point) {
if (point.datum != LatLonE.datum.OSGB36) throw new Error('Can only convert OSGB36 point to OsGrid');
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;
return new OsGridRef(E, N);
};
/**
* Converts Ordnance Survey grid reference easting/northing coordinate to (OSGB36) latitude/longitude
*
* @param {OsGridRef} gridref - Easting/northing to be converted to latitude/longitude.
* @returns {LatLonE} Latitude/longitude (in OSGB36) of supplied grid reference.
*
* @example
* var grid = new OsGridRef(651409, 313177);
* var p = OsGridRef.osGridToLatLon(grid); // p.toString(): 52°39′27″N, 001°43′04″E
*/
OsGridRef.osGridToLatLon = function(gridref) {
var E = gridref.easting;
var N = gridref.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*Math.PI/180, λ0 = -2*Math.PI/180; // NatGrid true origin
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;
return new LatLonE(φ.toDegrees(), λ.toDegrees(), LatLonE.datum.OSGB36);
};
/**
* Converts standard grid reference (eg 'SU387148') to fully numeric ref (eg [438700,114800]).
*
* @param {string} gridref - Standard format OS grid reference.
* @returns {OsGridRef} Numeric version of grid reference in metres from false origin, centred on
* supplied grid square.
*
* @example
* var grid = OsGridRef.parse('TG 51409 13177'); // grid: { easting: 651409, northing: 313177 }
*/
OsGridRef.parse = function(gridref) {
gridref = String(gridref).trim();
// 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 e = ((l1-2)%5)*5 + (l2%5);
var n = (19-Math.floor(l1/5)*5) - Math.floor(l2/5);
if (e<0 || e>6 || n<0 || n>12) return new OsGridRef(NaN, NaN);
// skip grid letters to get numeric part of ref, stripping any spaces:
gridref = gridref.slice(2).replace(/ /g,'');
// append numeric part of references to grid index:
e += gridref.slice(0, gridref.length/2);
n += gridref.slice(gridref.length/2);
// normalise to 1m grid, rounding up to centre of grid square:
switch (gridref.length) {
case 0: e += '50000'; n += '50000'; break;
case 2: e += '5000'; n += '5000'; break;
case 4: e += '500'; n += '500'; break;
case 6: e += '50'; n += '50'; break;
case 8: e += '5'; n += '5'; break;
case 10: break; // 10-digit refs are already 1m
default: return new OsGridRef(NaN, NaN);
}
return new OsGridRef(e, n);
};
/**
* Converts ‘this’ numeric grid reference to standard OS grid reference.
*
* @param {number} [digits=6] - Precision of returned grid reference (6 digits = metres).
* @returns {string} This grid reference in standard format.
*/
OsGridRef.prototype.toString = function(digits) {
digits = (typeof digits == 'undefined') ? 10 : digits;
var e = this.easting;
var n = this.northing;
if (isNaN(e) || isNaN(n)) return '??';
// 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;
};
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/** Extend Number object with method to trim whitespace from string
* (q.v. blog.stevenlevithan.com/archives/faster-trim-javascript) */
if (typeof String.prototype.trim == 'undefined') {
String.prototype.trim = function() {
return String(this).replace(/^\s\s*/, '').replace(/\s\s*$/, '');
};
}
/** Extend Number object with method to pad with leading zeros to make it w chars wide
* (q.v. stackoverflow.com/questions/2998784 */
if (typeof Number.prototype.pad == 'undefined') {
Number.prototype.pad = function(w) {
var n = this.toString();
while (n.length < w) n = '0' + n;
return n;
};
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
if (typeof module != 'undefined' && module.exports) module.exports = OsGridRef; // CommonJS
if (typeof define == 'function' && define.amd) define([], function() { return OsGridRef; }); // AMD