George MacKerron: code blog

When converting coordinates between WGS84 (GPS) and OSGB36 (UK National Grid), using OSTN02 can gain us a few metres in accuracy over the basic parametric transformation provided by PostGIS’s st_transform via PROJ.4.

Happily, Ordnance Survey now distribute an NTv2 version of the OSTN02 transformation, courtesy of the Defence Geographic Centre, which can be used by PROJ.4 and, therefore, PostGIS.

To set this up, first download and move the NTv2 file into place.

wget http://www.ordnancesurvey.co.uk/oswebsite/gps/docs/OSTN02_NTv2_DataFiles.zip
unzip OSTN02_NTv2_DataFiles.zip
sudo mv OSTN02_NTv2.gsb /usr/share/proj

Then tell your PostGIS-enabled database about it. I use a new SRID 927700 to distinguish this from the standard OSGB36 SRID, 27700. This entry differs from the standard entry by the replacement of the +towgs84 parameter with a +nadgrids parameter.

INSERT INTO public.spatial_ref_sys (srid, auth_name, auth_srid, srtext, proj4text) 
VALUES (
  927700,
  'EPSG',
  27700,
  'PROJCS["OSGB 1936 / British National Grid",GEOGCS["OSGB 1936",DATUM["OSGB_1936",SPHEROID["Airy 1830",6377563.396,299.3249646,AUTHORITY["EPSG","7001"]],AUTHORITY["EPSG","6277"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4277"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",49],PARAMETER["central_meridian",-2],PARAMETER["scale_factor",0.9996012717],PARAMETER["false_easting",400000],PARAMETER["false_northing",-100000],AUTHORITY["EPSG","27700"],AXIS["Easting",EAST],AXIS["Northing",NORTH]]',
  '+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs +nadgrids=OSTN02_NTv2.gsb'
);

Finally, test a transformation, using one of Ordnance Survey’s test coordinate sets.

mappiness=# SELECT st_asewkt(st_transform('SRID=4326;POINT(-0.077731 54.116851 86.778)', 927700));
-[ RECORD 1 ]----------------------------------------------------------
st_asewkt | SRID=927700;POINT(525745.680454343 470703.162768111 86.778)
 
mappiness=# SELECT st_asewkt(st_transform('SRID=4326;POINT(-0.077731 54.116851 86.778)', 27700));
-[ RECORD 1 ]-------------------------------------------------------------------
st_asewkt | SRID=27700;POINT(525744.782236142 470705.422552454 38.5565940965312)

The NTv2-derived northing differs from that of the built-in transformation by nearly 2m.

Limitations

The NTv2-derived coordinates also differ from the true OSTN02 results (525745.658m E, 470703.211m N) by a few mm, owing to limitations of the OSTN02/NTv2 conversion.

Furthermore, the NTv2 transformation passes the altitude/Z coordinate through unchanged, which makes it less accurate here than the built-in transformation (the true altitude derived from OSGM02 for this point is 41.217m).

Finally, PostGIS will happily transform any coordinates you throw at it to the new SRID. For those that are out of OSTN02 range, the results are likely to be horribly wrong (as OS caution, “The 10km off shore buffer is retained in the NTv2 version of OSTN02 so users must be aware that no transformation will take place beyond this limit. Unlike in a dedicated implementation of the original OSTN02 grid, a GIS system using the NTv2 version of OSTN02 is unlikely to warn the user that coordinates are outside the 10km boundary and will therefore not be transformed”).

An experimental solution

To work around the latter two limitations, I’m experimenting with a simple SQL function that wraps the transformation to OSGB36.

It throws away any misleading Z coordinates, so that untransformed values won’t be used. And it compares the results of the NTv2 transformation with those of a transformation to a new dummy SRID that lacks the +nadgrids parameter. If these match to within a metre, it implies that we’re outside the transformation range, since OSTN02 shifts are all much larger than this. In this case, the function just returns null.

INSERT INTO public.spatial_ref_sys (srid, auth_name, auth_srid, srtext, proj4text) 
VALUES (
  927701,
  'EPSG',
  927701,
  '',
  '+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs'
);
 
CREATE OR REPLACE FUNCTION ostn02_transform(geometry) RETURNS geometry AS $$
  SELECT CASE 
    WHEN st_asgeojson(st_transform($1, 927700), 0) <> st_asgeojson(st_transform($1, 927701), 0)
    THEN st_force_2d(st_transform($1, 927700))
    ELSE NULL 
  END;
$$ LANGUAGE SQL immutable strict;

Testing this gives us:

mappiness=# SELECT st_asewkt(ostn02_transform('SRID=4326;POINT(-0.077731 54.116851 86.778)'));
-[ RECORD 1 ]---------------------------------------------------
st_asewkt | SRID=927700;POINT(525745.680454343 470703.162768111)
 
mappiness=# SELECT st_asewkt(ostn02_transform('SRID=4326;POINT(-2.375291 56.175310 100)'));
-[ RECORD 1 ]
st_asewkt |

The first set of coordinates, which are in range, are transformed, minus altitude. For the second set, which are just out of range, null is returned.

An alternative

Note that if you need the OSTN02 transformation in other contexts, or you also need OSGM02 for deriving heights above sea level, I’ve written a C implementation of OSTN02 and OSGM02.