Click here to Skip to main content
13,257,276 members (47,657 online)
Click here to Skip to main content
Add your own
alternative version


56 bookmarked
Posted 26 Mar 2006

GPS - Deriving British Ordnance Survey Grid Reference from NMEA data

, 23 May 2006
Rate this:
Please Sign up or sign in to vote.
An article on a C# class to convert GPS derived NMEA data to the British Ordnance Survey Grid.

Sample Image - NMEA2OSGdemo-screenshot.jpg


The NMEA data exported by GPS units gives latitude and longitude and geometric distance above the WGS84 (GRS80) reference ellipsoid. The Ordnance Survey maps for Great Britain use grid references based on the Airy Spheroid (OSGB36) reference ellipsoid.

The following article discusses a C# class to convert GPS derived NMEA data to the British Ordnance Survey Grid.


This article assumes that the NMEA data is being derived using Jon Person's excellent NMEA interpreter class, as discussed in Writing Your Own GPS Applications: Part 2.

Using the code

1. Modify the NMEA interpreter class

In order to perform the transformation, some additional data is required in the form of ellipsoid height, which can be derived from the NMEA $GPGGA sentence. The $GPGGA sentence is not parsed in the original form of the NMEAinterpreter, so this must be added.

First off, we need to make a couple of small alterations to the NMEAinterpreter class; for an explanation of these changes, read the "Slight Error" comment.

In public bool ParseGPGSV(string sentence), where it says:

SignalToNoiseRatio = Convert.ToInt32(Words[Count * 4 + 2]);

it should say:

SignalToNoiseRatio = Convert.ToInt32(Words[Count * 4 + 3]);

Also, public string[] GetWords(string sentence) should be modified as follows:

public string[] GetWords(string sentence)
    //remove the * and checksum from the end of the sentence
    sentence = sentence.Substring(0, sentence.IndexOf("*")); 
    //now split it up
    return sentence.Split(',');

Now we need to add the parsing of the $GPGGA sentence. First, we will need to add a new delegate for the ellipsoid height event:

Below the line:

public delegate void PDOPReceivedEventHandler(double value);

add the line:

public delegate void EllipsoidHeightReceivedEventHandler(double value);

We also need to add a new event. Below the line:

public event PDOPReceivedEventHandler PDOPReceived;

add the line:

public event EllipsoidHeightReceivedEventHandler EllipsoidHeightReceived;

The public bool Parse(string sentence) routine needs a new case added for the $GPGGA sentence.


case "$GPGSA":
   return ParseGPGSA(sentence);      
 // Indicate that the sentence was not recognised
   return false;


case "$GPGSA":
   return ParseGPGSA(sentence); 
case "$GPGGA":
   return ParseGPGGA(sentence);
 // Indicate that the sentence was not recognised
  return false;

Now we need to add the parsing of the $GPGGA sentence.

There are numerous web pages which detail the structure of the various NMEA sentences. I got my information from If the first word of the GPGGA sentence is word 0, then word 9 contains the ellipsoid height.

public bool ParseGPGGA(string sentence) could be added below public bool ParseGPGSA(string sentence).

public bool ParseGPGGA(string sentence)
     // Divide the sentence into words
     string[] Words = GetWords(sentence);
     if (Words[9] != "")
         if (EllipsoidHeightReceived != null)
    return true;

So now, we have modified Jon Person's NMEAinterpreter class to give us the ellipsoid height. (In practice, you may want to pull out various other bits of data too.)

2. Create the NMEA2OSG class

The NMEA2OSG class will work in a similar way to the NMEAinterpreter class. Latitude, longitude, and ellipsoidal height will be passed to the NMEA2OSG class, which will then generate events producing the transformed data.

