Grid Reference Utilities

Navigation

Skip navigation.

Site search

Site navigation

Details and download

Introduction

In many countries, there are local grid reference systems used by local mapping, allowing easy grid references to be used instead of latitude/longitude. These script libraries provide ways to convert these into their latitude/longitude equivalents. In addition, they provide extra functionality for grid reference systems used in the British Isles, offering ways to convert them into different formats. They can also perform transforms between various different latitude/longitude Earth models, including those used outside the British Isles.

The mathematics used to convert between grid references and latitude/longitude (known as spherical or geodetic coordinates), and between different Earth models is extremely complicated. Minor mistakes can easily become a source of major errors. A great deal of care must be taken to follow the conversion algorithms precisely, and that is what these scripts do for you. Using these scripts allows you to simplify your own scripts, so they can concentrate on using the resulting values.

The scripts also offer methods devoted to geodesics; measuring distances between points, and calculating positions based on geodetic vectors. As well as map projections, there are methods to convert between other types of geographic coordinate systems, such as between spherical and cartesian. There are also methods for geoids; adjusting altitudes based on local gravitational differences.

There are two companion scripts, one for PHP 5.0.3+, and one for JavaScript, allowing you to use them on both the server and client side. The APIs offered by the scripts are almost entirely compatible with each other, receiving compatible inputs and producing compatible outputs. The scripts can be used either separately or together, to provide both client and server-side interaction models.

To download the script(s), and see the script license, use the links on the navigation panel at the top of this page.

Sample conversions

