How to convert UK NGR to latitude and longitude

This is a working example of how to convert a British National Grid Reference into Latitude and Longitude, both OSGB36 and WGS84 datums. You can type any NGR (with 2 letters and 2, 4, 6, 8 or 10 numbers, e.g. ST456789) into the green field below, and click on the green button. The NGR is first converted to Easting and Northing, the distance in metres from the National Grid reference point, and these are displayed in the yellow fields.

Enter UK NGR

Eastings
Northings

The Easting and Northing is then converted to Latitude and Longitude, still using the OSGB36 datum based on the Airey ellipsoid which is used for the National Grid. The results, in decimal degrees and in degrees, minutes and seconds, are displayed in the pink fields.

OSGB Latitude
OSGB Longitude
OSGB Latitude
OSGB Longitude

To convert OSGB36 latitude and longitude into WGS84 latitude and longitude several steps take place. Firstly the OSGB is translated to 3-D Cartesian Co-ordinates based on the OSGB36 model. These then undergo a 7-parameter Helmert transformation which rotates, translates and rescales the co-ordinates so as to apply to the WGS84 model. Then The new Cartesian figures are translated to WGS84 latitude and longitude. These figures are shown in blue below.

WGS84 Latitude
WGS84 Longitude
WGS84 Latitude
WGS84 Longitude

Note this live calculator accepts only NGRs from England, Wales, Scotland and their associated islands. The island of Ireland, and the Channel Islands (Jersey etc.), use different ellipsoids and datums to base their National Grids on. It is not difficult to change this example to cover these areas, the required parameters are shown elsewhere on this website.

How does it work?

The javascript code uses a few global variables to hold certain values:


// global conversions
var deg2rad = Math.PI / 180;
var rad2deg = 180.0 / Math.PI;
//globals
var north;    // osgb northings
var east;     // osgb eastings
var osgblon;  // osgb longitude degrees.decimals
var osgblat;  // osgb latitude degrees.decimals
var wgslat;   // wgs84 latitude degrees.decimals
var wgslon;   // wgs84 longitude degrees.decimals
var phip;     // osgb latitude radians
var lambdap;  // osgb longitude radians

Clicking on the green calculate button sends to NGR to the val_ngr() function. This checks the NGR for validity, mainly the right number of letters and numbers, expands it to show the full 1 metre reference which it displays in the green field. It then calls the other functions required.

The first function, called conv_ngr_to_ings(ngr) is shown below:


// function to derive Eastings and Northings from NGR
function conv_ngr_to_ings(ngr)
{
    north = Number(ngr.substring(7));
    east = Number(ngr.substring(2,7));
    var t1 = ngr.charCodeAt(0) - 65;
    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);

    // display results
    document.myform2.northings.value = north;
    document.myform2.eastings.value = east;
}

There's not much to say about this. It effectively converts the grid letters and numbers into two-dimensional distanced from the NGR origin. These are displayed in the yellow fields and stored in the globals east and north.

The next function called, ne2ll(north, east) does the maths to convert Easting and Northing to latitude and longitude:


// function converts NGR easting and nothing to lat, lon.
// input in metres
function ne2ll(north, east)
{
  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 * tlat4 ));
  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));
  lambdap = (lam0 + (Et * X) - ((Et * Et * Et) * XI) + (Math.pow(Et, 5) * XII) - (Math.pow(Et, 7) * XIIA));
  osgblat = phip * rad2deg;
  osgblon = lambdap * rad2deg;
  
  // display results
  if(osgblat < 0)
     document.myform2.osgblatlet.value = "S";
  else
     document.myform2.osgblatlet.value = "N";
 if(osgblon < 0)
     document.myform2.osgblonlet.value = "W";
   else
     document.myform2.osgblonlet.value = "E";
  document.myform2.osgblat.value = Math.abs(osgblat);
  document.myform2.osgblon.value = Math.abs(osgblon);
}  

This starts with a number of OSGB parameters, all commented. The trigonometry that follows is best explained by reading the Ordnance Survey document mentioned elsewhere. The resultant values phip and lambdap hold latitude and longitude in radians. These are converted to degrees, osgblat and osgblon. The final few lines decide latitudes below zero as South, longitudes below zero as east, and display the results in the red fields.

These results are OSGB36, so the next stage is to translate them to WGS84. The figures are first convert to Cartesian Co-ordinates, then transformed to WGS84, and then converted back to WGS84 latitude and longitude.

The first part of function convert_to_wgs(phip, lambdap) starts with a list of parameters for conversion to Cartesian. Then the transform() function is called.


function convert_to_wgs(phip, lambdap)
// constants and parameters for OSGB36 to WGS84 only (not Ireland or Channel Isles)
// inputs are latitude and longitude in radians
{
  var WGS84_AXIS = 6378137;
  var WGS84_ECCENTRIC = 0.00669438037928458;
  var OSGB_AXIS = 6377563.396;
  var OSGB_ECCENTRIC = 0.0066705397616;
  var height = 24.7;  // dummy height above ellipsoid in metres
  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);

The parameters passed to transform() are phip, lambdap, the OSGB semi-major axis and eccentricity, the height above the datum, the WGS84 semi-major axis and eccentricity, then the 7 transformation shifts for the Helmert transformation. Leaving this aside for the moment, the function continues as follows:


 // convert radians to degrees
  wgslat = geo.latitude * rad2deg;
  wgslon = geo.longitude * rad2deg; 
 
  // display result in decimal degrees
  if(wgslat < 0)
     document.myform2.wgslatlet.value = "S";
  else
     document.myform2.wgslatlet.value = "N";
  if(wgslon < 0)
     document.myform2.wgslonlet.value = "W";
  else
     document.myform2.wgslonlet.value = "E";
  document.myform2.wgslat.value = Math.abs(wgslat);
  document.myform2.wgslon.value = Math.abs(wgslon);
}

This simply converts the Latitude and Longitude, now WGS84, from radians to degrees and displays them in decimal form. The transform() function is as follows:


// Helmert transformation
// convert decimal latitude, longitude to cartesian coordinates
// lat, lon are radians
function transform(lat, lon, a, e, h, a2, e2, tx, ty, tz, rx, ry, rz, s)
{
  sf = s * 0.000001;
  v = a / (Math.sqrt(1 - (e * (Math.sin(lat) * Math.sin(lat)))));
  x = (v + h) * Math.cos(lat) * Math.cos(lon);
  y = (v + h) * Math.cos(lat) * Math.sin(lon);
  z = ((1 - e) * v + h) * Math.sin(lat);

  // convert rotations in seconds to radians
  xrot = (rx / 3600) * deg2rad;
  yrot = (ry / 3600) * deg2rad;
  zrot = (rz / 3600) * deg2rad;
  
  hx = x + (x * sf) - (y * zrot) + (z * yrot) + tx;
  hy = (x * zrot) + y + (y * sf) - (z * xrot) + ty;
  hz = (-1 * x * yrot) + (y * xrot) + z + (z * sf) + tz;

  // Convert back to lat, lon
  lon = Math.atan(hy / hx);
  p = Math.sqrt((hx * hx) + (hy * hy));
  lat = Math.atan(hz / (p * (1 - e2)));
  v = a2 / (Math.sqrt(1 - e2 * (Math.sin(lat) * Math.sin(lat))));
  errvalue = 1.0;
  lat0 = 0;
  while (errvalue > 0.001)
  {
    lat0 = Math.atan((hz + e2 * v * Math.sin(lat)) / p);  
    errvalue = Math.abs(lat0 - lat);
    lat = lat0;
  }
  h = p / Math.cos(lat) - v;
  var geo = { latitude: lat, longitude: lon };
  return(geo);
}

Quite a lot goes on here. Firstly the Cartesian x, y and z co-ordinates are calculated from the latitude and longitude. Then the Helmert transformation: tx, ty and tz (in metres) transform the co-ordinates. rx, ry and rz (seconds of arc) rotate the co-ordinates, and s rescales them. Then the new co-ordinates are converted back to WGS84 latitude and longitude. This is not such a simple process as you can see above, and the calculation of latitude is performed by an iterative method. All this is explained in the Ordnance Survey document.

Finally there are a few bits I missed out. Inside ne2ll() is a line:
var phid = InitialLat(north, n0, af0, phi0, n, bf0);

The whole function is below:


// called by ne2ll() function
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);  
}

There's another load of iteration in here, with another function Marc() called within the iteration loop. Marc() has a rather long calculation:


// called by function InitialLat()
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);
}

And finally finally there's a function called todms(lat, lon) which simply converts decimal latitude and longitude to degrees, minutes and seconds for the pink and blue fields towards the top of this page. You can see this if you examine the source of this page.

Roger Muggleton 4th Feb 2014.