At this point, I have to give thanks and credit to Roger Muggleton, the creator of the pages which detailed the necessary calculations (, and who kindly gave me permission to make use of his work.

Thanks are also due to Nick Daniels ( for his invaluable advice on C# coding.

To quote Roger Muggleton, 'This is quite involved math, and I am not going to try to explain all of it. It is the same as that used in the Ordnance Survey document, A Guide to Coordinate Systems in Great Britain, and converts latitude and longitude to a Transverse Mercator projection.'

The NMEAinterpreter class produces latitude and longitude in the following form:

52°09.1461"N         002°33.3717"W

In order to work mathematically with the latitude and longitude, they need to be converted into decimal values. So, the first function of our class is to parse the latitude and longitude and then pass them, along with the ellipsoid height, to the transformation function.

We want to convert the minutes to degrees (divide by 60), and take account of the fact that Southerly and Westerly values are considered to be negative.

public double deciLat;
public double deciLon;

public bool ParseNMEA(string Nlat, string Nlon, double height)
    //grab the bit up to the °
    deciLat = Convert.ToDouble(Nlat.Substring(0, Nlat.IndexOf("°")));
    deciLon = Convert.ToDouble(Nlon.Substring(0, Nlon.IndexOf("°")));

    //remove that bit from the string now we've used it and the ° symbol itself 
    Nlat = Nlat.Substring(Nlat.IndexOf("°") + 1);
    Nlon = Nlon.Substring(Nlon.IndexOf("°") + 1);

    //grab the bit up to the ", divide it by 60 to convert to degrees 
    //and add it to our double value
    deciLat += (Convert.ToDouble(Nlat.Substring(0, Nlat.IndexOf("\""))))/60;
    deciLon += (Convert.ToDouble(Nlon.Substring(0, Nlat.IndexOf("\""))))/60;

    //remove that part of the string now and just leave the compass direction
    Nlat = Nlat.Substring(Nlat.IndexOf("\"") + 1);
    Nlon = Nlon.Substring(Nlon.IndexOf("\"") + 1);

    // check for negative directions
    if (Nlat == "S") deciLat = 0 - deciLat;
    if (Nlon == "W") deciLon = 0 - deciLon;

    //now we can transform them
    return Transform(deciLat, deciLon, height);

3. Coping with different ellipsoids

Map makers use an ellipsoid whose surface closely matches the surface of the area to be mapped. In order to define how this ellipsoid and the real Earth are aligned, they use a concept called a Datum (or Terrestrial Reference Frame, TRF). This Datum ties actual points on the surface of the Earth to latitude and longitude points on the reference ellipsoid. Since Earth and Datum do not match exactly, this tie up distorts the math a little, and makes it almost impossible to convert one system to another. Conversions are worthwhile, as they may improve accuracy by one hundred meters or more, but they are not exact.

AiryOSGB36Used for the OS British National Grid
GRS80WGS84A system used to describe the whole planet, NMEA data is referenced to WGS84

There are many ways to convert data from one system to another, the most accurate being the most complex! Here, we will use a 7 parameter Helmert transformation.

The process is in three parts:

  1. Convert latitude and longitude to Cartesian coordinates (these also include height data, and have three parameters, X, Y, and Z).
  2. Transform to the new system by applying the 7 parameters and using a little math.
  3. Convert back to latitude and longitude.

We want to transform a GRS80 location to Airy, e.g., a GPS reading to the Airy Spheroid.

First, we need to convert our measurements from degrees to radians. This is accomplished by multiplying our latitude or longitude by PI / 180. This conversion factor is called deg2rad.

The following code converts latitude and longitude to Cartesian coordinates. The input parameters are: WGS84 latitude and longitude, axis is the GRS80/WGS84 major axis in meters, ecc is the eccentricity, and height is the height above the ellipsoid.

v = axis / (Math.Sqrt (1 - ecc * (Math.Pow (Math.Sin(lat), 2))));
x = (v + height) * Math.Cos(lat) * Math.Cos(lon);
y = (v + height) * Math.Cos(lat) * Math.Sin(lon);
z = ((1 - ecc) * v + height) * Math.Sin(lat);

The transformation requires the 7 parameters: xp, yp, and zp correct the coordinate origin, xr, yr, and zr correct the orientation of the axes, and sf deals with the changing scale factors.

sf = s * 0.000001;  // scale factor is in parts per million
xrot = (xr / 3600) * deg2rad;
yrot = (yr / 3600) * deg2rad;
zrot = (zr / 3600) * deg2rad;

hx = x + (x * sf) - (y * zrot) + (z * yrot) + xp;
hy = (x * zrot) + y + (y * sf) - (z * xrot) + yp;
hz = (-1 * x * yrot) + (y * xrot) + z + (z * sf) + zp;

Then, we convert back to latitude and longitude, where axis and ecc are those of the Airy Spheroid:

lon = Math.Atan(y / x);
p = Math.Sqrt((x * x) + (y * y));
lat = Math.Atan(z / (p * (1 - ecc)));
v = axis / (Math.Sqrt(1 - ecc * (Math.Sin(lat) * Math.Sin(lat))));
errvalue = 1.0;
lat0 = 0;
while (errvalue > 0.001)
  lat0 = Math.Atan((z + ecc * v * Math.Sin(lat)) / p);  
  errvalue = Math.Abs(lat0 - lat);
height = p / Math.Cos(lat) - v;

The above involves some iteration, hence the while loop. In this case, ecc is the eccentricity of the Airy Spheroid, and axis is its major axis.

The Transform method below will convert the WGS84 based latitude and longitude to the OSGB36 latitude and longitude; it then calls a method (LLtoNE) which will convert these to Ordinance Survey Eastings and Northings.

private static double deg2rad = Math.PI / 180;
private static double rad2deg = 180.0 / Math.PI;

private bool Transform(double WGlat, double WGlon, double height)
     //first off convert to radians
     double radWGlat = WGlat * deg2rad;
     double radWGlon = WGlon * deg2rad;
          //these are the values for WGS86(GRS80) to OSGB36(Airy)
     double a = 6378137;              // WGS84_AXIS
     double e = 0.00669438037928458;  // WGS84_ECCENTRIC
     double h = height;               // height above datum (from $GPGGA sentence)
     double a2 = 6377563.396;         // OSGB_AXIS
     double e2 = 0.0066705397616;     // OSGB_ECCENTRIC 
     double xp = -446.448;
     double yp = 125.157;
     double zp = -542.06;
     double xr = -0.1502;
     double yr = -0.247;
     double zr = -0.8421;
     double s = 20.4894;

     // convert to cartesian; lat, lon are in radians
     double sf = s * 0.000001;
     double v = a / (Math.Sqrt(1 - (e * (Math.Sin(radWGlat) * Math.Sin(radWGlat)))));
     double x = (v + h) * Math.Cos(radWGlat) * Math.Cos(radWGlon);
     double y = (v + h) * Math.Cos(radWGlat) * Math.Sin(radWGlon);
     double z = ((1 - e) * v + h) * Math.Sin(radWGlat);
     // transform cartesian
     double xrot = (xr / 3600) * deg2rad;
     double yrot = (yr / 3600) * deg2rad;
     double zrot = (zr / 3600) * deg2rad;
     double hx = x + (x * sf) - (y * zrot) + (z * yrot) + xp;
     double hy = (x * zrot) + y + (y * sf) - (z * xrot) + yp;
     double hz = (-1 * x * yrot) + (y * xrot) + z + (z * sf) + zp;
     // Convert back to lat, lon
     double newLon = Math.Atan(hy / hx);
     double p = Math.Sqrt((hx * hx) + (hy * hy));
     double newLat = Math.Atan(hz / (p * (1 - e2)));
     v = a2 / (Math.Sqrt(1 - e2 * (Math.Sin(newLat) * Math.Sin(newLat))));
     double errvalue = 1.0;
     double lat0 = 0;
     while (errvalue > 0.001)
         lat0 = Math.Atan((hz + e2 * v * Math.Sin(newLat)) / p);
         errvalue = Math.Abs(lat0 - newLat);
         newLat = lat0;
     //convert back to degrees
     newLat = newLat * rad2deg;
     newLon = newLon * rad2deg;
     //convert lat and lon (OSGB36) to OS 6 figure northing and easting
     LLtoNE(newLat, newLon);
     return true;

4. Convert latitude and longitude to a British National Grid Reference

We have already converted our measurements of latitude and longitude so that they are now based on the Airy Spheroid.

The Airy Spheroid was defined in 1830 by Britain's Astronomer Royal, Sir George Airy. This has a major axis (distance between the poles) of 6377563.396 meters, and a minor axis of 6356256.910 meters. In the equations that follow, we use a constant called the eccentricity, this is simply the difference between the squares of the axes, divided by the square of the major axis.

(major*major - minor*minor) / (major*major);

Using this formula with the Airy axes, we get an eccentricity of 0.006670539761597337.

Other values in this method are taken from the aforementioned Ordnance Survey publication: A Guide to Coordinate Systems in Great Britain.

This method converts the latitude and longitude to a 6 figure Northing and a 6 figure Easting. An Event is then raised containing these values.

At the end of the method, we call another method (NE2NGR) - this converts the 6 figure Northing and the 6 figure Easting produced to the National Grid Reference Convention.

LLtoNE also calls the following Marc' mathematical method:

private double Marc(double bf0, double n, double phi0, double phi)
     return bf0 * (((1 + n + ((5 / 4) * (n * n)) + ((5 / 4) * (n * n * n)))_
                * (phi - phi0))
        - (((3 * n) + (3 * (n * n)) + ((21 / 8) * (n * n * n))) *_
                (Math.Sin(phi - phi0)) * (Math.Cos(phi + phi0)))
        + ((((15 / 8) * (n * n)) + ((15 / 8) * (n * n * n))) *_
                (Math.Sin(2 * (phi - phi0))) * (Math.Cos(2 * (phi + phi0))))
        - (((35 / 24) * (n * n * n)) * (Math.Sin(3 * (phi - phi0))) *_<
                (Math.Cos(3 * (phi + phi0)))));

Here is the LLtoNE method:

//converts lat and lon (OSGB36) to OS 6 figure northing and easting
private void LLtoNE(double lat, double lon)
     double phi = lat * deg2rad;          // convert latitude to radians
     double lam = lon * deg2rad;          // convert longitude to radians
     double a = 6377563.396;              // OSGB semi-major axis
     double b = 6356256.91;               // OSGB semi-minor axis
     double e0 = 400000;                  // easting of false origin
     double n0 = -100000;                 // northing of false origin
     double f0 = 0.9996012717;            // OSGB scale factor on central meridian
     double e2 = 0.0066705397616;         // OSGB eccentricity squared
     double lam0 = -0.034906585039886591; // OSGB false east
     double phi0 = 0.85521133347722145;   // OSGB false north
     double af0 = a * f0;
     double bf0 = b * f0;
     // easting
     double slat2 = Math.Sin(phi) * Math.Sin(phi);
     double nu = af0 / (Math.Sqrt(1 - (e2 * (slat2))));
     double rho = (nu * (1 - e2)) / (1 - (e2 * slat2));
     double eta2 = (nu / rho) - 1;
     double p = lam - lam0;
     double IV = nu * Math.Cos(phi);
     double clat3 = Math.Pow(Math.Cos(phi), 3);
     double tlat2 = Math.Tan(phi) * Math.Tan(phi);
     double V = (nu / 6) * clat3 * ((nu / rho) - tlat2);
     double clat5 = Math.Pow(Math.Cos(phi), 5);
     double tlat4 = Math.Pow(Math.Tan(phi), 4);
     double VI = (nu / 120) * clat5 * ((5 - (18 * tlat2)) + tlat4 + (14 * eta2)_
         - (58 * tlat2 * eta2));
     double east = e0 + (p * IV) + (Math.Pow(p, 3) * V) + (Math.Pow(p, 5) * VI);
     // northing
     double n = (af0 - bf0) / (af0 + bf0);
     double M = Marc(bf0, n, phi0, phi);
     double I = M + (n0);
     double II = (nu / 2) * Math.Sin(phi) * Math.Cos(phi);
     double III = ((nu / 24) * Math.Sin(phi) * Math.Pow(Math.Cos(phi), 3))_
         * (5 - Math.Pow(Math.Tan(phi), 2) + (9 * eta2));
     double IIIA = ((nu / 720) * Math.Sin(phi) * clat5) * (61 - (58 * tlat2)_
         + tlat4);
     double north = I + ((p * p) * II) + (Math.Pow(p, 4) * III)_
         + (Math.Pow(p, 6) * IIIA);
     // make whole number values
     east = Math.Round(east);   // round to whole number
     north = Math.Round(north); // round to whole number
      // Notify the calling application of the change
     if (NorthingEastingReceived != null) NorthingEastingReceived(north, east);
     // convert to nat grid ref
     NE2NGR(east, north);

To reduce the number of figures needed to give a National Grid reference, the grid is divided into 100 km squares which each have a two-letter code. National Grid positions can be given with this code, followed by an Easting between 0 and 100 000 m and a Northing between 0 and 100 000 m.

NE2NGR takes Northing and Easting and raises an event with the National Grid Reference.

//convert 12 (6e & 6n) figure numeric to letter and number grid system
private void NE2NGR(double east, double north)
  double eX = east / 500000;
  double nX = north / 500000;
  double tmp = Math.Floor(eX) - 5.0 * Math.Floor(nX) + 17.0;
  nX = 5 * (nX - Math.Floor(nX));
  eX = 20 - 5.0 * Math.Floor(nX) + Math.Floor(5.0 * (eX - Math.Floor(eX)));
  if (eX > 7.5) eX = eX + 1;    // I is not used
  if (tmp > 7.5) tmp = tmp + 1; // I is not used
  string eing = Convert.ToString(east);
  string ning = Convert.ToString(north);
  int lnth = eing.Length;
  eing = eing.Substring(lnth - 5);
  lnth = ning.Length;
  ning = ning.Substring(lnth - 5);
  string ngr = Convert.ToString((char)(tmp + 65)) + Convert.ToString((char)_
     (eX + 65)) + " " + eing + " " + ning;
  // Notify the calling application of the change
  if (NatGridRefReceived != null) NatGridRefReceived(ngr);

And that's it!

5. Putting it to use

So now, we have created the class, there is a small demo project in C#, .NET 2.0 demonstrating its use by reading the NMEA from a GPS connected to a serial port and displaying the OS grid reference. It's available for download as

The project makes use of the .NET 2.0 SerialPort class to read the GPS. The NMEA strings are parsed via NMEAinterpreter. The lat and lon are parsed via NMEA2OSG. The resulting events cause the respective text boxes to be updated with the appropriate data.

I hope this article has been useful and/or informative.

Accuracy of the conversion should be within less than 4 meters (which is better than the inbuilt conversion of my Garmin Vista).

If you need to know more about the underlying mathematics, I suggest you go to the resources referred to in the article, as a deeper explanation is beyond me!

Points of interest

Oh, and I like to know the number of satellites in view and how many are in use by the GPS, so this functionality is also added to the NMEAinterpreter class.


This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)


About the Author

Software Developer University Of East Anglia, UK
United Kingdom United Kingdom
No Biography provided

You may also be interested in...

Comments and Discussions

QuestionGPGGA to Google Maps coordinates?! Pin
Mazen el Senih29-Jan-17 12:32
professionalMazen el Senih29-Jan-17 12:32 
QuestionConversion algorithm(s) - probably incomplete? Pin
RVW6-Jun-08 5:12
memberRVW6-Jun-08 5:12 
I wish I'd come across this excellent article sooner as it probably would have saved me a considerable amount of work. However, I didn't.

I've been working recently on a project which required the conversion from WGS84 Latitude and Longitude to OSGB36 Eastings and Northings and back again. During the analysis I did quite a bit of research on the conversion algorithms and, as a result, I'm not completely convinced that the methods described here provide the full picture.

In section 6.3 of the Ordnance Survey Guide to Coordinate Transformations it refers to the: National Grid Transformation OSTN02 (ETRS89-OSGB36) (ETRS89 by the way is Europe's version of WGS84).
It reads:
<i><b>To cope with the distortions in the OSGB36 TRF (transformation reference frame), different transformations are needed in different parts of the country. For this reason, the national standard datum transformation between OSGB36 and ETRS89 is not a simple Helmert datum transformation.</b> Instead, Ordnance Survey has developed a ‘rubber-sheet’ transformation which works with a transformation grid expressed in easting and northing coordinates. The grids of northing and easting shifts between ETRS89 and OSGB36 cover Britain with a grid resolution of one kilometre. From these grids a northing and easting shift for each point to be transformed is obtained by a bi-linear interpolation. This is called the National Grid Transformation OSTN02, and it is freely available in software packages from the Ordnance Survey GPS Website,</i>

Your article doesn't appear to mention the OSTN02 'rubber sheet' transformation that's also required in the calculation. So whilst you can convert Lat/Long to Easting/Northing using the methods described in this article I would advise readers not to rely on the values obtained if positional accuracy is very critical.

Nevertheless, as I said at the start it's an excellent article and I am pleased you took the time and effort to write it and, more importantly, share it.
AnswerRe: Conversion algorithm(s) - probably incomplete? Pin
Alex@UEA6-Jun-08 5:44
memberAlex@UEA6-Jun-08 5:44 
GeneralRe: Conversion algorithm(s) - probably incomplete? Pin
RVW6-Jun-08 5:57
memberRVW6-Jun-08 5:57 
GeneralRe: Conversion algorithm(s) - probably incomplete? Pin
Tom Owen (Quest)22-Sep-08 13:03
memberTom Owen (Quest)22-Sep-08 13:03 
GeneralRe: Conversion algorithm(s) - probably incomplete? Pin
phwip30-Mar-10 0:38
memberphwip30-Mar-10 0:38 
GeneralLittle remark Pin
Georg Kalus29-Jun-07 8:40
memberGeorg Kalus29-Jun-07 8:40 
GeneralRe: Little remark Pin
Alex@UEA16-Jul-07 3:21
memberAlex@UEA16-Jul-07 3:21 
GeneralOrdnance Survey Pin
Dave Cross23-May-06 8:09
memberDave Cross23-May-06 8:09 
GeneralRe: Ordnance Survey Pin
Alex@UEA23-May-06 22:20
memberAlex@UEA23-May-06 22:20 
GeneralLooks good Pin
Ed.Poore27-Mar-06 7:30
memberEd.Poore27-Mar-06 7:30 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Praise Praise    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.

Permalink | Advertise | Privacy | Terms of Use | Mobile
Web01 | 2.8.171114.1 | Last Updated 24 May 2006
Article Copyright 2006 by Alex@UEA
Everything else Copyright © CodeProject, 1999-2017
Layout: fixed | fluid