Movable Type Home Page

Movable Type Scripts


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...

Functional demo
  Latitude   Longitude
WGS-84  
OSGB36  

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:

  1. convert polar lat/long/height φ, λ, H to cartesian co-ordinates x, y, z (using the source ellipsoid parameters)
  2. apply 7-parameter Helmert transformation; this applies a 3-dimensional shift, rotation, and scale
  3. convert cartesian co-ordinates back to polar (using the destination ellipsoid 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′ = tx + (1+s).x rz.y + ry.z  
y′ = ty + rz.x + (1+s).y rx.z  
z′ = tzry.x + rx.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 learning anything from this page, you don’t want even to think of going there...).

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.


See below for the source code of the JavaScript implementation. These functions should be simple to translate into other languages if required, though the object literal notation (used for ellipsoid & transform parameters) is particularly concise in JavaScript.

The code below will convert between WGS84 and a number of other datums. To extend the list of conversions, you just need the ellipsoid and Helmert transform parameters.

Creative Commons License I offer these formulæ & scripts for free use and adaptation as my contribution to the open-source info-sphere from which I have received so much. You are welcome to re-use these scripts [under a simple attribution license, without any warranty express or implied] provided solely that you retain my copyright notice and a link to this page.

If you would like to show your appreciation and support continued development of these scripts, I would most gratefully accept donations.

If you need any advice or development work done, I am available for consultancy.

If you have any queries or find any problems, contact me at ku.oc.epyt-elbavom@oeg-stpircs

© 2006-2014 Chris Veness


/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */
/*  Convert lat/lon coordinates between different coordinate systems  (c) Chris Veness 2005-2014  */
/*   - www.movable-type.co.uk/scripts/coordtransform.js                                           */
/*   - www.movable-type.co.uk/scripts/latlong-convert-coords.html                                 */
/*                                                                                                */
/*  Usage: to eg convert WGS84 coordinate to OSGB coordinate:                                     */
/*   - var wgs84 = new LatLon(lat, lon, CoordTransform.datum.WGS84);                              */
/*   - var osgb = wgs84.convertDatum(CoordTransform.datum.OSGB36);                                */
/*                                                                                                */
/*  q.v. Ordnance Survey 'A guide to coordinate systems in Great Britain' Section 6               */
/*   - www.ordnancesurvey.co.uk/docs/support/guide-coordinate-systems-great-britain.pdf           */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */

var CoordTransform = {};   // CoordTransform namespace, representing static class

// Ellipsoid parameters
CoordTransform.ellipsoid = {
    WGS84:        { a: 6378137,     b: 6356752.3142,   f: 1/298.257223563 },
    GRS80:        { a: 6378137,     b: 6356752.314140, f: 1/298.257222101 },
    Airy1830:     { a: 6377563.396, b: 6356256.909,    f: 1/299.3249646   },
    AiryModified: { a: 6377340.189, b: 6356034.448,    f: 1/299.32496     },
    Intl1924:     { a: 6378388.000, b: 6356911.946,    f: 1/297.0         },
    Bessel1841:   { a: 6377397.155, b: 6356078.963,    f: 1/299.152815351 }
};

// Datums; with associated ellipsoid and Helmert transform parameters to convert from WGS84 into
// given datum
CoordTransform.datum = {
    WGS84: {
        ellipsoid: CoordTransform.ellipsoid.WGS84,
        transform: { tx:    0.0,    ty:    0.0,     tz:    0.0,    // m
                     rx:    0.0,    ry:    0.0,     rz:    0.0,    // sec
                      s:    0.0 }                                  // ppm
    },
    OSGB36: { // www.ordnancesurvey.co.uk/docs/support/guide-coordinate-systems-great-britain.pdf
        ellipsoid: CoordTransform.ellipsoid.Airy1830,
        transform: { tx: -446.448,  ty:  125.157,   tz: -542.060,  // m
                     rx:   -0.1502, ry:   -0.2470,  rz:   -0.8421, // sec
                      s:   20.4894 }                               // ppm
    },
    ED50: { // og.decc.gov.uk/en/olgs/cms/pons_and_cop/pons/pon4/pon4.aspx
        ellipsoid: CoordTransform.ellipsoid.Intl1924,
        transform: { tx:   89.5,    ty:   93.8,     tz:  123.1,    // m
                     rx:    0.0,    ry:    0.0,     rz:    0.156,  // sec
                      s:   -1.2 }                                  // ppm
    },
    Irl1975: { // maps.osni.gov.uk/CMS_UserFiles/file/The_irish_grid.pdf
        ellipsoid: CoordTransform.ellipsoid.AiryModified,
        transform: { tx: -482.530,  ty:  130.596,   tz: -564.557,  // m
                     rx:   -1.042,  ry:   -0.214,   rz:   -0.631,  // sec
                      s:   -8.150 }                                // ppm
    },
    TokyoJapan: { // www.geocachingtoolbox.com?page=datumEllipsoidDetails
        ellipsoid: CoordTransform.ellipsoid.Bessel1841,
        transform: { tx:  148,      ty: -507,       tz: -685,      // m
                     rx:    0,      ry:    0,       rz:    0,      // sec
                      s:    0 }                                    // ppm
    }
};


