GPS- Deriving British Ordnance Survey Grid Reference from NMEA data

Part 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 metres, and a minor axis of 6356256.910 metres. 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 .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' mathmatical 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! On the final page I talk about a small demo project (C# 2005).

You might also like...

Comments

Contribute

Why not write for us? Or you could submit an event or a user group in your area. Alternatively just tell us what you think!

Our tools

We've got automatic conversion tools to convert C# to VB.NET, VB.NET to C#. Also you can compress javascript and compress css and generate sql connection strings.

“C++: an octopus made by nailing extra legs onto a dog.” - Steve Taylor