Take this Irish grid reference (pointing at approximately the middle of Lough Neagh in Northern Ireland):
J0259874444
That can also be expressed as (assuming truncation):
302598,374444
Which is equivalent to these Irish Transverse Mercator coordinates:
E 702528m, N 874441m
Which is roughly equivalent to this point on the Irish modified Airy 1830 Earth model:
54.6077240961°, -6.4119822285°
Which can also be expressed as:
54°36'27.806746"N, 6°24'43.136023"W
(or 54°36.46344576'N, 6°24.71893371'W)
Which is roughly equivalent to these GPS coordinates:
54.6078031260°, -6.4129391413°
(also expressed as 54°36'28.091254"N, 6°24'46.580909"W)
Which is equivalent to these UTM coordinates:
29U 667082 6054225
(also expressed as 667082mE, 6054225mN, Zone 29, Northern Hemisphere)
(Which cannot be represented in UPS, since the UPS grid does not extend this far.)
Which is roughly equivalent to this point on the UK Airy 1830 Earth Model:
54.6077391122°, -6.4119726199°
(also expressed as 54°36'27.860804"N, 6°24'43.101432"W)
Which is equivalent to these UK grid coordinates (rounded):
115146,532581
Which can also be expressed as this UK grid reference (truncated):
NW1514632580

The JavaScript version can also be used to compute values without needing to reload the page. Test it here by modifying the contents of any of the boxes below:

Irish grid reference:
or with high precision
Irish grid coordinates:
or with 6 dp precision
ITM coordinates:
Point on the Irish modified Airy 1830 Earth model:
in DD format
in DD format with units
in DM format
in DMS format
GPS coordinates:
in DD format
in DD format with units
in DM format
in DMS format
UTM coordinates:
shorthand
longhand
UPS coordinates:
with latitude bands
with hemisphere
Point on the UK Airy 1830 Earth model:
in DD format
in DD format with units
in DM format
in DMS format
UK grid coordinates:
or with 6 dp precision
UK grid reference:
or with high precision

If your browser supports geolocation, you could even feed the script your (a toolbar or dialog may appear, normally at the top of the page, asking for permission).

My High Accuracy WGS84 GPS To British OSGB36 Grid Reference Conversion tool performs similar conversions (based on Britain only) using the higher accuracy methods provided by these libraries.

Copyright

These scripts were written by the author of the How To Create website; details at the bottom of the page. Script license is given on the sidebar at the top of the page. The algorithms and formulae used are publicly published methods by various geodesy agencies and authors, and links are given throughout the code.

These scripts are normally configured to use polynomial transformation for conversion to Irish Grid, and can optionally use the OSTN15 transformation to British National Grid, and OSGM15 geoid data for both British and Irish height datums. Coefficients and data for these transformations are available under BSD License. If you make use of these specific conversions within your software, then like this page, your software must include the following copyright statement: "© Copyright and database rights Ordnance Survey Limited 2016, © Crown copyright and and database rights Land & Property Services 2016 and/or © Ordnance Survey Ireland, 2016. All rights reserved.".

If you want to avoid having to include that copyright notice (you will still have to follow the How To Create license terms), you can choose to use Helmert transformations instead of OSTN15 and Irish Grid polynomial transformation, and public geoid data instead of OSGM15. The results will not be as accurate.

A note on conversion accuracy

Various stages of the conversions operate on floating point numbers. There are limitations to how computers handle floating point numbers, and some accuracy and precision is lost at various stages of the conversions. The overall effect of this is quite minimal (unless you start feeding the scripts coordinates that lie vast distances outside the expected grid areas), but do not expect the accuracy to be completely perfect. There will also be occasional minor differences between the accuracies of the PHP and JavaScript versions of the script. You can choose to have the scripts return very high precision results, but it is important to note that the precision can easily exceed the accuracy of the conversion methods themselves, and most definitely can exceed the measurement accuracy of most GPS devices.

When converting between British grid references and GPS coordinates, or between different Earth models, errors will be introduced by the conversion algorithms. The most basic conversion is done using Helmert transformations (one of the types of conversion suggested by Ordnance Survey), and these are known to be imperfect, due to the imperfect nature of Earth's own shape. The Helmert transformations between the Ordnance Survey or Ordnance Survey Ireland and GPS Earth models have an error of up to ±5 metres horizontally and vertically. These errors are cumulative, so converting from the OS to OSI model via the GPS model can produce an error of up to ±10 metres. Conversions between other models will depend on the accuracy of their ellipsoids, and the Earth's own shape at that point.

These scripts can also do much more accurate conversions, including the polynomial transformation to the Irish Grid with an accuracy of ±0.4 m for 95% of coordinates, and the OSTN15 transformation to British National Grid. OSTN15 is the method that defines the National Grid itself, when measurements are taken using the OS Net base stations, and is therefore considered to have the perfect accuracy for grid transformation (traditional surveying accuracy is within ±0.1 m, when tested against OSTN15). Polynomial transformation is done by default when converting to Irish Grid.

Note that even when performing high accuracy conversions, the scripts cannot take continental drift into account. Though normally small enough to be below the accuracy for a GPS device, continental drift does change the GPS coordinates for a specific point slowly over time; about 1 m per 40 years in Europe. In 2020, the error is a little under 1 metre. Differential GPS services in the British Isles (and much of Europe) are almost always tied to ETRS89 rather than GPS's normal WGS84. Conversion from ETRS89 coordinates to British National Grid can be done without any substantial error (floating point precision error is the only error worth mentioning) using OSTN15. Conversion between WGS84 and ETRS89 requires external data sources that keep track of the changes over time. While these scripts can help prepare numbers for use with those external services, they cannot determine the continental drift themselves.

Supported grid formats

The scripts are designed to work with map grid formats used by normal maps and mapping software in the British Isles, while having an API to enable conversions to compatible map grids. In addition, they support the most popular global grid reference formats. This means that by default, they support the British National Grid, the Irish Grid, the Irish Transverse Mercator, Universal Transverse Mercator, Universal Polar Stereographic, GPS coordinates, coordinates for the Earth models used by the UK and Irish national grids, and geoids govering the British Isles and the World.

Example uses

Converting between grid references and latitude/longitude GPS coordinates

This is made extremely easy by the scripts, using just two method calls. The first converts the grid reference into numeric coordinates, and the second converts those coordinates into longitude/latitude:

PHP
$grutoolbox = Grid_Ref_Utils::toolbox();
$uk_grid_reference = 'NN 16667 71285';
//convert to a numeric reference
$uk_grid_numbers = $grutoolbox->get_UK_grid_nums($uk_grid_reference);
//convert to global latitude/longitude
$gps_coords = $grutoolbox->grid_to_lat_long($uk_grid_numbers,$grutoolbox->COORDS_GPS_UK,$grutoolbox->HTML);
JavaScript
var grutoolbox = gridRefUtilsToolbox();
var ukGridReference = 'NN 16667 71285';
//convert to a numeric reference
var ukGridNumbers = grutoolbox.getUKGridNums(ukGridReference);
//convert to global latitude/longitude
var gpsCoords = grutoolbox.gridToLatLong(ukGridNumbers,grutoolbox.COORDS_GPS_UK,grutoolbox.HTML);

With UTM or UPS, this is just a single method call. You could optionally also use dd_to_dms/ddToDms to convert it into the desired format.

To perform the reverse conversion, use the lat_long_to_grid/latLongToGrid then get_UK_grid_ref/getUKGridRef. To work with heights, this must be done in two steps, covered below in the section on Helmert transformations.

Converting between UK and Irish grid references

The UK national grid reference system and the Irish national grid reference system do not line up conveniently at all, and although there are significant overlaps between them in various places, there is normally no easy way to translate grid references from one to the other. Within Northern Ireland, for example, grid references are usually expressed using the Irish national grid reference system. Using this library, it is easily possible to obtain the grid reference in the other system:

PHP
$grutoolbox = Grid_Ref_Utils::toolbox();
$irish_grid_reference = 'J0259874444';
//convert to a numeric reference
$irish_grid_numbers = $grutoolbox->get_Irish_grid_nums($irish_grid_reference);
//convert to global latitude/longitude
$gps_coords = $grutoolbox->grid_to_lat_long($irish_grid_numbers,$grutoolbox->COORDS_GPS_IRISH);
//convert to UK numeric reference
$uk_grid_numbers = $grutoolbox->lat_long_to_grid($gps_coords,$grutoolbox->COORDS_GPS_UK);
//convert to UK grid reference
print $grutoolbox->get_UK_grid_ref($uk_grid_numbers,5,$grutoolbox->HTML);
JavaScript
var grutoolbox = gridRefUtilsToolbox();
var irishGridReference = 'J0259874444';
//convert to a numeric reference
var irishGridNumbers = grutoolbox.getIrishGridNums(irishGridReference);
//convert to global latitude/longitude
var gpsCoords = grutoolbox.gridToLatLong(irishGridNumbers,grutoolbox.COORDS_GPS_IRISH);
//convert to UK numeric reference
var ukGridNumbers = grutoolbox.latLongToGrid(gpsCoords,grutoolbox.COORDS_GPS_UK);
//convert to UK grid reference
textNode.nodeValue = grutoolbox.getUKGridRef(ukGridNumbers,5,grutoolbox.TEXT);

A similar approach is used in reverse, or to convert between Irish Transverse Mercator and Irish national grid coordinates. The Irish Transverse Mercator coordinates can either be specified directly as numeric coordinates (so they do not need the extra conversion step), or can be converted into an alternative format using add_grid_units/addGridUnits. Converting between most different grids is done in this way, including UTM and UPS (for the small parts where they overlap), by converting the first reference to GPS coordinates, then converting from GPS coordinates to the other grid reference system.

Converting between Earth models using Helmert transformation

Conversions between different Earth models are only a little harder, as you have to select the appropriate source and destination ellipsoids, and transformation set. You will only need to use this if you are converting from another mapping system which does not use the WGS84/GRS80 Earth model (the scripts automatically apply it for you when converting between the UK/Irish grids and GPS coordinates). It is possible to create your own transformation sets, but to make life easier, the scripts come with transformation sets for the following transforms:

Transformation sets available by default
ConversionFromViaTo
UK to GPSAiry_1830OSGB36_to_WGS84WGS84
GPS to UKWGS84WGS84_to_OSGB36Airy_1830
Irish to GPSAiry_1830_modIreland65_to_WGS84WGS84
GPS to IrishWGS84WGS84_to_Ireland65Airy_1830_mod

Note that GPS software typically uses the WGS84 Earth model, which is a refined version of the GRS80 Earth model. The scripts use the values from WGS84, but you could also choose another ellipsoid, such as the GRS80 Earth model used by Europe's ETRS89. In practice, this makes no significant difference to the resulting values, given the accuracy of the transforms (the practical difference between the two ellipsoids is only about 0.15 mm). Some other ellipsoids produce substantially different coordinates, since they may be different sizes (to better represent the local shape of the Earth), or have different meridians (0 degree longitude).

For example, to convert between the Airy 1830 Earth model used by Ordnance Survey to the modified Airy 1830 Earth model used by Ordnance Survey Ireland, you can use the following code:

PHP
$grutoolbox = Grid_Ref_Utils::toolbox();
$source_coords = Array(54.607720,-6.411990);
//get the ellipsoids that will be used
$Airy_1830 = $grutoolbox->get_ellipsoid('Airy_1830');
$WGS84 = $grutoolbox->get_ellipsoid('WGS84');
$Airy_1830_mod = $grutoolbox->get_ellipsoid('Airy_1830_mod');
//get the transform parameters that will be used
$UK_to_GPS = $grutoolbox->get_transformation('OSGB36_to_WGS84');
$GPS_to_Ireland = $grutoolbox->get_transformation('WGS84_to_Ireland65');
//convert to GPS coordinates
$gps_coords = $grutoolbox->Helmert_transform($source_coords,$Airy_1830,$UK_to_GPS,$WGS84);
//convert to destination coordinates
print $grutoolbox->Helmert_transform($source_coords,$WGS84,$GPS_to_Ireland,$Airy_1830_mod,$grutoolbox->HTML);
JavaScript
var grutoolbox = gridRefUtilsToolbox();
var sourceCoords = [54.607720,-6.411990];
//get the ellipsoids that will be used
var Airy1830 = grutoolbox.getEllipsoid('Airy_1830');
var WGS84 = grutoolbox.getEllipsoid('WGS84');
var Airy1830Mod = grutoolbox.getEllipsoid('Airy_1830_mod');
//get the transform parameters that will be used
var UKToGPS = grutoolbox.getTransformation('OSGB36_to_WGS84');
var GPSToIreland = grutoolbox.getTransformation('WGS84_to_Ireland65');
//convert to GPS coordinates
var gpsCoords = grutoolbox.HelmertTransform(sourceCoords,Airy1830,UKToGPS,WGS84);
//convert to destination coordinates
element.innerHTML = grutoolbox.HelmertTransform(sourceCoords,WGS84,GPSToIreland,Airy1830Mod,grutoolbox.HTML);

It is worth noting that when converting between grid references and latitude/longitude GPS coordinates, heights are ignored completely, since the conversion formulae do not use them. If you wanted to convert heights at the same time, the conversions must be done in two steps, rather than using the COORDS_GPS_UK shortcut (which actually applies a Helmert transformation as well). The conversion between grid coordinates and latitude/longitude must be done using the correct Earth model for the grid coordinates, and a Helmert transformation can then be used to convert between that Earth model and the GPS Earth model, supplying the height during the Helmert transformation step. The conversion below is shown forwards in PHP, and in reverse in JavaScript:

PHP
$grutoolbox = Grid_Ref_Utils::toolbox();
$uk_grid_reference = 'NN 16667 71285';
$uk_height = 1345;
//convert to a numeric reference
$uk_grid_numbers = $grutoolbox->get_UK_grid_nums($uk_grid_reference);
//convert to global latitude/longitude
$uk_coords = $grutoolbox->grid_to_lat_long($uk_grid_numbers,$grutoolbox->COORDS_OS_UK);
//get the ellipsoids that will be used
$Airy_1830 = $grutoolbox->get_ellipsoid('Airy_1830');
$WGS84 = $grutoolbox->get_ellipsoid('WGS84');
//get the transform parameters that will be used
$UK_to_GPS = $grutoolbox->get_transformation('OSGB36_to_WGS84');
//convert to destination coordinates
$gps_coords = $grutoolbox->Helmert_transform($uk_coords,$uk_height,$Airy_1830,$UK_to_GPS,$WGS84);
//convert to degrees-minutes-seconds
print $grutoolbox->dd_to_dms($gps_coords,$grutoolbox->HTML) . ' height ' . $gps_coords[2];
//or it can just output the coordinates and height directly
print $grutoolbox->Helmert_transform($uk_coords,$uk_height,$Airy_1830,$UK_to_GPS,$WGS84,$grutoolbox->HTML);
JavaScript
var grutoolbox = gridRefUtilsToolbox();
var gpsString = '54°36\'27.792000"N, 6°24\'43.164000"W';
var gpsHeight = 1345;
//convert to decimal degrees
var gpsCoords = grutoolbox.dmsToDd(gpsString);
//get the ellipsoids that will be used
var WGS84 = grutoolbox.get_ellipsoid('WGS84');
var Airy1830 = grutoolbox.get_ellipsoid('Airy_1830');
//get the transform parameters that will be used
var GPSToUK = grutoolbox.get_transformation('WGS84_to_OSGB36');
//convert to destination coordinates
var ukCoords = grutoolbox.HelmertTransform(gpsCoords,gpsHeight,WGS84,GPSToUK,Airy1830);
var ukGridNumbers = grutoolbox.latLongToGrid(gpsCoords,grutoolbox.COORDS_OS_UK);
//convert to grid reference
textNode.nodeValue = grutoolbox.getUKGridRef(ukGridNumbers,5,grutoolbox.TEXT) + ' height ' + ukCoords[2];

Converting from other map systems

A large number of local mapping systems have been developed to serve particular countries or areas. In many cases these use the same approach as the mapping systems these scripts already deal with. They begin with a model of the Earth that attempts to reflect its ellipsoidal shape at the target location. They then use a Transverse Mercator projection (or polar stereographic projection) of the ellipsoid in order to produce the map. The map is aligned to one specific part of the Earth model (the true origin), and then usually offset to some other false grid origin, and scaled slightly using a scale factor. This is the datum of that map's projection.

You would need to work out (if needed) the distance from any grid reference on that grid system to its false origin. So if the maps use anything like myriad letters, you need to convert from that format to simple distances from the false origin - easting and northing coordinates. You then need to obtain the ellipsoid parameters used by the Earth model. Feed the ellipsoid parameters into the create_ellipsoid/createEllipsoid method to get an ellipsoid parameter set (if you are lucky enough to find that the map uses one of the Earth models already provided by the scripts, you can just use get_ellipsoid/getEllipsoid method to retrieve the appropriate ellipsoid parameter set). You then need to obtain the datum parameters used by the map projection (mentioned above). Feed all of those, along with the ellipsoid parameter set, into the create_datum/createDatum method to get a datum parameter set.

Then feed your numbers into either grid_to_lat_long/gridToLatLong or polar_to_lat_long/polarToLatLong (depending on what type of projection the map system uses), along with your custom datum parameter set, and it will return the geodetic coordinates of the point on that Earth model. If the mapping system uses the same Earth model as GPS (WGS84 or GRS80), then there is nothing more to do. If it uses a different ellipsoid, you need to also get the Helmert transformation parameters for converting between that Earth model and the WGS84/GRS80 Earth model. Feed those into create_transformation/createTransformation to get a custom Helmert transformation parameter set. Feed your geodetic coordinates, custom Earth model, Helmert transformation parameter set and the script's existing Earth model for WGS84 as the destination Earth model, into Helmert_transform/HelmertTransform, and it will return GPS coordinates for the point.

A similar process works in reverse. That may all sound complicated, but the details of Earth models and datums for most mapping systems are easily available, as are the Helmert transformation parameters, and these scripts take care of most of the rest for you. An example in JavaScript:

var grutoolbox = gridRefUtilsToolbox();
var sourceCoords = [12345,67890];
//get the ellipsoids that will be used
var earthModel = grutoolbox.createEllipsoid(6378206.4,6356583.8);
var datum = grutoolbox.createDatum(earthModel,0.99987,100000,200000,20,119);
var WGS84 = grutoolbox.getEllipsoid('WGS84');
//get the transform parameters that will be used
var modelToGPS = grutoolbox.createTransformation(-8,160,176,0,0,0,0);
//convert from grid coordinates to geodetic coordinates
var geoCoords = grutoolbox.gridToLatLong(sourceCoords,datum);
//convert to GPS coordinates
var gpsCoords = grutoolbox.HelmertTransform(geoCoords,earthModel,modelToGPS,WGS84);

And an example using a polar stereographic projection in PHP:

$grutoolbox = Grid_Ref_Utils::toolbox();
$source_coords = Array(12345,67890);
$earth_model = $grutoolbox->create_ellipsoid(6378206.4,6356583.8);
$datum = $grutoolbox->create_datum($earth_model,0.99987,100000,200000,20,119);
$WGS84 = $grutoolbox->get_ellipsoid('WGS84');
$model_to_gps = $grutoolbox->create_transformation(-8,160,176,0,0,0,0);
//convert from polar grid coordinates to geodetic coordinates
$geo_coords = $grutoolbox->polar_to_lat_long($source_coords,$datum);
$gps_coords = $grutoolbox->Helmert_transform($geo_coords,$earth_model,$model_to_gps,$WGS84);

Of course, it is possible that you need to convert from a mapping system that does not use a Transverse Mercator projection or basic polar stereographic projection. In these cases, you will need to either find another script that supports that mapping system, or write one yourself. I can recommend the EPSG guides, in particular guide 7's "Coordinate Conversions and Transformations including Formulas", as it has a very comprehensive listing of different mapping systems, and the algorithms needed to convert them to other coordinate systems.

Cartesian coordinates

GPS systems normally show their coordinates in latitude and longitude format. However, behind the scenes, they actually calculated the position using Cartesian X,Y,Z coordinates. These give the location of a point in metres (normally) from the centre of the ellipsoid/planet. X is the distance towards the prime meridian of the ellipsoid's datum (towards 0°), Y is the distance towards the 90° meridian, and Z is the distance towards the North Pole. This makes it very easy to perform calculations like the direct line of sight distance between two points, since the distances can be directly related to each other.

At all points around the globe, physical distances are measures with the same resolution. This contrasts with latitude and longitude, where longitude bunches up at the poles such that moving a few cm changes your longitude by 180°, while at the equator it would take 20000 km to change longitude by that much. However, it can seem very odd that moving horizontally across the Earth causes a shift in not just X and Y, but Z too, and the coordinates therefore can be difficult to comprehend. Several tools may operate in X,Y,Z coordinates, so these scripts provide methods to easily convert between them and latitude,longitude.

The conversion relies on knowing how big the ellipsoid is, in order to convert from coordinates based only on degrees, to coordinates that are given in metres. Therefore it needs to know which ellipsoid your original coordinates relate to. For regular GPS, this is WGS84. For differential GPS in Europe, it is normally ETRS89. The resulting XYZ coordinates also relate to that same ellipsoid, but only the datum part which says where the prime meridian is located, and therefore which direction X relates to.

PHP
$grutoolbox = Grid_Ref_Utils::toolbox();
$WGS84 = $grutoolbox->get_ellipsoid('WGS84');
$source_coords = Array(54.607720,-6.411990);
//height above the ellipsoid
$source_height = 66.341;
//convert to X,Y,Z
$xyz = $grutoolbox->lat_long_to_xyz($source_coords,$source_height,$WGS84);
//do something with the results
$some_other_xyz = ... etc. ...;
//convert back to latitude,longitude
$final_coords = $grutoolbox->xyz_to_lat_long($some_other_xyz,$WGS84);
JavaScript
var grutoolbox = gridRefUtilsToolbox();
var WGS84 = grutoolbox.getEllipsoid('WGS84');
var sourceCoords = [54.607720,-6.411990];
//height above the ellipsoid
var sourceHeight = 66.341;
//convert to X,Y,Z
var xyz = grutoolbox.latLongToXyz(sourceCoords,sourceHeight,WGS84);
//do something with the results
var someOtherXyz = ... etc. ...;
//convert back to latitude,longitude
finalCoords = grutoolbox.xyzToLatLong(someOtherXyz,WGS84);

Polynomial transformation

Polynomial transformation is a superior alternative to Helmert transformation. It can be used to transform latitude,longitude coordinates from one ellipsoid to another, but at the same time, it can take into account the way that gravity varies over an area and the distortions that occur within it, which is something that a Helmert transformation cannot do. Distortions occur because of changes in the density of the Earth's mantle over distances, and also because of the additional mass of nearby mountains. These can cause an area to to be level according to gravity, while a GPS might see it as being a slope. The result is that horizontal distances, which is what mapping cares about, can be more or less than the distance around an ellipsoid. In reality, these differences produce a few metres of error in various directions across in the British Isles.

The more coefficients that can be used, the more distortion fluctuations that can be represented. The scripts support third order polynomials, with 16 coefficient pairs. The primary use is to convert to Irish Grid, where the coefficients represent 183 different fluctuations that have been measured across the island. The conversion does not need to know which ellipsoid the coordinates relate to, but you need to make sure you are using them to convert in the correct direction from one ellipsoid to another.

The following example uses the built-in coefficients to translate from Irish grid references to ETRS89 coordinates using PHP. This is can be done much more simply by using the COORDS_GPS_IRISH type with grid_to_lat_long/gridToLatLong, which uses polynomial transformation automatically. However, it is shown here expanded into a few steps so you can see how the process works in cases where you might need to use the intermediate values. Of course, most GPS systems give results in WGS84 not ETRS89, so there will always be the continental drift error if those coordinates are converted directly. However, a 40 cm + 79 cm error (in 2020) is still far better than the 5 m + 79 cm error when using Helmert transformation. Note that if your software uses the built-in 'OSiLPS' coefficients, then you will need to check the copyright section above.

$grutoolbox = Grid_Ref_Utils::toolbox();
//convert Irish Grid reference to modified Airy 1830 coordinates
$irish_grid_reference = 'J0259874444';
$irish_grid_numbers = $grutoolbox->get_Irish_grid_nums($irish_grid_reference);
$ma1830_coords = $grutoolbox->grid_to_lat_long($irish_grid_numbers,$grutoolbox->COORDS_OSI);
//convert modified Airy 1830 coordinates to ETRS89 coordinates
$coefficients = $grutoolbox->get_polynomial_coefficients('OSiLPS');
$etrs89_coordinates = $grutoolbox->polynomial_transform($ma1830_coords,$coefficients);

This JavaScript example creates a custom set of third order polynomial coefficients, and converts in reverse using them. In this case, these coefficients might be used to apply or remove a correctional shift without transforming it between ellipsoids:

var grutoolbox = gridRefUtilsToolbox();
var sourceCoords = [54.607720,-6.411990];
//define coefficients
var coefficients = grutoolbox.createPolynomialCoefficients(
  3,
  true,
  54.721448,
  -6.860962,
  0.3,
  [[0.03,0.0011,0.012,-0.00263],[0.0131,0.002,-0.0324,-0.00121],[0.12,-1.75,0.211,1.3],[0.412,-2.188,7.456,-1.951]],
  [[0.0038,0.00029,0.00737,0.0175],[0.00753,-0.035,0.758,-0.622],[1.275,-0.25,-2.348,1.255],[1.112,0.7332,-0.8465,0.535]]
);
//shift the coordinates
var shiftedCoords = grutoolbox.reversePolynomialTransform(sourceCoords,coefficients);

It is worth noting that polynomial transformation ignores heights completely. With the Irish Grid, accurate conversion of heights is done using geoids instead, or they can be approximated with Helmert transformations.

Basic geodesics

Geodesics serve two purposes. The first is to work out where you would end up if you travelled in a given direction from a point for a certain distance around the curve of the Earth. The second is to work out the distance and direction between two points around the curve of the Earth. They are based only on horizontal distance around an ellipsoidal Earth model, and ignore heights. The scripts make this extremely easy.

The following example for PHP shows how to work out where you would end up after travelling from a point. If working with GPS coordinates, you will want the WGS84 ellipsoid.

$grutoolbox = Grid_Ref_Utils::toolbox();
$startpoint = array( 54.6078031260, -6.4129391413 );
print $grutoolbox->get_geodesic_destination( $startpoint, 1270.119, 137.24, $grutoolbox->get_ellipsoid('WGS84'), $grutoolbox->HTML );

The following example for JavaScript works out the distance and azimuth between two points:

var grutoolbox = gridRefUtilsToolbox();
var startpoint = [ 54.6078031260, -6.4129391413 ];
var endpoint = [ 51.884008219444446, -3.4364607166666667 ];
document.write( grutoolbox.getGeodesic( startpoint, endpoint, grutoolbox.getEllipsoid('WGS84'), grutoolbox.HTML ) );

Geodesics have certain limitations, and there are also ways they can be used to solve more complex problems. These are discussed in more detail in the more about geodesics section below.

Accepting input from users

Some of the methods provided by the scripts can accept string input, which can be used to interpret user-submitted input, which may or may not be in a valid format. In all cases, these methods accept an extra parameter (called deny_bad_*), which can be used to determine if the script was able to recognise it as the given coordinate type, and if not, try another one.

Note that there are some cases where dms_to_dd/dmsToDd may return a result, when it may be more appropriate to interpret the string as generic numeric grid reference, such as where the string contains '<number> ,<number> ' (with whitespace after the numbers) or '<number>m, <number>m'. As a result, the following example checks if it is a generic numeric grid reference first. Similarly, it may also accept '1U, 123 456 ' (with a trailing space), even though that may be more appropriately interpreted as a UTM grid reference, so that check is also placed beforehand.

If you are expecting UPS grid references as input, you can also use ups_to_lat_long/upsToLatLong to detect those. However, the normal UPS format is identical to the Irish grid format, with the only exception being that within the normal range, the numbers are usually two pairs of at least 6 digits, while with Irish grid references, they are normally only upto 5 digits. To detect this difference, the ups_to_lat_long/upsToLatLong method accepts an additional min_length parameter, which causes the method to only accept a string as a UPS grid reference if the easting and northing are at least 800000 (a number slightly below the minimum expected values within the UPS grids). This check must be made before the check for Irish grid references.

PHP
$coordsstring = $_GET['coords'];
if( function_exists('get_magic_quotes_gpc') && get_magic_quotes_gpc() ) { $coordsstring = stripslashes($coordsstring); }
$grutoolbox = Grid_Ref_Utils::toolbox();
if( $coords = $grutoolbox->parse_grid_nums($coordsstring,$grutoolbox->DATA_ARRAY,true) ) {
  //it was interpreted as a generic numeric grid reference
  //assume it was ITM and convert to lat/long
  $coords = $grutoolbox->grid_to_lat_long($coords,$grutoolbox->COORDS_GPS_ITM);
} elseif( $coords = $grutoolbox->utm_to_lat_long($coordsstring,null,$grutoolbox->DATA_ARRAY,true) ) {
  //it was interpreted as a UTM grid reference
  //the method converts it to GPS coordinates; nothing else is needed
} elseif( $coords = $grutoolbox->ups_to_lat_long($coordsstring,$grutoolbox->DATA_ARRAY,true,true) ) {
  //it was interpreted as a UPS grid reference
  //the method converts it to GPS coordinates; nothing else is needed
} elseif( $coords = $grutoolbox->dms_to_dd($coordsstring,$grutoolbox->DATA_ARRAY,true) ) {
  //it was interpreted as latitude/longitude coordinates
  //assume it was GPS coordinates; nothing else is needed
} elseif( $coords = $grutoolbox->get_UK_grid_nums($coordsstring,$grutoolbox->DATA_ARRAY,true) ) {
  //it was interpreted as a UK grid reference - convert to lat/long
  $coords = $grutoolbox->grid_to_lat_long($coords,$grutoolbox->COORDS_GPS_UK);
} elseif( $coords = $grutoolbox->get_Irish_grid_nums($coordsstring,$grutoolbox->DATA_ARRAY,true) ) {
  //it was interpreted as an Irish grid reference - convert to lat/long
  $coords = $grutoolbox->grid_to_lat_long($coords,$grutoolbox->COORDS_GPS_IRISH);
} else {
  //failed to recognise it as any format at all
}
JavaScript
var coordsString = document.getElementById('coordinateinput').value;
var coords, grutoolbox = gridRefUtilsToolbox();
if( coords = grutoolbox.parseGridNums(coordsString,grutoolbox.DATA_ARRAY,true) ) {
  //it was interpreted as a generic numeric grid reference
  //assume it was ITM and convert to lat/long
  coords = grutoolbox.gridToLatLong(coords,grutoolbox.COORDS_GPS_ITM);
} else if( coords = grutoolbox.utmToLatLong(coordsString,null,grutoolbox.DATA_ARRAY,true) ) {
  //it was interpreted as a UTM grid reference
  //the method converts it to GPS coordinates; nothing else is needed
} else if( coords = grutoolbox.upsToLatLong(coordsString,grutoolbox.DATA_ARRAY,true,true) ) {
  //it was interpreted as a UPS grid reference
  //the method converts it to GPS coordinates; nothing else is needed
} else if( coords = grutoolbox.dmsToDd(coordsString,grutoolbox.DATA_ARRAY,true) ) {
  //it was interpreted as latitude/longitude coordinates
  //assume it was GPS coordinates; nothing else is needed
} else if( coords = grutoolbox.getUKGridNums(coordsString,grutoolbox.DATA_ARRAY,true) ) {
  //it was interpreted as a UK grid reference - convert to lat/long
  coords = grutoolbox.gridToLatLong(coords,grutoolbox.COORDS_GPS_UK);
} else if( coords = grutoolbox.getIrishGridNums(coordsString,grutoolbox.DATA_ARRAY,true) ) {
  //it was interpreted as an Irish grid reference - convert to lat/long
  coords = grutoolbox.gridToLatLong(coords,grutoolbox.COORDS_GPS_IRISH);
} else {
  //failed to recognise it as any format at all
}

Try it here; input any one of the formats mentioned above:

You will need JavaScript support to test this script.

Note that deny_bad_* does not enforce a very strict restriction on the format to ensure that it is a perfectly valid sequence. There are cases where you could feed invalid data to the method, and it will accept it, if it can interpret it as a potential match for the given type of data that it was expecting. For example, the string '796.345( 5632.222.g N, 999 E' is not a valid format for latitude/longitude, but the scripts will be able to extract meaning from it - deny_bad_coords will therefore not cause the method to return false in that case. The purpose of deny_bad_* is to help your script make sense out of the many potential types of data that users may feed the scripts, not to complain to them if they don't give a perfect string format.

The parse_grid_nums/parseGridNums method is an exception here, as it has a strict_nums parameter that allows you to force it to accept only a limited number of formats. This exists only for the purpose of determining if a set of coordinates are given with units that make them most likely to be ITM grid coordinates, or if the coordinates are simply generic coordinates that could apply to any of the numeric grid reference formats. You can choose to use this functionality at your own discretion.

Likewise, the dms_to_dd/dmsToDd method has the allow_unitless parameter that makes it accept simple number pairs, separated by a comma. These number pairs cannot be distinguished from grid coordinates, so if you want to allow both, you will need to try parsing as grid coordinates first, and reject it if the grid coordinates are below 180, or negative. It is important to note, however that '56.78,12.34' is valid for both grid coordinates and GPS coordinates, so you will need to be careful about which numbers will be considered valid for each of them. Note that dms_to_dd/dmsToDd will accept any numbers, even those too large to be valid, and will wrap the longitude around to the right range. This means it cannot be used to detect whether the numbers are too large.

Users may still input values that extend far beyond the expected grids. While the scripts will force longitude to wrap-around when converting from Transverse Mercator or polar stereographic grid to latitude/longitude, they do not attempt to prevent out-of-bounds latitudes, and will return whatever values the conversion algorithms produce. This can lead to unexpectedly high GPS coordinates, for example, potentially enough to travel several quadrillion times around the World. You may wish to check the returned values are within the expected limits.

The exceptions here are shift_grid/shiftGrid, reverse_shift_grid/reverseShiftGrid and apply_geoid/applyGeoid, which also have a deny_bad_coordinates parameter. In these cases, it will refuse to accept coordinates which are outside the grid, because there is no data covering that area. This allows the situation to be detected, so that warnings can be shown. The methods will also return false if the data could not be loaded, for whatever reason.

Similarly, lat_long_to_utm/latLongToUtm and lat_long_to_ups/latLongToUps have a deny_bad_coordinates parameter to allow you to detect when the coordinates are outside the coordinate system region, so you can switch to the other system instead. These coordinate systems are designed to work together. UTM can represent coordinates (badly) inside the polar regions using only a band letter, but UPS cannot represent coordinates outside the polar regions at all. In this case, the parameter is not used for detecting invalid input strings.

Geoids

The Earth's shape is not perfect, and gravity is influenced both by the mass of the crust, and the density of the mantle. Mountains have a significant affect on gravity and increase it noticeably, but the mantle has the most substantial effect. Diverging tectonic plate boundaries increase gravity the most, with the strongest effects being near Iceland and Indonesia. The most significant low points are in the World's oceans, with the Indian Ocean and Sri Lanka having the lowest gravity of any part of the planet. The effect of this is to alter the mean sea level at that point, with the strongest gravity raising the sea level significantly. Ignoring tides, the sea level at Sri Lanka is 192 metres lower than nearby in Indonesia, compared with what would be expected on a perfectly smooth spheroid. This is why in the UK, a GPS height is about 50 metres higher than the local mapping (orthometric) height; sea level in the UK is 50 metres higher than the World average. It also means that almost all mountains are actually shorter compared with the surrounding terrain in orthometric height than they are in ellipsoidal height.

A geoid represents the vertical gravity effect over an area of the planet. A geoid can be thought of as the sea level, if the land were ignored, tides were removed, and ocean currents were stopped. By subtracting the geoid undulation (the difference between the geoid height and the ellipsoid height), you can get an idea of the local terrain height. It still normally will not take the fluctuations caused by tidal currents and ocean currents into account, so it is not perfect, but still better than nothing. It gives a better idea of whether you would feel like you were walking uphill, and whether streams would flow in one direction or the other, while ellipsoid height can give the wrong impression (within the UK, the ellipsoid height would suggest that rivers cannot slowly flow west, but the Rivers Severn and Clyde would disagree). The OSGM15 geoid is actually calibrated to use the local mapping sea level datums, so it does give real local mapping heights.

Geoids are calculated by precise gravimetric scans of the area, typically using satellites. These are then normally converted into mathematical formulas. These or the original data can be used to export a grid of data points covering the area. A common file format for these data points is the .grd (grid) file format, and these scripts are designed to work with software that interprets the .grd file format. These files can be very large. EGM96, a very popular global geoid, has a low density version called ww15mgh.grd, which is 9 MB. Higher density data covering just Ireland is just over 1 MB.

The scripts expect data to be provided by an external database, with your own scripts communicating with the database and providing the data when needed. The data is expected to be provided in a regular latitude-longitude grid of points, typically supplied by the .grd file format, and the scripts will ask for whichever data points they need for the current conversion. Firstly, the scripts need to know the specification for the grid, which can all be found in a .grd file's header. When calling the apply_geoid/applyGeoid method, it must be given a record fetching callback function which it can use to request data. If it needs data, it will pass the callback function an array of the data points it needs, giving the latitude, longitude, latitude index, longitude index, and record index of the data point within the grid. This allows the record fetching callback function to choose whichever method it would like, to identify the data points to return. The callback must return the data points that it got, in an array of geoid grd data point records; for that, it will need the data point value and the record index, for each point.

Processing the .grd file is fairly simple due to it being a very basic format. The first line contains 6 numbers which are used to make the geoid grd specification. Every line after that contains several data points, which can be assigned an increasing record index starting from 0, beginning at the northwest corner of the grid, and progressing eastwards, one row at a time. If you wanted to work with latitude and longitude indexes, or actual coordinates, the first 6 numbers are the minimum latitude, maximum latitude, minimum longitude, maximum longitude, latitude spacing and longitude spacing. It seems like fairly simple mathematics to work out how many latitude rows there are, and how many longitude columns (watch out for the fencepost problem). However, note that files may use truncated recurring numbers as the spacing, and so the numbers may not divide perfectly, and you need to work out the intent rather than the actual specified numbers. If it helps, I have also released a .grd file parser for PHP, which can help to populate your database, and help you to visualise the data.

The following example shows how to get from ellipsoid height to geoid height, using the predefined EGM96_ww15mgh geoid:

function load_geoid_records( $name, $records ) {
  global $grutoolbox;
  $whereclause = array();
  foreach( $records as $onerecord ) {
    //this example database uses recordindex
    $whereclause[] = "recordindex = '" . $onerecord['recordindex'] . "'";
  }
  $whereclause = implode( ' OR ', $whereclause );
  $results = database_query_command( "SELECT * FROM `" . addslashes($name) . "`  WHERE $whereclause" );
  if( database_row_count($results) ) {
    $resultarray = array();
    while( $row = database_fetch_array($results) ) {
      $resultarray[] = $grutoolbox->create_geoid_record( $row['recordindex'], $row['shift'] );
    }
    return $resultarray;
  } else {
    //database data did not match the expected grid
    return array();
  }
}
print $grutoolbox->apply_geoid( 54.6078, -6.4129, 66.341, true, $grutoolbox->get_geoid_grd('EGM96_ww15mgh'), 'load_geoid_records', $grutoolbox->HTML );

The scripts will cache every returned value automatically, so that subsequent requests that use the same data points will not need to load the data dynamically. It is also possible to load parts of the data first, create geoid grd data point records for them, and use cache_geoid_records/cacheGeoidRecords to cache them. Once data has been cached, it continues to use up memory, which can become a problem if large amounts of data gets loaded. If this is likely to be an issue, the delete_geoid_cache/deleteGeoidCache method may be used to delete the cached data for one geoid or all geoids.

Because the files are generally too large to transmit in full as Web page data, it is expected that JavaScript will need to use dynamic requests to a server, and as such, it will probably need to operate asynchronously. The JavaScript version of the script therefore offers a return callback feature that can be used to run the code asynchronously if needed. If it is not provided, the scripts will operate synchronously and return the values directly, just like the PHP version. Synchronous operation may be suitable for very small geoid grids, for JavaScript applications, for scripts that cache the data first, or for scripts that run in a JavaScript (Web) worker using synchronous XMLHttpRequest to load the results. Note that while it is currently possible to use synchronous XMLHttpRequest from a normal JavaScript thread, it is strongly discouraged as it locks up the page while waiting for a response, and browsers may stop supporting it in future.

If the return callback function is provided, then the method will call that function with the result instead, and asynchronous operation will be enabled. The record fetching callback function also gets given a third parameter; a database callback function. When the results have loaded, the database callback function must be called with the results, instead of returning the results directly. Make sure that you do not delete the geoid cache until the final result has been returned, since it will be operating asynchronously and may still need to use the cache. (You may note that this relies on callbacks instead of promises, async functions or the convenient await syntax. This allows the script to work with popular legacy browsers such as Internet Explorer.)

Basically, for synchronous operation, call the method, and it returns the results. For asynchronous operation, call the method and pass it a callback function. In response, it will pass you a callback function when it asks for data. You use that callback to pass back the data, and it will deliver the results to your callback.

The following example shows how to use both synchronous and asynchronous operation with JavaScript. The functions can also be implemented as anonymous functions, if needed:

function loadGeoidRecordsSynchronous( name, records ) {
  //implement this dtabase connection yourself
  ... do something to get the data synchronously ...
  ... parse the response, such as JSON.parse ...
  var newarray = [];
  ... for each record in the results ...
    newarray.push( grutoolbox.createGeoidRecord( oneResult.index, oneResult.shift ); );
  ... then ...
  //synchronous operation, return directly
  return newarray;
}
var result = grutoolbox.applyGeoid(
  54.6078, -6.4129, 66.341,
  false,
  grutoolbox.getGeoidGrd('OSGM15_Belfast'),
  loadGeoidRecordsSynchronous,
  grutoolbox.DATA_ARRAY
);

function loadGeoidRecordsAsynchronous( name, records, dataCallback ) {
  //implement this dtabase connection yourself
  ...do something to get the data asynchronously, such as asynchronous XMLHttpRequest ...
  dataFetcher.onload = function (responseText) {
    ... parse the response, such as JSON.parse ...
    var newarray = [];
    ... for each record in the results ...
      newarray.push( grutoolbox.createGeoidRecord( oneResult.index, oneResult.shift ); );
    ... then ...
    //asynchronous operation, use the callback
    dataCallback(newarray);
  }
}
function handleResults( result ) {
  ... data has been processed, use the result ...
  ... geoid cache can be deleted if needed ...
}
var anotherResult = grutoolbox.applyGeoid(
  54.6078, -6.4129, 66.341,
  false,
  grutoolbox.getGeoidGrd('OSGM15_Malin'),
  loadGeoidRecordsAsynchronous,
  grutoolbox.DATA_ARRAY,
  false,
  handleResults
);

Bilinear interpolated grid shifting

This is the highest accuracy method for shifting coordinates between earth models, taking into account all of the distortions caused by gravity fluctuations, and optionally applying geoid corrections at the same time. Unlike polynomial transformation, it has no limits to how many different corrections it can apply, and it can be given as much or as little detail as needed, to cater to the specific needs of a particular situation. Grid shifts can be used to convert between one idealised Earth model to another, or to correct inaccuracies in an Earth model, or do both at once.

The primary use for this is the simultaneous OSTN15 and OSGM15 transformation used by the British National Grid, to convert between ETRS89 GPS coordinates and OSGB36 National Grid references and orthometric heights, bypassing the need to convert to Airy 1830 coordinates first. This transformation is so accurate that it is now used along with the OS Net GNSS base stations to define the National Grid itself. When tested against land-based surveying techniques, there was found to be a 10 cm RMS distortion, most of which comes down to limitations in land based surveying techniques. This means that when you start with ETRS89 coordinates, and convert to OSGB36 using OSTN15, no approximations are used in the transformation itself, and the transformation is 100% compatible with the National Grid. Of course, your GPS readings can still have their own precision errors or accuracy errors. The grid corrections are given at a resolution of one data point per km. The data is available as part of the Ordnance Survey OSTN15 developer pack, in a 40 MB file called OSTN15_OSGM15_DataFile.txt, which is in CSV format. It is perhaps best to use a proper database to store that data, due to its size, and due to the speed at which a database can retrieve results without using too much memory. If it helps, I have also released an OSTN15 data file parser, which can help to populate your database from the Ordnance Survey data file.

Just like when dealing with geoids, the scripts expect data to be provided by an external database, with your own scripts communicating with the database and providing the data when needed. The data is expected to be provided in a regular easting,northing grid of points, and the scripts will ask for whichever data points they need for the current conversion. Firstly, the scripts need to know the shift grid set; the specification for the grid. The OSTN15 shift grid set is provided by default, so you don't have to look up the details. When calling the shift_grid/shiftGrid and reverse_shift_grid/reverse_shiftGrid methods, they must be given a record fetching callback function which it can use to request data. If they need data, they will pass the callback function an array of the data points they need, giving the easting and northing for each. If the shift grid set specifies the number of columns, then the callback function will also be given the record index of the data point within the grid, which the fetching callback can use to look up records. However, this is not mandatory, as a shift grid may have uneven numbers of columns, and cover an irregularly shaped area. In all cases, the easting and northing will always be provided, and will be used internally to recognise data points. For this reason, shift grids should, if at all possible, use integers for the coordinates and spacing - they are given to a resolution of 1 metre and it is very unlikely that a grid will need spacing that is not an exact number of metres. This is not mandatory though, and it might work for some floating point numbers. However, it will almost certainly run into floating point precision problems recognising records, if you try to use grids with spacings like 0.333... The callback must return the data points that it got, in an array of bilinear interpolation shift records; for that, it will need each data point's east, north and vertical shift values, as well as their easting and northing.

When reversing the transformation, it is important to note that the positions of the data points relates to the unshifted coordinates, while the coordinates that you supply will be the shifted coordinates. To reverse the transformation, the scripts have to try using the current position's shift, apply it, and then re-check the forward shift from that new position, until it finds the point which shifts to the right position. During this process, it might need to load more and more records, if the shift takes it outside the original grid square. In the worst case, it might need to call the record fetching callback function 4 times for different sets of coordinates (or more if the shifts are larger than the grid spacing), since it cannot predict this in advance. This will make it slower than a forwards transformation.

To get from ETRS89 GPS coordinates to OSGB36 grid references, you need to firstly project the coordinates from latitude,longitude, into ETRS89 grid coordinates using the OSGB36 datum. This positions them in the wrong place, since the two grids use different datums. Then you need to apply the OSTN15 shift to move the coordinates into the right location for the OSGB36 grid. To make this easier for you, the scripts provide the ETRS89+OSGB36 datum as a datum called 'OSTN15'. The record fetching callback function can also return height datum information, which allows you to display a warning if the shifting relied on more than one datum (in which case the height results can be nonsense). The following PHP example shows how to get from ETRS89 latitude and longitude to OSGB36 grid reference:

function load_shift_records( $name, $records ) {
  global $grutoolbox;
  $whereclause = array();
  foreach( $records as $onerecord ) {
    //this example database uses recordindex
    //$whereclause[] = "easting = '" . $onerecord['easting'] . "' AND northing = '" . $onerecord['northing'] . "'";
    $whereclause[] = "recordindex = '" . $onerecord['recordindex'] . "'";
  }
  $whereclause = implode( ' OR ', $whereclause );
  $results = database_query_command( "SELECT * FROM `" . addslashes($name) . "`  WHERE $whereclause" );
  if( database_row_count($results) ) {
    $resultarray = array();
    while( $row = database_fetch_array($results) ) {
      $resultarray[] = $grutoolbox->create_shift_record( $row['easting'], $row['northing'], $row['se'], $row['sn'], $row['sg'], $row['datum'] );
    }
    return $resultarray;
  } else {
    //database data did not match the expected grid
    return array();
  }
}
$etrs89_gps = array(51.8840136123,-3.4364538123);
$etrs89_height = 938.946;
$etrs89_grid = $grutoolbox->lat_long_to_grid( $etrs89_gps, $grutoolbox->get_datum('OSTN15'), $grutoolbox->DATA_ARRAY );
print $grutoolbox->shift_grid( $etrs89_grid, $etrs89_height, $grutoolbox->get_shift_set('OSTN15'), 'load_shift_records', $grutoolbox->HTML );

This would perform the reverse conversion, using the same load_shift_records function:

$osgb36_grid = array( 12345, 67890 );
$osgb36_height = 938.946;
$etrs89_grid = $grutoolbox->reverse_shift_grid( $os_grid, $osgb36_height, $grutoolbox->get_shift_set('OSTN15'), 'load_shift_records', $grutoolbox->DATA_ARRAY );
print $grutoolbox->grid_to_lat_long( $etrs89_grid, $grutoolbox->get_datum('OSTN15'), $grutoolbox->HTML );
print ', height ' . $etrs89_grid[2];

Just like with geoids, the scripts will cache every returned value automatically, so that subsequent requests that use the same data points will not need to load the data dynamically. It is also possible to load parts of the data first, create bilinear interpolation shift records for them, and use cache_shift_records/cacheShiftRecords to cache them. Once data has been cached, it continues to use up memory, which can become a problem if large amounts of data gets loaded. If this is likely to be an issue, the delete_shift_cache/deleteShiftCache method may be used to delete the cached data for one shift grid or all shift grids.

Because the files are generally too large to transmit in full as Web page data, far larger than a typical geoid, JavaScript will almost certainly need to use dynamic requests to a server, and it will probably need to operate asynchronously. The JavaScript version of the script therefore offers a return callback feature that can be used to run the code asynchronously if needed. If it is not provided, the scripts will operate synchronously and return the values directly, just like the PHP version. Synchronous operation may be suitable for very small shift grids, for JavaScript applications, for scripts that cache the data first, or for scripts that run in a JavaScript (Web) worker using synchronous XMLHttpRequest to load the results. Note that while it is currently possible to use synchronous XMLHttpRequest from a normal JavaScript thread, it is strongly discouraged as it locks up the page while waiting for a response, and browsers may stop supporting it in future.

If the return callback function is provided, then the method will call that function with the result instead, and asynchronous operation will be enabled. The record fetching callback function also gets given a third parameter; a database callback function. When the results have loaded, the database callback function must be called with the results, instead of returning the results directly. Make sure that you do not delete the shift cache until the final result has been returned, since it will be operating asynchronously and may still need to use the cache. This is especially important when reversing the shift transformation, as it will need to use the shift cache several times for each operation. (You may note that this relies on callbacks instead of promises, async functions or the convenient await syntax. This allows the script to work with popular legacy browsers such as Internet Explorer.)

Basically, for synchronous operation, call the method, and it returns the results. For asynchronous operation, call the method and pass it a callback function. In response, it will pass you a callback function when it asks for data. You use that callback to pass back the data, and it will deliver the results to your callback.

The following example shows how to use both synchronous and asynchronous operation with JavaScript. The functions can also be implemented as anonymous functions, if needed:

function loadShiftRecordsSynchronous( name, records ) {
  //implement this dtabase connection yourself
  ... do something to get the data synchronously ...
  ... parse the response, such as JSON.parse ...
  var newarray = [];
  ... for each record in the results ...
      newarray.push( grutoolbox.createShiftRecord( oneResult.easting, oneResult.northing, oneResult.eastshift, oneResult.northshift, oneResult.heightshift, oneResult.datum ); );
  ... then ...
  //synchronous operation, return directly
  return newarray;
}
var etrs89Gps = [51.8840136123,-3.4364538123];
var etrs89Height = 938.946;
var etrsGrid = grutoolbox.latLongToGrid( etrs89Gps, grutoolbox.getDatum('OSTN15'), grutoolbox.DATA_ARRAY );
var result = grutoolbox.shiftGrid(
  etrs89Gps, etrs89Height,
  grutoolbox.getShiftSet('OSTN15'),
  loadShiftRecordsSynchronous,
  grutoolbox.DATA_ARRAY
);

function loadShiftRecordsAsynchronous( name, records, dataCallback ) {
  //implement this dtabase connection yourself
  ...do something to get the data asynchronously, such as asynchronous XMLHttpRequest ...
  dataFetcher.onload = function (responseText) {
    ... parse the response, such as JSON.parse ...
    var newarray = [];
    ... for each record in the results ...
      newarray.push( grutoolbox.createShiftRecord( oneResult.easting, oneResult.northing, oneResult.eastshift, oneResult.northshift, oneResult.heightshift, oneResult.datum ); );
    ... then ...
    //asynchronous operation, use the callback
    dataCallback(newarray);
  }
}
function handleResults( result ) {
  ... data has been processed, use the result ...
	var etrsLatLong = grutoolbox.gridToLatLong( etrs89Gps, grutoolbox.getDatum('OSTN15'), grutoolbox.DATA_ARRAY );
	var etrsHeight = result[2];
  ... shift cache can be deleted if needed ...
}
var reversed = grutoolbox.reverseShiftGrid(
  result,
  grutoolbox.getShiftSet('OSTN15'),
  loadGeoidRecordsAsynchronous,
  grutoolbox.DATA_ARRAY
  false,
  false,
  handleResults
);

For those who are interested, the way that reverse iteration is possible both in synchronous and asynchronous operation, is to use a scoped function which calls itself recursively intead of using a loop for iteration (this allows an asynchronous callback to run the scoped function after fetching data). This has the limitation that recursive functions can cause stack overflows. In this case, it is limited to 100 iterations, and rarely uses more than 5. This is well within the stack sizes of even most legacy browsers - modern browsers have a stack size of several thousand.

More about geodesics

The scripts offer methods that can solve two problems: adding a geodetic vector to a point to calculate the resulting position, and finding the distance and bearing of the shortest route between two points around the curve of the World - a geodesic. These are known as the Direct Problem and Inverse Problem respectively. There are several methods of doing each of these, some of which only work on spheres, even though the World is better represented as a spheroid. The most complete approach (used by GeographicLib, if you wanted to try it) is a very complex and extremely lengthy one developed by Charles Karney (2012) which uses several methods at once to obtain a nearly perfect result in all cases (it can also tell you where two vecors collide, and how far north a vector will travel). However, there is a much more simple (relatively speaking) method that can work in the vast majority of cases, with minimal errors (typically less than 5 mm), designed to be used in places where computation efficiency is important. That is what these scripts use.

Geodesics are calculated using a modified version of the Vincenty formulae, inspired by the works and modifications of Chris Veness, but reworked to cater to a few more situations. The modifications are intended to cope with the less common situations where the Vincenty iteration does not converge, such as points on opposite sides of the World from each other, or points that are positioned extremely close to each other. These are detected either by the the trigonometry functions returning values too small to be added to larger numbers (a limitation of floating point arithmetic), or by the formulas requiring too many iterations.

The error handling for problematic coordinates tries to correct the returned distances and azimuths, with varying degrees of accuracy. For points on opposite longitudes from each other, or for points where at least one is situated on a pole, it can calculate the correct azimuth. For points that are nearly opposite, it gives the most likely direction, but with a very tiny error (which normally gets lost in rounding precision anyway). For points extremely close to each other, it returns a poor approximation of the azimuth between them (since at those tiny distances, trigonometry functions cannot return anything useful), within 45°. However, this error gets worse closer to the poles, and is completely meaningless right beside the poles. This is really quite an unimportant correction given that the distance between the points is so close to 0 that it gets rounded to 0 by the rounding precision anyway.

When error handling has been used, an extra parameter will be returned in any data arrays, saying how the problem was detected, the type of separation between the points, and the confidence in how well it has approximated the distances and azimuths. The error handling approximates the returned distance, in many cases just giving the antipodal distance for points that are not quite antipodal, or 0 for points that are very close to each other, so the mere presence of this return parameter indicates that the distance data is not entirely trustworthy. You can use this to display notices when the distances may be inaccurate, or when the confidence in the azimuths is particularly bad. When returning text or HTML strings, this is signified with an approximation symbol in the returned string.

The azimuth - which is also known as a bearing in British English - is the direction of the geodetic vector as seen from the start or end points. At the poles, this may seem like nonsense because every direction from the North Pole is southwards, and vice versa. However, with a geodesic, the azimuth will be the direction compared with the line of longitude which is given for a point. This means that although 90°,0° and 90°,123° are the same point (both are the North Pole), they will produce different azimuths, as the first compares the direction with that of the prime meridian, and the second compares the direction with that of the 123° meridian. For 90°,0°, an azimuth of 0° would continue "forwards" down the 180° meridian, while for 90°,123°, an azimuth of 123° would point down the 180° meridian. For -90°,0°, an azimuth of 0° would point up the prime meridian. Polar longitude is relevant with geodesics.

The distance between two mountain tops

Geodetic distances are the distances around the ellipsoid being used, which is a rough approximation of the shape of the World, ignoring slopes and imperfections in the shape of the World. Since they are being measured around a curve, the distances are further for points that are higher than that ellipsoid, either because they are high above sea level, or because they are on a part of the World where sea level is higher than the ellipsoid being used. The distance of a geodesic must therefore be treated as only a rough approximation, and not the distance that you would need to travel in order to follow it on the actual World's surface. Following the slopes and bumps of the surface would add a lot more distance, and a direct line of sight between points is shorter than the curved distance between those same points. This should put that <= 5 mm error into perspective; no matter how perfect your formula is, the answer can never be more than a simple guide. Even for navigation at sea where the surface is relatively flat and consistent, the ellipsoid is still only an approximation compared with a geoid (gravitational representation of the World), but there are no perfect formulae for determining the shortest route between points on a geoid. Geodesics still have their uses - as long as you accept the limitations - such as approximating the length of a GPS track by adding the geodesic distances between pairs of points, or approximating the shortest direction between distant islands.

Mountain tops are typically raised some distance above the WGS84 ellipsoid (some will be under it, depending on how poorly WGS84 represents the area). Assuming you know the height of your points compared with the reference ellipsoid (such as a GPS altitude, which is the height above the WGS84 ellipsoid), you can use the height above the ellipsoid as a circle radius and a simple 2πr formula to calculate how much further the distance should be over a full circle. For points 100 m above the ellipsoid, the difference over a full circle is 628 m. Assuming the World is a sphere, to get the distance adjustment for your arc, you would need to multiply it by the length of your arc divided by the length of a full great circle. For two hills 10 km apart that rise 100 m above a sphere based on the WGS84 major axis, the tops are over 15 cm further apart than the bases which sit perfectly on the sphere. This is how far an aircraft would have to fly to get between the tops if it maintained a constant altitude above the sphere. The direct line of sight between the two hills is under a tenth of a millimetre longer than 10 km, yet another way to measure distances. To an aircraft, this would feel like it was losing altitude for half the journey, then gaining altitude for the other half.

However, since the World is better represented as an oblate spheroid, the result will be slightly wrong, and the exact result should be different at different points and orientations on the ellipsoid. It is, however, possible to get the distance by supplying a different ellipsoid. If the two points sit 100 m above the WGS84 ellipsoid, supply an ellipsoid with the major and minor axes 100 m larger than WGS84's. The resulting geodesic will be the correct length for the curved path between those points.

If the points are at different heights, then simple Pythagorean theorem (hypotenuse squared = adjacent squared + opposite squared) seems tempting to get the length of the diagonal between them, but around a curve this needs an extra step. Work with the height half way between the points (rather than the higher or lower point's altitude) as the height of the new spheroid surface. The resulting geodesic distance is the length of the adjacent, and the difference in height between the points is the opposite side.

It is also possible to work out the line of sight distance between the points, but not using geodesics. With each point, use lat_long_to_xyz/latLongToXYZ to get cartesian coordinates of the two points. Then use three dimensional Pythagoras theorem to get the distance (square root of XDifference squared plus YDifference squared plus ZDifference squared). Of course, this will not actually help to work out if the two points can actually be seen from each other around the curve of the World.

Great ellipse

It is easily possible to use the two geodesics methods to determine a great ellipse distance. If you are starting with two points, and you want a great ellipse that intersects both of them, use get_geodesic/getGeodesic to get the initial azimuth. Once you have a start point and azimuth, use get_geodesic_destination/getGeodesicDestination to get the result of a geodetic vector in that direction, of longer than half the distance around the ellipsoid. You can either specify the length 30000000 for most global ellipsoids, or for any oblate spheroid, a length of more than π times the major axis radius (over half way around the equator) and less than 4 times the major axis radius (the shortest distance across an infinitely thin ellipsoid and back) will work perfectly. Then use get_geodesic/getGeodesic to get from the resulting point back to the start of the ellipse. The initial length plus the final length add up to the length of the great ellipse.

$lat = 53;
$lon = -2;
$dir = 12;
$grutoolbox = Grid_Ref_Utils::toolbox();
$WGS84 = $grutoolbox->get_ellipsoid('WGS84');
$length1 = ( ( M_PI + 4 ) / 2 ) * max( $WGS84['a'], $WGS84['b'] );
$p2 = $grutoolbox->get_geodesic_destination( $lat, $lon, $length1, $dir, $WGS84, $grutoolbox->DATA_ARRAY );
$length2 = $grutoolbox->get_geodesic( $p2, $lat, $lon, $WGS84, $grutoolbox->DATA_ARRAY );
$greatellipselength = $length1 + $length2[0];

Changelogs

Changes in version 3.0, 23/12/2020

Changes in version 2.1.2 for PHP, 19/10/2020

Changes in version 2.1.1 for PHP, 16/9/2020

Changes in version 2.1, 9/5/2010

Changes in version 2.0 for PHP, 6/5/2010

Changes in version 1.1 and 1.2 for PHP

Changes in version 1.0 for PHP, 10/4/2010

API documentation

The API documentation has its own page.

This site was created by Mark "Tarquin" Wilton-Jones.
Don't click this link unless you want to be banned from our site.