OSTN02 in NTv2 Transformation in ArcGIS Desktop

In its original form OSTN02 consists of a plane 700 km by 1250 km grid of 1 km resolution with an eastings shift and northing shift at each grid node.

The plane grid form of OSTN02 is not easily compatible with some GIS systems.  Some of these systems utilise a transformation data format known as NTv2.  NTv2 format is a gridded format but the grid is realised on a latitude/longitude graticule and the shifts are expressed as differences in latitude and longitude.” OS

Ordnance Survey working with Defence Geographic Centre have produced an NTV2 version of their OSTN02 transformation.

To implement the transformation and make it available in ArcGIS Desktop:

1. Download the data file from OS website: http://www.ordnancesurvey.co.uk/business-and-government/help-and-support/navigation-technology/os-net/ostn02-ntv2-format.html

2. Using Windows Explorer, create a new folder called greatbritain (if it doesn’t already exist) using lowercase letters in (On a 64 bit machine this will be in Program Files (x86) ):

  • 9.x – C:\Program Files\ArcGIS\pedata\ntv2\
  • 10.0 – C:\Program Files\ArcGIS\Desktop10.0\pedata\ntv2\

3. Copy the NTv2 grid file(s) to the greatbritain folder. Note: There is no limit to how many  .gsb files you can store in this folder.

4. Download GEOGTRAN.txt file. Paste that file in some folder on the C drive like C:\PEOBJ (please do not put the GEOGTRAN file to the Program Files Folder).Add an Environment Variable of PEOBJEDITHOME to the System Variables with the value of the folder where the GEOGTRAN file will reside  (e.g. C:\PEOBJ ).

When you download and copy the file into the e.g. C:\PEOBJ – please remove the file extension (txt) leaving the file without extension.

The content of GEOGTRAN file is as follow:

GEOGTRAN, 208339, "OSGB_1936_To_WGS_1984_OSTN02_NTV2",
PE_GCS_OSGB_1936, PE_GCS_WGS_1984, PE_MTH_NTV2, PE_PAR_NAME_DATASET,
"Dataset_greatbritain/OSTN02_NTv2"
GEOGTRAN, 208338, "OSGB_1936_To_ETRS_1989_OSTN02_NTV2",
PE_GCS_OSGB_1936, PE_GCS_ETRS_1989, PE_MTH_NTV2,
PE_PAR_NAME_DATASET, "Dataset_greatbritain/OSTN02_NTv2"

5. Open Arc Map and set the Layers projection to Geographic projection OSGB36

6. Add any data projected in OSGB36 (geographic) or cartesian British National Grid.

7. Add ESRI_STREETMAP_World_2D or any other data projected in WGS1984 – do not apply any trasformation:

As you can see there will be the data shift of approximately 100m.

8. Right click on Layers then go to Properties and click Transformation button and select one of the new transformation methods from the drop-down list – in our case it will be OSGB_1936_To_WGS_1984_OSTN02_NTV2:

9. Apply the selected transformation – as you can see the WGS84 data has been properly re-projected to OSGB36 using OSTN02 in NTv2 Transformation method:

The above method is working for ArcGIS Desktop 9.2 to 10.0 – it has been tested and is working fine. If you have any questions please do not hesitate to ask.

Please remember – if you have data in British National Grid (Easting, Northing), you need to project it first to OSGB36 (world projection) to use the method described above.

Python – Define Workspace for SDE Connection

One of the possiblilities is to provide the full path to the sde connection file. You can check it creating the connection in ArcCatalog then right clicking and selecting Properties.


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 arcpy
>>> arcpy.env.workspace="C:\Users\osedok\AppData\Roaming\ESRI\
Desktop10.0\Arc
Catalog\PABLO to TEST_SDE_SERVER.sde"
>>> dataset_list=arcpy.ListDatasets()
>>> for feature_dataset in dataset_list:
...     print feature_dataset
...
MASTERMAP.ITN_Network
PABLO.Administrative
PABLO.Default_Blank_Dataset
PABLO.Default_Dataset
PABLO.Development
>>>

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):

OS Coordinate Transformer - proper OSTN02 transformation

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.

ArcGISServer 10.1 Transformation

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:

ArcGISServer 10.0.1 Transformation - wrong results

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 libraryDownload

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.