/**
 * A LatLon (polar) point has latitude & longitude values and height, on a specified datum
 */
function LatLon(lat, lon, datum, height) {
    if (typeof height == 'undefined') height = 0;

    this.lat = lat;       // in degrees
    this.lon = lon;       // in degrees
    this.datum = datum;
    this.height = height; // height above ellipsoid in metres
}

/**
 * A cartesian point has x, y, z values (in metres)
 */
function Point(x, y, z) {
    this.x = x;
    this.y = y;
    this.z = z;
}


/**
 * Converts 'this' lat/lon coordinate to new coordinate system (to/from WGS84)
 */
LatLon.prototype.convertDatum = function(toDatum) {
    var oldLatLon = this;

    if (oldLatLon.datum == CoordTransform.datum.WGS84) {
        // converting from WGS84
        var transform = toDatum.transform;
    }
    if (toDatum == CoordTransform.datum.WGS84) {
        // converting to WGS84; use inverse transform (don't overwrite original!)
        var transform = {};
        for (var param in oldLatLon.datum.transform) {
            transform[param] = -oldLatLon.datum.transform[param];
        }
    }
    if (typeof transform == 'undefined') throw new Error('Can only convert to/from WGS84');

    // convert polar to cartesian
    var cartesian = oldLatLon.toCartesian();

    // apply transform
    cartesian = cartesian.helmertTransform(transform);

    // convert cartesian to polar
    var newLatLon = cartesian.toPolar(toDatum);

    return newLatLon;
}


/**
 * Converts 'this' point from polar (lat/lon) coordinates to cartesian coordinates
 */
LatLon.prototype.toCartesian = function() {
    var φ = this.lat.toRad(), λ = this.lon.toRad(), H = this.height;
    var a = this.datum.ellipsoid.a, b = this.datum.ellipsoid.b;

    var sinφ = Math.sin(φ);
    var cosφ = Math.cos(φ);
    var sinλ = Math.sin(λ);
    var cosλ = Math.cos(λ);

    var eSq = (a*a - b*b) / (a*a);
    var ν = a / Math.sqrt(1 - eSq*sinφ*sinφ);

    var x = (ν+H) * cosφ * cosλ;
    var y = (ν+H) * cosφ * sinλ;
    var z = ((1-eSq)*ν + H) * sinφ;

    var point = new Point(x, y, z);
    return point;
}

/**
 * Converts 'this' point from cartesian coordinates to polar (lat/lon) coordinates on specified datum
 */
Point.prototype.toPolar = function(datum) {
    var x = this.x, y = this.y, z = this.z;

    var a = datum.ellipsoid.a, b = datum.ellipsoid.b;

    var eSq = (a*a - b*b) / (a*a);
    var p = Math.sqrt(x*x + y*y);
    var phi = Math.atan2(z, p*(1-eSq)), phiP = 2*Math.PI;

    var precision = 4 / a;  // results accurate to around 4 metres
    while (Math.abs(phi-phiP) > precision) {
        ν = a / Math.sqrt(1 - eSq*Math.sin(phi)*Math.sin(phi));
        phiP = phi;
        phi = Math.atan2(z + eSq*ν*Math.sin(phi), p);
    }

    var lambda = Math.atan2(y, x);
    var H = p/Math.cos(phi) - ν;

    var point = new LatLon(phi.toDeg(), lambda.toDeg(), datum, H);
    return point;
}

/**
 * Applies Helmert transform to 'this' point using transform parameters t
 */
Point.prototype.helmertTransform = function(t)   {
    var x1 = this.x, y1 = this.y, z1 = this.z;

    var tx = t.tx, ty = t.ty, tz = t.tz;
    var rx = (t.rx/3600).toRad();  // normalise seconds to radians
    var ry = (t.ry/3600).toRad();
    var rz = (t.rz/3600).toRad();
    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;

    var point = new Point(x2, y2, z2);
    return point;
}


/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */

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

/** Converts numeric degrees to radians */
if (typeof Number.prototype.toRad == 'undefined') {
    Number.prototype.toRad = function() {
        return this * Math.PI / 180;
    }
}

/** Converts radians to numeric (signed) degrees */
if (typeof Number.prototype.toDeg == 'undefined') {
    Number.prototype.toDeg = function() {
        return this * 180 / Math.PI;
    }
}

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */