This web page contains worked examples using javascript code. If your browser
displays a security warning, e.g. as Internet Explorer in Windows XP will do,
you can safely take the action to show the active content. Nothing damaging
should happen. If you are unsure, you can view the html source to verify the
safety of doing this! The javascript is all standard and makes no reference
to anything external to this page.
Although there is a lot of involved maths, this is not a difficult task. There are some problems, however, firstly due to the different ways latitude and longitude can be defined. All mapping systems have to use a model of the shape of planet Earth, which is not a perfect sphere. The distance between the poles is less that the distance across the equator, making a shape called an ellipsoid. But the ellipsoid is not a perfect fit either.
Most accurate maps show only part of the surface of the Earth, usually a single country or part of a country. So the usual way of getting round the problem is to define an ellipsoid whose surface has a good alignment with reality over the area the maps is to cover. The Ordnance Survey maps of the United Kingdom use three such reference ellipsoids, one for the Island of Great Britain (England, Scotland and Wales), another for the Island of Ireland, and finally a third ellipsoid to deal with the Channel Islands. This article primarily describes Great Britain, but the same maths can be applied to most of the other two areas.
The traditional ellipsoid adopted for Great Britain is called the Airy Spheroid. This was defined in 1830 by Britain's Astronomer Royal, Sir George Airy. This has a semi-major axis of 6377563.396 metres, and a semi-minor axis of 6356256.910 metres. In the equations that follow we use a constant called the eccentricity, this is simply the diffence between the squares of the axes, divided by the square of the major axis.
function ecc(major, minor)
{
return (major*major - minor*minor) / (major*major);
}
Using the formula above with the Airy axes we get an eccentricity of
Let us assume that that our measurements of latitude and longitude are based on the Airy spheroid. This may not be true, but I'll cover that later on. Firstly we need to convert our measurements from degrees to radians. There is accomplished by multiplying our latitude or longitude by PI / 180. I shall call this conversion factor deg2rad. A few other constants are also defined:
function LLtoNE(lat, lon)
{
var deg2rad = Math.PI / 180;
var rad2deg = 180.0 / Math.PI;
var phi = lat * deg2rad; // convert latitude to radians
var lam = lon * deg2rad; // convert longitude to radians
a = 6377563.396; // OSGB semi-major axis
b = 6356256.91; // OSGB semi-minor axis
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;
// easting
var slat2 = Math.sin(phi) * Math.sin(phi);
var nu = af0 / (Math.sqrt(1 - (e2 * (slat2))));
var rho = (nu * (1 - e2)) / (1 - (e2 * slat2));
var eta2 = (nu / rho) - 1;
var p = lam - lam0;
var IV = nu * Math.cos(phi);
var clat3 = Math.pow(Math.cos(phi),3);
var tlat2 = Math.tan(phi) * Math.tan(phi);
var V = (nu / 6) * clat3 * ((nu / rho) - tlat2);
var clat5 = Math.pow(Math.cos(phi), 5);
var tlat4 = Math.pow(Math.tan(phi), 4);
var VI = (nu / 120) * clat5 * ((5 - (18 * tlat2)) + tlat4 + (14 * eta2) - (58 * tlat2 * eta2));
east = e0 + (p * IV) + (Math.pow(p, 3) * V) + (Math.pow(p, 5) * VI);
// northing
var n = (af0 - bf0) / (af0 + bf0);
var M = Marc(bf0, n, phi0, phi);
var I = M + (n0);
var II = (nu / 2) * Math.sin(phi) * Math.cos(phi);
var III = ((nu / 24) * Math.sin(phi) * Math.pow(Math.cos(phi), 3)) * (5 - Math.pow(Math.tan(phi), 2) + (9 * eta2));
var IIIA = ((nu / 720) * Math.sin(phi) * clat5) * (61 - (58 * tlat2) + tlat4);
north = I + ((p * p) * II) + (Math.pow(p, 4) * III) + (Math.pow(p, 6) * IIIA);
east = Math.round(east); // round to whole number
north = Math.round(north); // round to whole number
nstr = String(north); // convert to string
estr = String(east); // ditto
}
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);
}
This is quite involved maths, and I am not going to try to explain it! This is the
same as that used in the Ordnance Survey document
A Guide to Coordinate Systems in Great Britain and converts latitude and longitude to a Transverse Mercator projection.
These projections are often used for paper (flat) maps. You can find some explanation in the OS document.
The Javascript maths functions are fairly obvious, I hope! You will have to scroll
far to the right to see some of the code. Some of the expressions are quite long-winded,
so I have split them up a little to avoid too many confusing brackets and Math functions.
We can now try this with some test figures, latitude 52° 39’ 27.2531" north,
1° 43’ 4.5177" east. These figures are taken from a worked example in the
Ordnance Survey document.
lat = 52.65757; // decimal version of latitude
lon = 1.71791; // decimal version of longitude
LLtoNE(lat, lon);
document.write("And the results are:");
document.write("Northings: ",nstr, " metres.");
document.write("Eastings: ", estr, " metres.");
The NGR northings and eastings are shown along the axes of any Ordnance Survey map. Normally the NGR is shown as letters and numbers. The letters replace the 100km digit(s) of the eastings and northings. This could easily be determined from a simple look-up table, which would need to handle the ninety one 100km 'squares'. So here is a weird formula, source unknown, which does the job in less space:
function NE2NGR(east, north)
{
var eX = east / 500000;
var nX = north / 500000;
var tmp = Math.floor(eX)-5.0 * Math.floor(nX)+17.0;
nX = 5 * (nX - Math.floor(nX));
eX = 20 - 5.0 * Math.floor(nX) + Math.floor(5.0 * (eX - Math.floor(eX)));
if (eX > 7.5)
eX = eX + 1;
if (tmp > 7.5)
tmp = tmp + 1;
var eing = String(east);
var ning = String(north);
var lnth = eing.length;
eing = eing.substring(lnth - 5, lnth);
lnth = ning.length;
ning = ning.substring(lnth - 5, lnth);
ngr = String.fromCharCode(tmp + 65) + String.fromCharCode(eX + 65) + " " + eing + " " + ning;
return ngr;
}
So for a worked example using the same latitude and longitude, we run the code:
lat = 52.65757;
lon = 1.71791;
LLtoNE(lat, lon);
document.write("And the results are:");
document.write("Northings: ",nstr, " metres.");
document.write("Eastings: ", estr, " metres.");
ngr = NE2NGR(east, north);
document.write("NGR:", ngr);
This is a one metre grid reference, it can be reduced to the standard 100m reference by omitting the last two figures of each five figure group: TG 514 131.
That is basically how it is done. There are a few short cuts, for example I have not taken into account the fact that the NGR grid has a scale factor that varies as you move laterally across the map. You can find appropriate code to deal with this issue in the OS document linked to earlier. I haven't bothered here, in any case the error involved is quite small.
The OS Irish Grid uses a different reference ellipsoid, the 'Modified Airy'. This has a major axis of 6377340.189 metres, and a minor axis of 6356034.447 metres. The NGR origin is obviously different, and the offsets are also different. So in the constants section we need to alter the code as below:
a = 6377340.189; // OSI semi-major
b = 6356034.447; // OSI semi-minor
e0 = 200000; // OSI easting of false origin
n0 = 250000; // OSI northing of false origin
f0 = 1.000035; // OSI scale factor on central meridian
e2 = 0.00667054015; // OSI eccentricity squared
lam0 = -0.13962634015954636615389526147909; // OSI false east
phi0 = 0.93375114981696632365417456114141; // OSI false north
The rest of that procedure LL2NE is the same. There is one simple change to the procedure NE2NGR. The first letter is omitted, thus:
ngr = String.fromCharCode(eX + 65) + " " + eing + " " + ning;
The Channel Islands are slightly different, and use the International 1924 ellipsoid. The major axis is 6378388 metres, the minor is 6356752.3141 metres. Thus the constants are modified:
a = 6378388.000; // INT24 ED50 semi-major
b = 6356911.946; // INT24 ED50 semi-minor
e0 = 500000; // CI easting of false origin
n0 = 0; // CI northing of false origin
f0 = 0.9996; // INT24 ED50 scale factor on central meridian
e2 = 0.0067226700223333; // INT24 ED50 eccentricity squared
lam0 = -0.0523598775598; // CI false east
phi0 = 0 * deg2rad; // CI false north
The rest of that procedure LL2NE is the same. The letters of the resultant NGR are calculated in a different way. A new procedure to replace NE2NGR is required:
function NE2CINGR(east, north)
{
if (north > 5500000)
{
var ngr = "WA " + estr.substring(1, 6) + " " + nstr.substring(2, 7);
}
else if (north < 5500000)
{
var ngr = "WV " + estr.substring(1, 6) + " " + nstr.substring(2, 7);
}
return ngr;
}
Here is a simple demonstation using northings of 5567890 and eastings of 542345:
north = 5567890;
east = 542345;
ngr = NE2CINGR(east, north);
document.write("NGR is: ", ngr);
As I described earlier, map makers use an ellipsoid whose surface closely matches
the surface of the area to be mapped. In order to define how this ellipsoid and the
real Earth are aligned, they use a concept called a Datum (or Terrestrial Reference Frame, TRF).
This Datum ties actual points on the surface of the Earth to latitude and longitude points
on the reference ellipsoid. Since the Earth and Datum do not match exactly, this tie up
distorts the maths a little, and makes it almost impossible to convert one system to another.
Conversions are worthwhile, as they may improve accuracy by one hundred metres or more, but
they are not exact.
Firstly a table showing some ellipsoids and Datums in use in the UK.
| Ellipsoid | Datum | Description |
| Airy | OSGB36 | Used for the OS British National Grid |
| Modified Airy | Ireland 65 | Used for the OS Irish National Grid |
| International 1924 | ED50 | Used for the Channel Islands |
| GRS80 | WGS84 | A system used to describe the whole planet, closely related to the GPS system |
There are many ways to convert data from one system to another, the most
accurate being the most complex! For this example I shall use a 7 parameter Helmert
transformation.
The process is in three parts: (1) Convert latitude and longitude to Cartesian
coordinates (these also include height data, and have three parameters, X, Y and Z).
(2) Transform to the new system by applying the 7 parameters and using a little maths.
(3) Convert back to latitude and longitude.
For the example we shall transform a GRS80 location to Airy, e.g. a GPS reading to
the Airy spheroid.
The following code converts latitude and longitude to Cartesian coordinates. The input parameters are: WGS84 latitude and longitude, axis is the GRS80/WGS84 major axis in metres, ecc is the eccentricity, and height is the height above the ellipsoid.
v = axis / (Math.sqrt (1 - ecc * (Math.pow (Math.sin(lat), 2))));
x = (v + height) * Math.cos(lat) * Math.cos(lon);
y = (v + height) * Math.cos(lat) * Math.sin(lon);
z = ((1 - ecc) * v + height) * Math.sin(lat);
The transformation requires the 7 parameters: xp, yp and zp correct the coordinate origin, xr, yr and zr correct the orientation of the axes, and sf deals with the changing scale factors.
sf = s * 0.000001; // scale factor is in parts per million
xrot = (xr / 3600) * deg2rad;
yrot = (yr / 3600) * deg2rad;
zrot = (zr / 3600) * deg2rad;
hx = x + (x * sf) - (y * zrot) + (z * yrot) + xp;
hy = (x * zrot) + y + (y * sf) - (z * xrot) + yp;
hz = (-1 * x * yrot) + (y * xrot) + z + (z * sf) + zp;
Then we simply convert back to latitude and longitude, where axis and ecc are those of the Airy spheroid:
var lon = Math.atan(y / x);
var p = Math.sqrt((x * x) + (y * y));
var lat = Math.atan(z / (p * (1 - ecc)));
var v = axis / (Math.sqrt(1 - ecc * (Math.sin(lat) * Math.sin(lat))));
var errvalue = 1.0;
var lat0 = 0;
while (errvalue > 0.001)
{
lat0 = Math.atan((z + ecc * v * Math.sin(lat)) / p);
errvalue = Math.abs(lat0 - lat);
lat=lat0;
}
var height = p / Math.cos(lat) - v;
The above involves some iteration, hence the while loop. In this case, ecc is the eccentricity of the Airy spheroid, and axis is its major axis. For convenience, I have joined these three bits of code together into one function:
function transform(lat, lon, a, e, h, a2, e2, xp, yp, zp, xr, yr, zr, s)
{
// convert to cartesian; lat, lon are radians
sf = s * 0.000001;
h = h * 1;
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);
// transform cartesian
xrot = (xr / 3600) * deg2rad;
yrot = (yr / 3600) * deg2rad;
zrot = (zr / 3600) * deg2rad;
hx = x + (x * sf) - (y * zrot) + (z * yrot) + xp;
hy = (x * zrot) + y + (y * sf) - (z * xrot) + yp;
hz = (-1 * x * yrot) + (y * xrot) + z + (z * sf) + zp;
// 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, height: h }; // object to hold lat and lon
return(geo);
}
Which can be called like this:
WGS84_AXIS = 6378137;
WGS84_ECCENTRIC = 0.00669438037928458;
OSGB_AXIS = 6377563.396;
OSGB_ECCENTRIC = 0.0066705397616;
lat = 51.500; // 51.5 degrees north
lon = -0.500; // 0.5 degrees west (negative for west)
phip = lat * deg2rad;
lambdap = lon * deg2rad;
var geo = transform(phip, lambdap, WGS84_AXIS, WGS84_ECCENTRIC, height, OSGB_AXIS, OSGB_ECCENTRIC,
-446.448, 125.157, -542.06, -0.1502, -0.247, -0.8421, 20.4894);
lat = geo.latitude * rad2deg;
lon = geo.longitude * rad2deg;
And now we run the code:
To reverse the transformation process, i.e. to convert OSGB36 to WGS84 you simply swap round the positions of the OSGB and WGS84 axis and eccentricity, and reverse the sign of the seven transformation parameters, like this:
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);
return(geo);
And now we run the code using the result from the last time:
Note the inaccuracy is pretty small, I've displayed the full Javascript output although in practice the error represents a few millimetres. Note also that the transform function also uses the height output parameter in calculating the results above. The height used in this function is the height above the datum. OSGB36 has no concept of height, a different datum based on the average sea level at Newlyn is used to display height on OS maps.
Parameters for transformation from the Irish and the Channel Island grids to WGS84 are shown below. The constants for the semi-major axis and eccentricity can be found further up this page.
var geo = transform(phip, lambdap, IRISH_AXIS, IRISH_ECCENTRIC, height, WGS84_AXIS, WGS84_ECCENTRIC,
482.53, -130.596, 564.557, -1.042, -0.214, -0.631, -8.15);
var geo = transform(phip, lambdap, INT24_AXIS, INT24_ECCENTRIC, height, WGS84_AXIS, WGS84_ECCENTRIC,
-83.901, -98.127, -118.635, 0, 0, 0, 0);
Well, that's about it. Most of the work on this page is carried out by functions in the head section (and also in the body) on this html page. They are slightly different from the functions shown in yellow for presentational purposes, but the maths is the same. All results from this code are shown in green. Note that most of the code is lifted from working calculators elsewhere on this site, hence the use of 'var' may not be consistent! There may be too many brackets here and there, as this code has been converted from other programming languages two or three times!
If you need any more information, the Ordnance Survey document mentioned
previously will probably help you. Information on the OS Irish Grid
can be found at
www.osi.ie/en/alist/geodetic-publications.aspx.
Information on the Channel Islands is harder to come by. Some of my information has
come from an OS map of Jersey, and the rest from hours spent searching the net.
For any more help, please get in touch. My email address can be
found elsewhere on this site.
Roger Muggleton. October 2005
Updated July 2008, August/September 2008, June 2011, November 2017.