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:
Hence:
| x′ = tx + (1+s).x − rz.y + ry.z | |
| y′ = ty + rz.x + (1+s).y − rx.z | |
| z′ = tz − ry.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 ellipse & transform parameters) is particularly concise in JavaScript.
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-2012 Chris Veness
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* Coordinate transformations, lat/Long WGS-84 <=> OSGB36 (c) Chris Veness 2005-2012 */
/* - www.movable-type.co.uk/scripts/coordtransform.js */
/* - www.movable-type.co.uk/scripts/latlon-convert-coords.html */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
var CoordTransform = {}; // CoordTransform namespace, representing static class
// ellipse parameters
CoordTransform.ellipse = {
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.910, 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 }
};
// helmert transform parameters from WGS84 to other datums
CoordTransform.datumTransform = {
toOSGB36: { tx: -446.448, ty: 125.157, tz: -542.060, // m
rx: -0.1502, ry: -0.2470, rz: -0.8421, // sec
s: 20.4894 } } // ppm
/**
* Convert lat/lon point in OSGB36 to WGS84
*
* @param {LatLon} pOSGB36: lat/lon in OSGB36 reference frame
* @return {LatLon} lat/lon point in WGS84 reference frame
*/
CoordTransform.convertOSGB36toWGS84 = function(pOSGB36) {
var eAiry1830 = CoordTransform.ellipse.Airy1830;
var eWGS84 = CoordTransform.ellipse.WGS84;
var txToOSGB36 = CoordTransform.datumTransform.toOSGB36;
var txFromOSGB36 = {}; // negate the 'to' transform to get the 'from'
for (var param in txToOSGB36) txFromOSGB36[param] = -txToOSGB36[param];
var pWGS84 = CoordTransform.convertEllipsoid(pOSGB36, eAiry1830, txFromOSGB36, eWGS84);
return pWGS84;
}
/**
* Convert lat/lon point in WGS84 to OSGB36
*
* @param {LatLon} pWGS84: lat/lon in WGS84 reference frame
* @return {LatLon} lat/lon point in OSGB36 reference frame
*/
CoordTransform.convertWGS84toOSGB36 = function(pWGS84) {
var eWGS84 = CoordTransform.ellipse.WGS84;
var eAiry1830 = CoordTransform.ellipse.Airy1830;
var txToOSGB36 = CoordTransform.datumTransform.toOSGB36;
var pOSGB36 = CoordTransform.convertEllipsoid(pWGS84, eWGS84, txToOSGB36, eAiry1830);
return pOSGB36;
}
/**
* Convert lat/lon from one ellipsoidal model to another
*
* q.v. Ordnance Survey 'A guide to coordinate systems in Great Britain' Section 6
* www.ordnancesurvey.co.uk/oswebsite/gps/docs/A_Guide_to_Coordinate_Systems_in_Great_Britain.pdf
*
* @private
* @param {LatLon} point: lat/lon in source reference frame
* @param {Number[]} e1: source ellipse parameters
* @param {Number[]} t: Helmert transform parameters
* @param {Number[]} e1: target ellipse parameters
* @return {Coord} lat/lon in target reference frame
*/
CoordTransform.convertEllipsoid = function(point, e1, t, e2) {
// -- 1: convert polar to cartesian coordinates (using ellipse 1)
var lat = point.lat().toRad();
var lon = point.lon().toRad();
var a = e1.a, b = e1.b;
var sinPhi = Math.sin(lat);
var cosPhi = Math.cos(lat);
var sinLambda = Math.sin(lon);
var cosLambda = Math.cos(lon);
var H = 24.7; // for the moment
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;
// -- 2: apply helmert transform using appropriate params
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;
// -- 3: 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);
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */