Following on from the previous page, here is some Javascript to enable conversion from OSGB eastings and northings to latitude and longitude. The procedure is similar to the previous process, with the exception that an iterative process is used to calculate the latitude. In practice, only three or four iterations are required to achieve a good result.
Once again the trigonometry requires all parameters to be in radians, so rad2deg and deg2rad are used for the necessary conversions to and from degees. So here are the functions that do the job (note that the Marc function has an expression spread over four lines, so watch out if you are converting to some other language!):
function NEtoLL(east, north)
{
// converts NGR easting and nothing to lat, lon.
// input metres, output radians
var nX = Number(north);
var eX = Number(east);
a = 6377563.396; // OSGB semi-major
b = 6356256.91; // OSGB semi-minor
e0 = 400000; // OSGB easting of false origin
n0 = -100000; // OSGB northing of false origin
f0 = 0.9996012717; // OSGB scale factor on central meridian
e2 = 0.0066705397616; // OSGB eccentricity squared
lam0 = -0.034906585039886591; // OSGB false east
phi0 = 0.85521133347722145; // OSGB false north
var af0 = a * f0;
var bf0 = b * f0;
var n = (af0 - bf0) / (af0 + bf0);
var Et = east - e0;
var phid = InitialLat(north, n0, af0, phi0, n, bf0);
var nu = af0 / (Math.sqrt(1 - (e2 * (Math.sin(phid) * Math.sin(phid)))));
var rho = (nu * (1 - e2)) / (1 - (e2 * (Math.sin(phid)) * (Math.sin(phid))));
var eta2 = (nu / rho) - 1;
var tlat2 = Math.tan(phid) * Math.tan(phid);
var tlat4 = Math.pow(Math.tan(phid), 4);
var tlat6 = Math.pow(Math.tan(phid), 6);
var clatm1 = Math.pow(Math.cos(phid), -1);
var VII = Math.tan(phid) / (2 * rho * nu);
var VIII = (Math.tan(phid) / (24 * rho * (nu * nu * nu))) * (5 + (3 * tlat2) + eta2 - (9 * eta2 * tlat2));
var IX = ((Math.tan(phid)) / (720 * rho * Math.pow(nu, 5))) * (61 + (90 * tlat2) + (45 * Math.pow(Math.tan(phid), 4) ));
var phip = (phid - ((Et * Et) * VII) + (Math.pow(Et, 4) * VIII) - (Math.pow(Et, 6) * IX));
var X = Math.pow(Math.cos(phid), -1) / nu;
var XI = (clatm1 / (6 * (nu * nu * nu))) * ((nu / rho) + (2 * (tlat2)));
var XII = (clatm1 / (120 * Math.pow(nu, 5))) * (5 + (28 * tlat2) + (24 * tlat4));
var XIIA = clatm1 / (5040 * Math.pow(nu, 7)) * (61 + (662 * tlat2) + (1320 * tlat4) + (720 * tlat6));
var lambdap = (lam0 + (Et * X) - ((Et * Et * Et) * XI) + (Math.pow(Et, 5) * XII) - (Math.pow(Et, 7) * XIIA));
var geo = { latitude: phip, longitude: lambdap };
return(geo);
}
function Marc(bf0, n, phi0, phi)
{
var Marc = bf0 * (((1 + n + ((5 / 4) * (n * n)) + ((5 / 4) * (n * n * n))) * (phi - phi0))
- (((3 * n) + (3 * (n * n)) + ((21 / 8) * (n * n * n))) * (Math.sin(phi - phi0)) * (Math.cos(phi + phi0)))
+ ((((15 / 8) * (n * n)) + ((15 / 8) * (n * n * n))) * (Math.sin(2 * (phi - phi0))) * (Math.cos(2 * (phi + phi0))))
- (((35 / 24) * (n * n * n)) * (Math.sin(3 * (phi - phi0))) * (Math.cos(3 * (phi + phi0)))));
return(Marc);
}
function InitialLat(north, n0, af0, phi0, n, bf0)
{
var phi1 = ((north - n0) / af0) + phi0;
var M = Marc(bf0, n, phi0, phi1);
var phi2 = ((north - n0 - M) / af0) + phi1;
var ind = 0;
while ((Math.abs(north - n0 - M) > 0.00001) && (ind < 20)) // max 20 iterations in case of error
{
ind = ind + 1;
phi2 = ((north - n0 - M) / af0) + phi1;
M = Marc(bf0, n, phi0, phi2);
phi1 = phi2;
}
return(phi2);
}
We can try this with some figures from the last page:
east = 651409;
north = 313177;
geo = NEtoLL(east, north);
lat = geo.latitude * rad2deg;
lon = geo.longitude * rad2deg;
document.write("And the results are:");
document.write("Latitude ", lat);
document.write("Longitude ", lon);
It's easy enough to get this to work for Northern Ireland and the Channel
Islands. All the constants are shown on the previous page. And to transform
from the local grid (OSGB36 in this case) to WGS84 is simply a matter of
reversing the geoid constants and changing the sign of the Helmert parameters:
Instead of
var geo = transform(phip, lambdap, WGS84_AXIS, WGS84_ECCENTRIC, height, OSGB_AXIS, OSGB_ECCENTRIC,
-446.448, 125.157, -542.06, -0.1502, -0.247, -0.8421, 20.4894);
We should call:
var geo = transform(phip, lambdap, OSGB_AXIS, OSGB_ECCENTRIC, height, WGS84_AXIS, WGS84_ECCENTRIC,
446.448, -125.157, 542.06, 0.1502, 0.247, 0.8421, -20.4894);
Finally a simple function to convert an OS NGR into eastings amd northings. Initially it is easiest to pad out the NGR to the 1 metre reference, by padding with zeros. Irish NGRs start with only one letter, not two, so add a leading space first. Make sure letters are upper case. Then you can do the padding with code like this:
if (ngr.length == 4)
t = ngr.substring(0,3) + "0000" + ngr.substring(3) + "0000";
if (ngr.length == 6)
t = ngr.substring(0,4) + "000" + ngr.substring(4) + "000";
if (ngr.length == 8)
t = ngr.substring(0,5) + "00" + ngr.substring(5) + "00";
if (ngr.length == 10)
t = ngr.substring(0,6) + "0" + ngr.substring(6) + "0";
if (ngr.length == 12)
t = ngr;
ngr = t;
This doesn't check for NGR validity, such as too few or too many letters or numbers, extra spaces, odd number of digits, illegal characters, etc. An easy way is to remove all spaces, add a leading space if there are an odd number of characters, then run the code above. The regular expressions in the following code should validate the result. Then convert using a function like the one below. Note that the NGR is tested, using regular expressions, to determine whether it is British, Irish, or the Channel Islands. The latter always begin with WA or WV, and as I stated earlier, the Irish NGR is tested by looking for a space and one letter. Any two letter NGR is assumed to be British, although this who range is not used, the first letter is in the range H to T. You should add another test after the first three 'if' tests, to abort if (grid == ""), i.e. the NGR is not valid.
function conv_ngr_to_ings(ngr)
{
var myRegExp = /^[A-Z][A-Z][0-9]/;
if (myRegExp.test(ngr))
grid = "British";
myRegExp = /^ [A-Z][0-9]/;
if (myRegExp.test(ngr))
grid = "Irish";
myRegExp = /^W[AV][0-9]/;
if (myRegExp.test(ngr))
grid = "Channel Islands";
if (grid == "Channel Islands")
{
var eci = "5" + ngr.substring(2,7);
if (ngr.charAt(1) == "V")
var nci = "54" + ngr.substring(7,12);
else
var nci = "55" + ngr.substring(7,12);
east = Number(eci);
north = Number(nci);
}
if ((grid == "British") || (grid == "Irish"))
{
north = Number(ngr.substring(7));
east = Number(ngr.substring(2,7));
var t1 = ngr.charCodeAt(0) - 65;
if (t1 < 0) // probably a space, so Irish
t1 = 18; // S assumed for Irish (83-65)
if (t1 > 8)
t1 = t1 -1;
var t2 = Math.floor(t1 / 5);
north = north + 500000 * (3 - t2);
east = east + 500000 * (t1 - 5 * t2 - 2);
t1 = ngr.charCodeAt(1) - 65;
if (t1 > 8)
t1 = t1 - 1;
t2 = Math.floor(t1 / 5);
north = north + 100000 * ( 4 - t2);
east = east + 100000 * ( t1 - 5 * t2);
}
}
This can be called like:
ngr = "TQ1234567890";
conv_ngr_to_ings(ngr);
document.write("NGR ", ngr);
document.write("Easting ", east);
document.write("Northing ", north);
I think that just about concludes the examples. I'll leave it to you
to figure out how to convert decimal degrees to degrees, minutes and seconds.
Once again, the main procedure is taken from the OS document.
If there is anything I've missed out, please get in touch.
Roger Muggleton. 24th November 2005
Corrected and updated 5th September 2008
Tracking removed November 2017