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

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)
{
sentence = sentence.Substring(0, sentence.IndexOf("*"));
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.

Replace:

case "$GPGSA":
return ParseGPGSA(sentence);
default:
return false;

with:

case "$GPGSA":
return ParseGPGSA(sentence);
case "$GPGGA":
return ParseGPGGA(sentence);
default:
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)`

.

public bool ParseGPGGA(string sentence)
{
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.

public double deciLat;
public double deciLon;
public bool ParseNMEA(string Nlat, string Nlon, double height)
{
deciLat = Convert.ToDouble(Nlat.Substring(0, Nlat.IndexOf("Â°")));
deciLon = Convert.ToDouble(Nlon.Substring(0, Nlon.IndexOf("Â°")));
Nlat = Nlat.Substring(Nlat.IndexOf("Â°") + 1);
Nlon = Nlon.Substring(Nlon.IndexOf("Â°") + 1);
deciLat += (Convert.ToDouble(Nlat.Substring(0, Nlat.IndexOf("\""))))/60;
deciLon += (Convert.ToDouble(Nlon.Substring(0, Nlat.IndexOf("\""))))/60;
Nlat = Nlat.Substring(Nlat.IndexOf("\"") + 1);
Nlon = Nlon.Substring(Nlon.IndexOf("\"") + 1);
if (Nlat == "S") deciLat = 0 - deciLat;
if (Nlon == "W") deciLon = 0 - deciLon;
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.

Ellipsoid | Datum | Description |

Airy | OSGB36 | Used for the OS British National Grid |

GRS80 | WGS84 | A 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:

- Convert latitude and longitude to Cartesian coordinates (these also include height data, and have three parameters, X, Y, and Z).
- Transform to the new system by applying the 7 parameters and using a little math.
- 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 `lat`

itude and `lon`

gitude, `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; 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);
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.

private static double deg2rad = Math.PI / 180;
private static double rad2deg = 180.0 / Math.PI;
private bool Transform(double WGlat, double WGlon, double height)
{
double radWGlat = WGlat * deg2rad;
double radWGlon = WGlon * deg2rad;
double a = 6378137; double e = 0.00669438037928458; double h = height; double a2 = 6377563.396; double e2 = 0.0066705397616; 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;
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);
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;
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;
}
newLat = newLat * rad2deg;
newLon = newLon * rad2deg;
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:

private void LLtoNE(double lat, double lon)
{
double phi = lat * deg2rad; double lam = lon * deg2rad; double a = 6377563.396; double b = 6356256.91; double e0 = 400000; double n0 = -100000; double f0 = 0.9996012717; double e2 = 0.0066705397616; double lam0 = -0.034906585039886591; double phi0 = 0.85521133347722145; double af0 = a * f0;
double bf0 = b * f0;
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);
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);
east = Math.Round(east); north = Math.Round(north);
if (NorthingEastingReceived != null) NorthingEastingReceived(north, east);
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.

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; if (tmp > 7.5) tmp = tmp + 1;
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;
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.