Convert co-ordinates between WGS-84 and OSGB36

Before I started writing these geodesy scripts, I assumed a lat/long point was a lat/long point, and that was it. Since then, I have discovered that – if you’re being accurate – you need to know what datum you are working within.

In the UK, a latitude/longitude point in WGS84 (as generally used by GPS systems or any world-wide reference system) can be over 100 metres from the same latitude/longitude point in OSGB36 (as used in all Ordnance Survey mapping). The exact difference varies according to location. If you are using Vincenty formulae for calculating distances between points on an ellipsoidal model of the earth, these differences can become relevant.

Note: this is why the Greenwich Meridian shows up about 100 metres from the 0° meridian on GPS units: try 51°28′39″N, 000°00′00″W as an OSGB36 point...

To convert between datums, a ‘Helmert transformation’ is used. The Ordnance Survey explain the details in section 6 (and the annexes) of their Guide to coordinate systems in Great Britain.

The procedure is:

convert polar lat/long/height φ, λ, H to cartesian co-ordinates x, y, z (using the source ellipse parameters)
apply 7-parameter Helmert transformation; this applies a 3-dimensional shift, rotation, and scale
convert cartesian co-ordinates back to polar (using the destination ellipse parameters)

1 convert polar co-ordinates to cartesian:

e² = (a²−b²) / a² (eccentricity of source ellipsoid )
ν = a / √(1−e².sin²φ) (transverse radius of curvature)
x = (ν+H).cosφ.cosλ
y = (ν+H).cosφ.sinλ
z = ((1−e²).ν+H).sinφ

2 apply Helmert transform

The Helmert transform is given by:

transform matrix

Hence:

x′ = t_x + (1+s).x − r_x.y + r_y.z
y′ = t_y + r_z.x + (1+s).y − r_x.z
z′ = t_z − r_y.x + r_x.y + (1+s).z

3 convert cartesian co-ordinates back to polar

e² = (a²−b²) / a² (eccentricity of destination ellipsoid)
p = √(x² + y²)
φ = atan2(z, p.(1−e²)) (initial value); φ′ = 2.π
while (φ−φ′ > 4m) {
	ν = a / √(1−e².sin²φ)
	φ′ = φ
	φ = atan2(z + e².ν.sinφ, p)
}
λ = atan2(y, x)
H = p / cosφ − ν

A single Helmert transformation is accurate to within about 4–5 metres (for greater accuracy, a ‘rubber-sheet’ style transformation, which takes into account distortions in the ‘terrestrial reference frame’ (TRF), must be used: if you’re reading this page, you don’t want to go there...).

The JavaScript should be fairly straightforward to translate into other languages if required – the use of arrays and objects may need to be tweaked for your application/language.

For other scripts for calculating distances, bearings, etc between latitude/longitude points, see my Lat/Long page. I have also done a script for converting between (OSGB36) lat/long and OS grid references.

Principal JavaScript functions are shown below; use ‘View Source’ to see the full JavaScript implementation. You are welcome to re-use these scripts [without any warranty express or implied] provided you retain my copyright notice and when possible a link to my website (under a LGPL license). If you have any queries or find any problems, please contact me.
© 2006 Chris Veness



// ellipse parameters
var e = { WGS84:    { a: 6378137,     b: 6356752.3142, f: 1/298.257223563 },
          Airy1830: { a: 6377563.396, b: 6356256.910,  f: 1/299.3249646   } };

// helmert transform parameters
var h = { WGS84toOSGB36: { tx: -446.448,  ty:  125.157,   tz: -542.060,   // m
                           rx:   -0.1502, ry:   -0.2470,  rz:   -0.8421,  // sec
                           s:    20.4894 },                               // ppm
          OSGB36toWGS84: { tx:  446.448,  ty: -125.157,   tz:  542.060,
                           rx:    0.1502, ry:    0.2470,  rz:    0.8421,
                           s:   -20.4894 } };

                 
function convertOSGB36toWGS84(p1) {
  var p2 = convert(p1, e.Airy1830, h.OSGB36toWGS84, e.WGS84);
  return p2;
}


function convertWGS84toOSGB36(p1) {
  var p2 = convert(p1, e.WGS84, h.WGS84toOSGB36, e.Airy1830);
  return p2;
}


