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 ellipsoid 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 ellipsoid parameters)
1 convert (geodetic) latitude/longitude 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′ = t_{x} + (1+s) · x − r_{z} · 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 latitude/longitude
While the conversion from geodetic to cartesian is straightforward, converting cartesian to geodetic is a complicated problem.
There are a number of possible approaches: this uses Bowring’s 1985 method (The Accuracy of Geodetic Latitude and Height Equations – Survey Review Vol 28, 218, Oct 1985), which provides approximately μm precision on the earth ellipsoid, in a concise formulation. If you have more demanding requirements, you could compare e.g. Fukushima (2006), Vermeille (2011), Karney (2012), and others.
e² = (a²−b²) / a² | 1^{st} eccentricity of ellipsoid | |
ε² = (a²−b²) / b² | 2^{nd} eccentricity of ellipsoid | |
ν = a · √(1 − e² · sin²φ) | length of normal terminated by minor axis | |
p = √(x² + y²) | distance from minor axis | |
R = √(p² + z²) | polar radius | |
tanβ = (b·z)/(a·p) · (1 + ε² · b/R) | parametric latitude | (17) |
tanφ = (z + ε²·b·sin³β) / (p − e²·a·cos³β) | geodetic latitude | (18) |
tanλ = y / x | longitude | |
h = p·cosφ + z·sinφ − a²/ν | height above ellipsoid | (7) |
Accuracy: the Ordnance Survey say a single Helmert transformation between WGS84 and OSGB36 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 (‘OSTN02’): this is a different level of complexity, if you require such a level of accuracy, I expect you already know more about OSTN02 than me.
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 JavaScript source code, also available on GitHub. Note I use Greek letters in variables representing maths symbols conventionally presented as Greek letters: I value the great benefit in legibility over the minor inconvenience in typing.
With its untyped C-style syntax, JavaScript reads remarkably close to pseudo-code: exposing the algorithms with a minimum of syntactic distractions. 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.
I have extended base JavaScript object prototypes with trim(), toRadians(), and toDegrees() methods: I don’t see great likelihood of conflicts; trim() is now in the language, and degree/radian conversion is ubiquitous.
January 2015: I have refactored the scripts to inter-operate better: as the spherical model and the ellipsoidal
model cannot be used at the same time, they constructor for both is now simply LatLon
; also, the
height
(above ellipsoid surface) parameter has been removed from the constructor.
I offer these scripts for free use and adaptation to balance my debt to the open-source info-verse. You are welcome to re-use these scripts [under an MIT licence, 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-2015 Chris Veness