How to Convert from Latitude and Longitude to Ordnance Survey Grid References

I have often naively wondered how to "convert from latitude and longitude to Ordnance Survey grid references" and recently finally found out. Though this question sounds fairly reasonable, it actually hides a large degree of naivety. Not only is it not actually possible to make an exact conversion between the two, the question also fails to realise that 'latitude' and 'longitude' have many wildly differing yet very sensible definitions. This article aims to be a hand-wavy explanation of why everything can't just be nice and simple while also working towards answering my initial question.

This article contains many simplifications to keep things short. Though I hope this is not a misleading summary if you're interested in getting a deep understanding of what is going on I strongly recommend not reading this and instead referring to the excellently written (beginner friendly) source materials on which this article is based.

Latitude and Longitude

To start lets define what we mean by latitude and longitude. Imagine that the picture below represents the earth and the star is some point on the earth. The latitude of the point is the angle between it and the equator (the thick west-east line) and the longitude is the angle between it and the 'prime meridian' (the thick north-south line). Finally, there is also some measure of altitude which you might think of as being a height above sea level.

Unfortunately, life is not as simple as it first appears. The most obvious problem is how we define the prime meridian. Though we can easily define the equator as being the half-way point between the poles (which we can pin-point because the earth is spinning), there is no such luck for the prime meridian. Instead we have to choose an arbitrary point to be 0°. This is made especially difficult by the fact that points on the Earth's surface are in constant motion. As you can imagine, the precise definition of this point differs wildly depending on who you ask.

The other issue we have is that the Earth is not a perfect sphere, in fact, it is slightly elongated and distorted. Though an ellipsoid (a slightly elongated sphere) is a pretty good approximation, it leaves us with even more questions. How big is the ellipsoid? How distorted is it? Where is the centre? Where is the equator?

Why so many Ellipsoids?

There are many standard ellipsoids in use, and each has its own benefits and specific advantages. For example, on the (greatly exaggerated) cross section of the earth below there exists a globally best-fitting ellipsoid and a locally best-fitting one which is more accurate for a specific area of the Earth.

As a practical example, most GPS receivers use the WGS-84 (World Geodetic System, 1984) ellipsoid which is a good approximation to the shape of the Earth as a whole. Meanwhile, Ordnance Survey uses the Airy 1830 ellipsoid which fits the UK much better than WGS-84 but is a poor match for the rest of the world. Both are wholly sensible systems but each for different tasks and so standardising on a single, universal, ellipsoid doesn't really make sense.

So, with all that said, it should now be clear that when giving a position as a latitude and longitude it is very important to explicitly state what ellipsoid we're talking about and also what arbitrary assumptions we're making.

But wait... What about altitude?

So, we're making progress, but what do we mean by "height above average sea level"? When you're hundreds of miles in-land how should you compensate for the curvature of the earth when comparing your height (do we use our chosen ellipsoid)? Alternatively we could measure how strong the force of gravity is and assume that points with the same gravitational field strength are at the same altitude. Sounds good, but guess what? That shape is just as complex as the shape of the earth's surface and as a bonus, completely different.

If things couldn't get any worse, what on earth is average sea level anyway? If you compare the height of the sea around Cornwall in the south of England you'll find that the sea level is some 80cm below the global average - almost a meter!

Oh yes! This stuff is surprisingly tricky. For people who really care about altitude measurements there is a plethora of ellipsoids and more complex models of the Earth's shape, known as geoids, to play with. If you're like me, however, you probably don't care about having really accurate altitudes and so you can often just get away with using your original ellipsoid for measuring altitude and still get a semi-reasonable answers. Sometimes, at least. Phew!

Converting between Ellipsoids

So, given that GPS receivers and the Ordnance Survey use different ellipsoids, we need to be able to convert between latitudes and longitudes given in one ellipsoid to latitudes and longitudes in another. Unfortunately, all those arbitrary definitions we had to make about ellipsoids have introduced error in their measurement and definition. As such it means that it is not possible to convert from one coordinate system to another without incurring some loss of accuracy.

If you imagine two 1D coordinate systems 'A' and 'B', they will each assign different coordinates to the same physical locations ('red', 'green' and 'blue') due to both measurement error in the point's locations and differing models of the Earth's shape. As a result, there is no exact linear transformation between values in one coordinate system to another. The closest we can get still inevitably introduces some error.

Though more complex transformations can get closer to reality, they can never get all the way there. For the transformation between GPS' WGS-84 coordinates to those used by the Ordnance Survey, however, the linear transformation yields an error of no more than around 5 meters. Since this is better than most personal GPS receivers can manage anyway, a linear transformation is acceptable.

The conventional way to do this starts with converting the latitude, longitude and ellipsoidal height of a point (that is the distance the point is from the centre of the ellipsoid) to a 3D Cartesian (x,y,z) coordinate. This coordinate is then transformed using a standard transformation matrix (i.e. simple linear translate/scale/rotate). People who work with coordinate systems use the jargon term "Helmert transform" to describe this process. The new Cartesian coordinate is then mapped back to a latitude and longitude.