function convert(p1, e1, t, e2) {
  // -- convert polar to cartesian coordinates (using ellipse 1)

  p1.lat = p1.lat.toRad(); p1.lon = p1.lon.toRad(); 

  var a = e1.a, b = e1.b;

  var sinPhi = Math.sin(p1.lat), cosPhi = Math.cos(p1.lat);
  var sinLambda = Math.sin(p1.lon), cosLambda = Math.cos(p1.lon);
  var H = p1.height;

  var eSq = (a*a - b*b) / (a*a);
  var nu = a / Math.sqrt(1 - eSq*sinPhi*sinPhi);

  var x1 = (nu+H) * cosPhi * cosLambda;
  var y1 = (nu+H) * cosPhi * sinLambda;
  var z1 = ((1-eSq)*nu + H) * sinPhi;


  // -- apply helmert transform using appropriate params
  
  var tx = t.tx, ty = t.ty, tz = t.tz;
  var rx = t.rx/3600 * Math.PI/180;  // normalise seconds to radians
  var ry = t.ry/3600 * Math.PI/180;
  var rz = t.rz/3600 * Math.PI/180;
  var s1 = t.s/1e6 + 1;              // normalise ppm to (s+1)

  // apply transform
  var x2 = tx + x1*s1 - y1*rz + z1*ry;
  var y2 = ty + x1*rz + y1*s1 - z1*rx;
  var z2 = tz - x1*ry + y1*rx + z1*s1;


  // -- convert cartesian to polar coordinates (using ellipse 2)

  a = e2.a, b = e2.b;
  var precision = 4 / a;  // results accurate to around 4 metres

  eSq = (a*a - b*b) / (a*a);
  var p = Math.sqrt(x2*x2 + y2*y2);
  var phi = Math.atan2(z2, p*(1-eSq)), phiP = 2*Math.PI;
  while (Math.abs(phi-phiP) > precision) {
    nu = a / Math.sqrt(1 - eSq*Math.sin(phi)*Math.sin(phi));
    phiP = phi;
    phi = Math.atan2(z2 + eSq*nu*Math.sin(phi), p);
  }
  var lambda = Math.atan2(y2, x2);
  H = p/Math.cos(phi) - nu;

  return new LatLon(phi.toDeg(), lambda.toDeg(), H);
}


// ---- the following are duplicated from LatLong.html ---- //


/*
 * construct a LatLon object: arguments in numeric degrees & metres
 *
 * note all LatLong methods expect & return numeric degrees (for lat/long & for bearings)
 */
function LatLon(lat, lon, height) {
  if (arguments.length < 3) height = 0;
  this.lat = lat;
  this.lon = lon;
  this.height = height;
}

/*
 * represent point {lat, lon} in standard representation
 */
LatLon.prototype.toString = function() {
  return this.lat.toLat() + ', ' + this.lon.toLon();
}


// extend String object with method for parsing degrees or lat/long values to numeric degrees
//
// this is very flexible on formats, allowing signed decimal degrees, or deg-min-sec suffixed by 
// compass direction (NSEW). A variety of separators are accepted (eg 3º 37' 09"W) or fixed-width 
// format without separators (eg 0033709W). Seconds and minutes may be omitted. (Minimal validation 
// is done).

String.prototype.parseDeg = function() {
  if (!isNaN(this)) return Number(this);                 // signed decimal degrees without NSEW

  var degLL = this.replace(/^-/,'').replace(/[NSEW]/i,'');  // strip off any sign or compass dir'n
  var dms = degLL.split(/[^0-9.]+/);                     // split out separate d/m/s
  for (var i in dms) if (dms[i]=='') dms.splice(i,1);    // remove empty elements (see note below)
  switch (dms.length) {                                  // convert to decimal degrees...
    case 3:                                              // interpret 3-part result as d/m/s
      var deg = dms[0]/1 + dms[1]/60 + dms[2]/3600; break;
    case 2:                                              // interpret 2-part result as d/m
      var deg = dms[0]/1 + dms[1]/60; break;
    case 1:                                              // decimal or non-separated dddmmss
      if (/[NS]/i.test(this)) degLL = '0' + degLL;       // - normalise N/S to 3-digit degrees
      var deg = dms[0].slice(0,3)/1 + dms[0].slice(3,5)/60 + dms[0].slice(5)/3600; break;
    default: return NaN;
  }
  if (/^-/.test(this) || /[WS]/i.test(this)) deg = -deg; // take '-', west and south as -ve
  return deg;
}
// note: whitespace at start/end will split() into empty elements (except in IE)


// extend Number object with methods for converting degrees/radians

Number.prototype.toRad = function() {  // convert degrees to radians
  return this * Math.PI / 180;
}

Number.prototype.toDeg = function() {  // convert radians to degrees (signed)
  return this * 180 / Math.PI;
}

// extend Number object with methods for presenting bearings & lat/longs

Number.prototype.toDMS = function(dp) {  // convert numeric degrees to deg/min/sec
  if (arguments.length < 1) dp = 0;      // if no decimal places argument, round to int seconds
  var d = Math.abs(this);  // (unsigned result ready for appending compass dir'n)
  var deg = Math.floor(d);
  var min = Math.floor((d-deg)*60);
  var sec = ((d-deg-min/60)*3600).toFixed(dp);
  // fix any nonsensical rounding-up
  if (sec==60) { sec = (0).toFixed(dp); min++; }
  if (min==60) { min = 0; deg++; }
  if (deg==360) deg = 0;
  // add leading zeros if required
  if (deg<100) deg = '0' + deg; if (deg<10) deg = '0' + deg;
  if (min<10) min = '0' + min;
  if (sec<10) sec = '0' + sec;
  return deg + '\u00B0' + min + '\u2032' + sec + '\u2033';
}

Number.prototype.toLat = function(dp) {  // convert numeric degrees to deg/min/sec latitude
  return this.toDMS(dp).slice(1) + (this<0 ? 'S' : 'N');  // knock off initial '0' for lat!
}

Number.prototype.toLon = function(dp) {  // convert numeric degrees to deg/min/sec longitude
  return this.toDMS(dp) + (this>0 ? 'E' : 'W');
}

Movable Type Scripts

Convert co-ordinates between WGS-84 and OSGB36