Click here to Skip to main content
15,860,943 members
Articles / Programming Languages / C#

GPS - Deriving British Ordnance Survey Grid Reference from NMEA data

Rate me:
Please Sign up or sign in to vote.
4.62/5 (23 votes)
23 May 2006CPOL7 min read 106.7K   3K   58   11
An article on a C# class to convert GPS derived NMEA data to the British Ordnance Survey Grid.

Sample Image - NMEA2OSGdemo-screenshot.jpg

Introduction

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.

Background

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:

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

it should say:

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

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

C#
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:

C#
public delegate void PDOPReceivedEventHandler(double value);

add the line:

C#
public delegate void EllipsoidHeightReceivedEventHandler(double value);

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

C#
public event PDOPReceivedEventHandler PDOPReceived;

add the line:

C#
public event EllipsoidHeightReceivedEventHandler EllipsoidHeightReceived;

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

Replace:

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

with:

C#
case "$GPGSA":
   return ParseGPGSA(sentence); 
case "$GPGGA":
   return ParseGPGGA(sentence);
default:
 // 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 http://www.commlinx.com.au/NMEA_sentences.htm. 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).

C#
public bool ParseGPGGA(string sentence)
{
     // Divide the sentence into words
     string[] Words = GetWords(sentence);
     if (Words[9] != "")
     {
         if (EllipsoidHeightReceived != null)
             EllipsoidHeightReceived(double.Parse(Words[9]));
     }
    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 (http://www.carabus.co.uk/), and who kindly gave me permission to make use of his work.

Thanks are also due to Nick Daniels (nickdaniels.com) 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.

C#
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.

EllipsoidDatumDescription
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.

C#
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.

C#
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:

C#
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);
  lat=lat0;
}
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.

C#
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.

C#
(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:

C#
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:

C#
//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.

C#
//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 NMEA2OSGdemo.zip.

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.

License

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


Written By
Software Developer University Of East Anglia, UK
United Kingdom United Kingdom
This member has not yet provided a Biography. Assume it's interesting and varied, and probably something to do with programming.

Comments and Discussions

 
QuestionGPGGA to Google Maps coordinates?! Pin
Mazen el Senih29-Jan-17 11:32
professionalMazen el Senih29-Jan-17 11:32 
QuestionConversion algorithm(s) - probably incomplete? Pin
richardw486-Jun-08 4:12
richardw486-Jun-08 4:12 
AnswerRe: Conversion algorithm(s) - probably incomplete? Pin
Alex@UEA6-Jun-08 4:44
Alex@UEA6-Jun-08 4:44 
GeneralRe: Conversion algorithm(s) - probably incomplete? Pin
richardw486-Jun-08 4:57
richardw486-Jun-08 4:57 
GeneralRe: Conversion algorithm(s) - probably incomplete? Pin
TMH Owen22-Sep-08 12:03
TMH Owen22-Sep-08 12:03 
GeneralRe: Conversion algorithm(s) - probably incomplete? Pin
phwip29-Mar-10 23:38
phwip29-Mar-10 23:38 
GeneralLittle remark Pin
Georg Kalus29-Jun-07 7:40
Georg Kalus29-Jun-07 7:40 
GeneralRe: Little remark Pin
Alex@UEA16-Jul-07 2:21
Alex@UEA16-Jul-07 2:21 
Glad you liked it, thank you for the tip on setting culture

Alex
GeneralOrdnance Survey Pin
Dave Cross23-May-06 7:09
professionalDave Cross23-May-06 7:09 
GeneralRe: Ordnance Survey Pin
Alex@UEA23-May-06 21:20
Alex@UEA23-May-06 21:20 
GeneralLooks good Pin
Ed.Poore27-Mar-06 6:30
Ed.Poore27-Mar-06 6: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.