Conversion of British National Grid (WKID:27700) to WGS84(WKID:4326) and then to Web Mercator (WKID:3857)

Terms:

1. British National Grid – OSGB36 datum, ellipsoid: Airy 1830 ellipsoid (WKID: 27700, EPSG:27700)
2. Great Britain OSGB36 – OSGB36 datum, ellipsoid: Airy 1830 ellipsoid (WKID:4277, EPSG:4277)

“The Airy 1830 ellipsoid is a bit smaller than Geodetic Reference System 1980 (GRS80) and a slightly different shape. Also, it is not geocentric as GRS80 is: it is designed to lie close to the Geoid beneath the British Isles. Hence, only a tiny fraction of the surface of the ellipsoid has ever been used – the part lying beneath Great Britain. The rest is not useful. So the Airy ellipsoid differs from GRS80 in size, shape, position and orientation, and this is generally true of any pair of geodetic ellipsoids.”- OS

 

3. WGS84 – Datum ellipsoid: GRS80 (WKID: 4326, EPSG:4326) – coordinate system used by GPS system.
4. ETRS89 – European Terrestrial Reference System 1989 – Datum ellipsoid: GRS80 – A high-accuracy version of WGS84. (Zone 28 WKID: 25828, Zone 29 WKID: 25829, Zone 30 WKID: 25830)

 

5. WebMerkator (Auxiliary Sphere) – used by Bing Maps, Google Maps an Open Street Map:

 

    – WKID: 102113 – First used then replaced by 102100 – (deprecated)

    – WKID: 102100, EPSG:3785 – (deprecated)

    – WKID:3857, EPSG:3857 since v.10 of ArcGIS Server – replaced ESRI WKID:102100 with this (current)

 

I decided to write this article as in my work (local authority) the 90% of the spatial data is created using British National Grid coordinate system.

Working with ArcGIS server there was never problem to display the data properly on the top of OS Mastermap or rasters, however few weeks ago I decided to do some mash up using the ArcGIS Java Script API and the basemap dijit and the problems jumped out.
I had to use the British National Grid data over the basemap having the WebMercator projection or a basemap coming from the basemap gallery (You can see it working here).

According to the ArcGIS JavaScript API help, the map’s spatial reference is based on the first layer added to the map and it cannot be overridden (source).
Exactly the same problem you can see when you try to export the features using ArcGIS Server REST API to KML. If you display that data over the WGS84 base map like Google maps you will notice the shift of approximately 100m (datum shift).

See the live example of the problem

The ‘datum shift’ was exactly the origin of the problem. The fact that British National Grid is Transverse Mercator projection (with origin at 49°N, 2°W) based on the Airy 1830 ellipsoid using and WGS84 is based on different GRS80 ellipsoid means the coordinates cannot be transformed directly between those two projections.

An intermediate step is needed for that transformation and before you start converting coordinates I would strongly recommend to read the document below. It’s very simple and clear explanation why we need to transform the coordinates and what is currently possible.

Ordnance Survey Guide – Must read

So, after the above lecture you should have quite good understanding of the problem :).

Let’s have a look at some sample point data.

Sample British National Grid point:

Easting: 420000
Northing: 160000

