GPS- Deriving British Ordnance Survey Grid Reference from NMEA data

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 the Earth and Datum do not match exactly, this tie up distorts the maths a little, and makes it almost impossible to convert one system to another. Conversions are worthwhile, as they may improve accuracy by one hundred metres 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:

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

Firstly 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 metres, 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);
      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 WGS84 based latitude and longitude to OSGB36 latitude and longitude it then calls a method ( LLtoNE) which will convert these to Ordnance 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;
    }

Now that we have corrected values, we can transform them to the OS grid- see the next page.