Latitude and Longitude to Grid References

OK, so we've converted from the latitude and longitude used by GPS receivers to the one used by Ordnance Survey (albeit with some loss of accuracy). The next, and thankfully final, task is to turn this into a grid reference.

Ordnance Survey grid references look like "SQ 21409 59110" and consist of three parts. The figure below (reproduced from A guide to coordinate systems in Great Britain) is called the 'National Grid' and shows how the letters correspond to an area of Great Britain. Each Square covers an area of 100km by 100km. The numbers then give the number of meters east and north (respectively) of the bottom-left corner of the named grid square.

The astute reader will point out that since the Earth is not flat and so these 'Eastings' and 'Northings', as they're known, don't make any sense. Since maps are flat (and not curved like the Earth or our ellipsoids), we must decide how to 'project' our points from these curvy surfaces to the flat world of the map. There are many projections out there but the Ordnance Survey have chosen a particular Transverse Mercator projection which causes relatively little distortion to most areas in Great Britain. This projection is exactly what is used to get the Eastings and Northings used in Ordnance Survey grid references.

Summary

So, to sum up, converting between "GPS and Ordnance Survey coordinates" is a pretty involved task filled with lots of surprising gotchas. For the most common case and when you don't care too much about coordinates incorrect by multiple meters, the process can be described as follows:

1. Convert the WGS-84 latitude, longitude and altitude into Cartesian coordinates.
2. Use a Helmert transform before converting back to an Airy 1830 latitude, longitude and altitude to account for various differences in the way these two systems are defined and measured. (Yes, this is the step with the most hand-waving.)
3. Project the points on the Airy 1830 ellipsoid onto a 2D plane using a Transverse Mercator projection to give us Eastings and Northings.
4. Re-write this as a Ordnance Survey grid references.

If this sounded pretty fun, I highly recommend reading the source behind this article. It is far more detailed and has very minimal hand-waving by comparison. It also contains numerous fascinating historical anecdotes which show just the incredible achievements of early mapping efforts and the even more incredible capabilities of satellite assisted technologies.

Sources

The more interested reader is enthusiastically pointed towards A guide to coordinate systems in Great Britain. This booklet an absolutely fascinating, approachable and remarkably thorough introduction to this topic and is the source for much of information in this article.

Appendix: Implementation in C

If you're interested in seeing a C implementation of this process or simply want a library which does it all for you, I have ported/re-written Chris Veness' implementation of the conversion process to C:

Git Repo: https://github.com/mossblaser/os_coord

The example program is reproduced below:

#include <stdio.h>
#include <stdlib.h>

#include "os_coord.h"
#include "os_coord_math.h"
#include "os_coord_data.h"
#include "os_coord_transform.h"
#include "os_coord_ordinance_survey.h"

int
main(int argc, char *argv[])
{
    if (argc != 4) {
        fprintf(stderr, "%s: Usage %s [latitude] [longitude] [ellipsoidal height]\n"
                      , argv[0]
                      , argv[0]
                      );
        return -1;
    }

    // Parse arguments
    double gps_latitude           = strtod(argv[1], NULL);
    double gps_longitude          = strtod(argv[2], NULL);
    double gps_ellipsoidal_height = strtod(argv[3], NULL);

    // Convert to radians
    os_lat_lon_t   home_ll_wgs84 = { .lat=DEG_2_RAD(gps_latitude)
                                   , .lon=DEG_2_RAD(gps_longitude)
                                   , .eh=gps_ellipsoidal_height
                                   };

    // Convert from the WGS-84 ellipsoid to 3D Cartesian coordinates
    os_cartesian_t home_c_wgs84 = os_lat_lon_to_cartesian(home_ll_wgs84, OS_EL_WGS84);

    // Perform the Helmert transform to (approximately) correct for differences in
    // datum
    os_cartesian_t home_c_airy30 = os_helmert_transform(home_c_wgs84, OS_HE_WGS84_TO_OSGB36);

    // Convert back to latitude, longitude and ellipsoidal height to give
    // coordinates on the Airy 1830 ellipsoid.
    os_lat_lon_t home_ll_airy30 = os_cartesian_to_lat_lon(home_c_airy30, OS_EL_AIRY_1830);

    // Project the coordinates to the National Grid's 2D projection
    os_eas_nor_t home_en_airy30 = os_lat_lon_to_tm_eas_nor(home_ll_airy30, OS_TM_NATIONAL_GRID);

    // Turn these Eastings and Northings into a grid references
    os_grid_ref_t home_grid_ref = os_eas_nor_to_grid_ref(home_en_airy30, OS_GR_NATIONAL_GRID);

    // Print the result
    if (home_grid_ref.code[0] != '\0') {
        printf("%s %05.0f %05.0f (Altitude: %0.1fm)\n", home_grid_ref.code, home_grid_ref.e, home_grid_ref.n, home_grid_ref.h);
        return 0;
    } else {
        fprintf(stderr, "%s: Coordinate not covered by National Grid\n", argv[0]);
        return -1;
    }
}