The correct conversion using OS coordinate transformer service (http://gps.ordnancesurvey.co.uk/convert.asp):

Looks like there will be no problem with this conversion at Version 10.1 of ArcGIS Server the REST Interface of Geometry Server is updated to include Datum Transformation. at 10.1 the OSTN02 NTv2 Transformation is available as EPSG:5339, see below.

However if you need to have a transformation, which one does between WGS 1984 and British National Grid, then SOAP is the only way before the release of Version 10.1. The result from ArcGis 9.3.1 transformation is shown below:

The transformation is wrong, as you can see there is a significant difference in the output X and Y as the transformation does not include the Datum change.

Manual solution:

See the live example

The definitive transformation from ETRS89 that is published by the OSGB is called the National Grid Transformation OSTN02. This models the detailed distortions in the 1936–1962 retriangulation, and achieves backwards compatibility in grid coordinates to sub-metre accuracy.

However for the purposes of displaying not sensitive data data over the web we can use the second method – the Helmert transformation.

For the transformation from  OSGB36 to WGS84 in Britain, using a single Helmert transformation will give errors of up to 5 metres, depending where in the country the points of interest are. One solution to this is to compute a special set of Helmert parameters for a particular region – this is known as a ‘local transformation’.
Alternatively, more complex transformation types that model TRF distortion can be used. An example
of this is the National Grid Transformation OSTN02.

To summarise: For a simple datum change of latitude and longitude coordinates from datum A to datum B, first convert to Cartesian coordinates, taking all ellipsoid heights as zero and using the ellipsoid parameters of datum A; then apply a Helmert transformation from datum A to datum B using equation:

finally convert back to latitude and longitude using the ellipsoid parameters of datum B, discarding the datum B ellipsoid height.

In our case as we have initial coordinates in British National Grid we need to do some more steps:

1. Convert  easting and northing to Lat, Lon based on OSGB36 datum,
2. Convert polar lat/long/height φ, λ, H to cartesian co-ordinates x, y, z (using the source ellipse parameters),
3. apply 7-parameter Helmert transformation; this applies a 3-dimensional shift, rotation, and scale,
4. convert cartesian co-ordinates back to polar (using the destination ellipse parameters – GRS80).

Convert British National Grid to WGS84 – Java Script Solution

To do all these conversion we can use Java Script functions written by Chris Veness – his two documents:

1.Convert co-ordinates between WGS-84 and OSGB36 and
2.Convert between Latitude/Longitude & OS National Grid Reference points

are presenting the methodology and logic used during the coordinate conversion.

[Download slightly modified conversion script  – cords_convert.js]

The further conversion from WGS84 to  WebMerkator (WKID: 102100) can be achieved using the ArcGIS Server, JavaScript API method or using this short Java Script function:

function LatLonWGS84ToWebMercator(wgs84lat, wgs84lon)
{
if((Math.abs(wgs84lon)>180|| Math.abs(wgs84lat)>90)){return;}
var rMajor = 6378137; //Equatorial Radius, WGS84
var num = whs84lon * 0.017453292519943295;
var x = rMajor *num;
var a = wgs84lat * 0.017453292519943295;
var mercatorX = x;
var mercatorY = 3189068.5 * Math.log
((1.0 +Math.sin(a))/(1.0-Math.sin(a)));
return { 'MercatorX': mercatorX, 'MercatorY':mercatorY};
}

Convert British National Grid to WGS84 – Python Script Solution

You need to install the pyproj library – Download

import pyproj

# Define two projections, one for the British National Grid and one for WGS84 (Lat/Lon)
# You can use the full PROJ4 definition or the EPSG identifier (PROJ4 uses a file that matches the two)
#bng = Proj(“+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs towgs84=’446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894′”)
#wgs84 = Proj(‘+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs’)

bng = pyproj.Proj(init='epsg:27700')
wgs84 = pyproj.Proj(init='epsg:4326')

lon,lat = pyproj.transform(bng,wgs84,275331.897,657213.866)

# lon,lat should be now -3.989913896812542, 55.792093458315854
# which is the centre of the roundabout in Motherwell (Merry Street) :
# http://maps.google.com/maps?f=q&source=s_q&hl=en&geocode=&q=55.792093458315854,-3.989913896812542

C:\Python26\ArcGIS10.0>python
Python 2.6.5 (r265:79096, Mar 19 2010, 21:48:26) [MSC v.1500 32 bit (Intel)]on
win32
Type "help", "copyright", "credits" or "license" for more information.
>>> import pyproj
>>> bng=pyproj.Proj(init='epsg:27700')
>>> wgs84 = pyproj.Proj(init='epsg:4326')
>>> lon,lat = pyproj.transform(bng,wgs84,439725,557002)
>>> lat
54.906163255053876
>>> lon
-1.3819797470583637
>>> lon,lat = pyproj.transform(bng,wgs84,275331.897,657213.866)
>>> lon
-3.989913896812542
>>> lat
55.792093458315854
>>>

 

I hope you have found this article useful, please note I am not coordinate system expert and was only adopting some solutions and written this article as there is not much about that so you have few useful links in one place.

Any comments or questions are welcome.