Geo::Coordinates::OSTN02 - An implementation of the OSGB's OSTN02 transformation
use Geo::Coordinates::OSTN02 qw/OSGB36_to_ETRS89 ETRS89_to_OSGB36/; ($x, $y, $z) = OSGB36_to_ETRS89($e, $n, $h); ($e, $n, $h) = ETRS89_to_OSGB36($x, $y, $z);
The purpose and use of these modules is described briefly in the companion Geo::Coordinates::OSGB modules. In essence they implement in Perl the Ordnance Survey's OSTN02 transformation which is part of the current definition of the British National Grid. You are strongly recommended to read the official public guides and other introductory documents that are published by the Ordnance Survey of Great Britain with a wealth of other technical information about the OSTN02 transformation.
The following functions can be exported from the
None is exported by default.
Transforms from normal OSGB36 grid references to a pseudo-grid based on the WGS84/ETRS89 geoid model, which can then be translated to lat/long coordinates using
grid_to_ll() with the 'WGS84' parameter.
my $elevation = '100'; # in metres my ($e, $n) = parse_grid_trad('TQ304293'); my ($x, $y, $z) = OSGB36_to_ETRS89($e, $n, $elevation); my ($lat, $lon) = grid_to_ll($x, $y, 'ETRS89'); # or 'WGS84'
Elevation will default to 0 metres if you omit it. The routine should only be called in a list context.
Transforms WGS84/ETRS89 pseudo-grid coordinates into OSGB36 grid references.
my ($lat, $lon, $height) = (51.5, -1, 10); # Somewhere in London my ($x, $y) = ll_to_grid($lat, $lon, 'ETRS89'); # or 'WGS84' my ($e, $n, $elevation) = ETRS89_to_OSGB36($x, $y, $height);
The routine should only be called in a list context.
Since 2002 the British Ordnance Survey have defined the UK national grid not as a projection from their own model of the earth (the Airy 1830 geoid, revised 1936, known as OSGB36), but as a simple table of calculated differences from a projection based on the European standard geoid ETRS89 (which is for Europe a functional equivalent of the international WGS84 geoid). This revision is known as OSGM02 and the transformation is called OSTN02.
The idea is that you project your WGS84 latitude and longitude coordinates into WGS84 (or ETRS89) pseudo-grid references, and then look up the appropriate three dimensional correction in the OS table, and add the results to your grid reference to give you a standard British National Grid reference that you can use with British OS maps. Going back the other way is slightly more complicated, as you have to work backwards, but a simple iteration will do the job. This package implements these lookups and adjustments. You give it a three dimensional grid reference (easting, northing, and altitude, all in metres) and the package returns it corrected from one system to the other.
The problem in the implementation is that the table is huge, and most of it is empty as the OS have chosen to leave the transformation undefined for areas that are off shore. So this package only works properly for grid references that are actually on one or more OS maps. The complete table (including all the 0 lines) contains nearly 1 million lines with six data points and a key. In text form as supplied by the OS that is about 36M bytes of table. By leaving out the lines where the transformation is undefined, omitting a couple of redundant fields, and storing everything as hex strings, this module brings the amount of data down to just under 6M bytes, which loads in about 0.1 sec on my test system. It would be possible to compress the data down to 3M bytes by storing it as packed decimals, but then it would be difficult to include inline in this module, as it would break every time I edited it.
The data is stored below, after the __DATA__ line. Each line is 17 bytes long and represents the transformation parameters for an individual grid square of 1km by 1km. Each line contains five fields all expressed as hexadecimal integers.
Start Length Meaning 0 5 The index of the square 5 4 The x-offset in mm (easting) 9 4 The y-offset in mm (northing) 13 4 The z-offset in mm (elevation)
The index is worked out by multiplying the northing (in km) by 701, adding the easting (in km) and converting to hexadecimal. We add leading zeros so that each key is exactly five bytes long. The maximum value for easting in the OSGB grid is 700km and for northing 1250km, so the maximum possible key value is 876950 = 0xD6196. This fits neatly into 5 bytes.
To keep the other numbers small and positive the data values given for the offsets are actually the amount that they exceed the respective smallest values in the data set. Currently these minima are x: 86.275m, y: -81.603m, and z: 43.982m. So when we read a line from the data we have to convert to decimal, divide by 1000 to get back to metres, and add these minima to the values.
When you load the OSTN02 module, the first thing it does is to load all 309,798 lines into an array called @ostn_data from the internal <DATA> file handle. This is why it takes a little while to load the module, but once loaded it's all very fast. I have also tried using an external file, but this creates issues with distributing the module, and does not offer any increase in speed. Tests with Tie::File improved the load time, but drastically slowed down the rest of the module.
When we need to find the data values for a given grid reference, we work out the appropriate grid square by truncating the easting and northing to the next lowest whole kilometer, and then get the data values for each corner of the square and interpolate between them to find the values for the spot within the square.
The routine that gets these values from the OSTN02 data array is called
_get_ostn_ref. This takes an easting and northing (in km) as arguments and returns a reference to an array of three numbers. These numbers will all be zero if the easting and northing are outside the model.
_get_ostn_ref is the only subroutine that actually touches the OSTN02 data. In order to save space, I've removed the 500,000 zero entries from the model data and so we have to do a binary search instead of simply reading the appropriate line from the array.
This works pretty quickly, the only slow bit is loading the array at the beginning, but it is much faster and needs *much* less memory than loading all the data into the hash. (This would be simpler, and is what the original version did, but it was too slow to be usable and meant that the tests failed on many smaller systems as part of CPAN testing). We do still use a hash, but only to cache lines that we've already found. Read the code for details. This haash-cache only gives a tiny speed up in general, so I might remove it in future versions.
Back in 2008, loading the array took about 1.2 seconds on my Windows machine (a 2.8G Hz Intel Pentium M processor with 2G byte of memory) and just under 0.5 seconds on my Linux system (a 2.8G Hz Intel Pentium 4 processor with 512M bytes of memory). I think that probably said more about speed of the disks I had (and probably the efficiency of Perl under Linux) back then. Today on my Mac Mini the time is down to about 0.1 seconds. Once the data is loaded, calling the routines is reasonably quick.
Please report any to email@example.com
See the test routines included in the distribution.
Copyright (C) 2002-2015 Toby Thurston
OSTN02 transformation data is freely available but remains Crown Copyright (C) 2002
This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
Toby Thurston -- 23 Sep 2015
See Geo::Coordinates::OSGB for the main routines.