Go Down

Topic: GPS Lon/Lat to British OS Grid (Read 444 times) previous topic - next topic

Dan0

has anyone attempted or had any success in converting GPS coordinates to British OS Grid?

here are some articles I found whilst researching

https://www.codeproject.com/Articles/13577/GPS-Deriving-British-Ordnance-Survey-Grid-Referenc

https://www.movable-type.co.uk/scripts/latlong-os-gridref.html



thanks

Dan

jremington

Not a task for the standard AVR-based Arduino, because it does not support the 64 bit double data type (required for accuracy).

srnet

It would be interesting to see how close a standard Arduino could get mind.

With a GPS capable of around 6M accuracy that is just about the limit of navigating and reading using the UKs 1:25000 maps.

Nope, I ma not offering to do the calc code.
http://www.50dollarsat.info/
http://www.loratracker.uk/

Dan0

thanks for the replies

Not a task for the standard AVR-based Arduino, because it does not support the 64 bit double data type (required for accuracy).
I guess thats why nobody has done it then :)

couka

Not a task for the standard AVR-based Arduino, because it does not support the 64 bit double data type (required for accuracy).
Nothing nickgammon's BigNumber library shouldn't be able to solve :)
youtube.com/DieBastler1234
Don't send me technical questions via PM. They will be deleted unanswered.

jonp

Yes it can be done pretty well: not with military precision but certainly within 10metres which is good enough to use on a hike.  The LLtoNE function below outputs OS eastings, OS northings each in 6 digits (the first from each defines the 100km square usually represented by a two letter code in the map key eg 'NZ' in the north east of England, the last 5 from each give a 10 digit grid ref), and the elevation based on the Airy Ellipsoid used for heights on OS maps.  Its float inputs must be the longitude, latitude and mean sea level elevation.  If you're using data from the GGA NMEA sentence from a cheap Ublox module for example, you'll need to change the format to pure fractions of degrees from their degree and minutes format (they output in format longitude dddmm.mmmm, latitude ddmm.mmmmm) and make them negative if west or south.  To get height above sea level, just add the 'antenna altitude' and 'geoidal separation' fields.  Doing all of this will use up about 8-10k, so quite a bit of your Arduino's program memory. 

I'm trying to pare it down at the moment to use on a tiny85.

I can't take credit for the maths - that comes from http://www.haroldstreet.org.uk/other/convertosgblatlong/



#include <avr/io.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>


//definitions for geometric conversions
#define deg2rad 0.017453292519943295 //(PI / 180)
#define rad2deg 57.29577951308232087 //(180/ PI)

#define a 6377563.396       // OSGB semi-major axis
#define b 6356256.91        // OSGB semi-minor axis
#define e0 400000           // OSGB easting of false origin
#define n0 -100000          // OSGB northing of false origin
#define f0 0.9996012717     // OSGB scale factor on central meridian
#define e2 0.0066705397616  // OSGB eccentricity squared
#define lam0 -0.034906585039886591  // OSGB false east
#define phi0 0.85521133347722145    // OSGB false north
#define af0   6375020.48098897069   //(a * f0)
#define bf0   6353722.49048791244   //(b * f0)
#define n 0.0016732202503250876 //(af0 - bf0) / (af0 + bf0)

#define WGS84_AXIS 6378137 // a
#define WGS84_ECCENTRIC 0.00669438037928458 //e
#define OSGB_AXIS 6377563.396 //a2
#define OSGB_ECCENTRIC 0.0066705397616   //e2
#define _xp -446.448  //OSGB/Airy datums/parameters
#define _yp 125.157
#define _zp -542.06
#define xrot -0.000000728190149026 //_xr -0.1502; (_xr / 3600) * deg2rad;
#define yrot -0.000001197489792340 //_yr -0.247; (_yr / 3600) * deg2rad;
#define zrot -0.000004082616008623 //_zr -0.8421; (_zr / 3600) * deg2rad;
#define _sf 0.0000204894  // s=20.4894 ppm


float Marc(float phi)   // used in LLtoNE function below
{
   float Marc = bf0 * (((1 + n + ((5 / 4) * (n * n)) + ((5 / 4) * (n * n * n))) * (phi - phi0))
      - (((3 * n) + (3 * (n * n)) + ((21 / 8) * (n * n * n))) * (sin(phi - phi0)) * (cos(phi + phi0)))
      + ((((15 / 8) * (n * n)) + ((15 / 8) * (n * n * n))) * (sin(2 * (phi - phi0))) * (cos(2 * (phi + phi0))))
      - (((35 / 24) * (n * n * n)) * (sin(3 * (phi - phi0))) * (cos(3 * (phi + phi0)))));
   return (Marc);
}

void LLtoNE(
   float latConv,  // in degrees
   float lonConv,   // in degrees
   float heightConv,  // in meters 
   uint32_t * const p_os_northings, //pointer to output variable for OS northings
   uint32_t * const p_os_eastings, //pointer to output variable for OS eastings
   signed int * const p_airy_elevation //pointer to output variable for elevation (based on Airy elipsoid used for heights OS maps)
   )
{
   latConv*= deg2rad;      // convert latitude to radians
   lonConv*= deg2rad;      // convert longitude to radians

   // Convert WGS84/GRS80 into OSGB36/Airy
   // convert to cartesian
   float v = WGS84_AXIS / (sqrt(1 - (WGS84_ECCENTRIC *(sin(latConv) * sin(latConv)))));
   float x = (v + heightConv) * cos(latConv) * cos(lonConv);
   float y = (v + heightConv) * cos(latConv) * sin(lonConv);
   float z = ((1 - WGS84_ECCENTRIC) * v + heightConv) * sin(latConv);
   // transform cartesian
   float hx = x + (x * _sf) - (y * zrot) + (z * yrot) + _xp;
   float hy = (x * zrot) + y + (y * _sf) - (z * xrot) + _yp;
   float hz = (-1 * x * yrot) + (y * xrot) + z + (z * _sf) + _zp;
   // Convert back to lat, lon
   lonConv = atan(hy / hx);
   float p = sqrt((hx * hx) + (hy * hy));
   latConv = atan(hz / (p * (1 - OSGB_ECCENTRIC)));
   v = OSGB_AXIS / (sqrt(1 - OSGB_ECCENTRIC * (sin(latConv) * sin(latConv))));
   float errvalue = 1.0;
   float lat1 = 0;
   while (errvalue > (1/1024))
   {
      lat1 = atan((hz + OSGB_ECCENTRIC * v * sin(latConv)) / p);
      errvalue = abs(lat1 - latConv);
      latConv = lat1;
   }
   *p_airy_elevation = p / cos(latConv) - v;



   // Convert OSGB36/Airy into OS grid eastings and northings
   // easting
   float slat2 = sin(latConv) * sin(latConv);
   float nu = af0 / (sqrt(1 - (e2 * (slat2))));
   float rho = (nu * (1 - e2)) / (1 - (e2 * slat2));
   float eta2 = (nu / rho) - 1;
   float pp = lonConv - lam0;
   float IV = nu * cos(latConv);
   float clat3 = pow(cos(latConv), 3);
   float tlat2 = tan(latConv) * tan(latConv);
   float V = (nu / 6) * clat3 * ((nu / rho) - tlat2);
   float clat5 = pow(cos(latConv), 5);
   float tlat4 = pow(tan(latConv), 4);
   float VI = (nu / 120) * clat5 * ((5 - (18 * tlat2)) + tlat4 + (14 * eta2) - (58 * tlat2 * eta2));
   *p_os_eastings = e0 + (pp * IV) + (pow(pp, 3) * V) + (pow(pp, 5) * VI);
   // northing
   float M = Marc(latConv);
   float I = M + (n0);
   float II = (nu / 2) * sin(latConv) * cos(latConv);
   float III = ((nu / 24) * sin(latConv) * pow(cos(latConv), 3)) * (5 - pow(tan(latConv), 2) + (9 * eta2));
   float IIIA = ((nu / 720) * sin(latConv) * clat5) * (61 - (58 * tlat2) + tlat4);
   *p_os_northings = I + ((pp * pp) * II) + (pow(pp, 4) * III) + (pow(pp, 6) * IIIA);
}



AWOL

...and that's why we use code tags

jonp

apologies for not putting in code tags - first post to any forum!

AWOL

And yet you managed to find the TTY tags.

srnet

#9
Jan 28, 2019, 03:12 pm Last Edit: Jan 28, 2019, 03:16 pm by srnet
I converted the code to pass variables rather than pointers and took a reference location from the Digital OS map.

//Reference location North East corner of Cardiff Castle Wall, taken from OS Memory Map
//Is UK Grid ST 18138 76715

const float TestLatitude = 51.48343;        
const float TestLongitude = -3.18030;
const float TestAltitude = 48;

The code calculates the grid ref as;

(3) 18137  (1) 76714

Which is within 1M of the location indicated by Latitude and longitude.

So the normal position error from a GPS derived location will predominate.


Well done !



http://www.50dollarsat.info/
http://www.loratracker.uk/

Dan0

interesting, will have a little go of that later

tasmod

What reference does it use to numerically indicate the letters?  From the ST posted above it shows as 3 and 1 to represent the letters in the 12 digit string.
Given that 5 letters cover the whole of GB for the first letter.  How does it handle the A to Z (minus I) of the second letter say?

srnet

Have a look at the grid square index of the UK and why ST is (3) (1) should be obvious.
http://www.50dollarsat.info/
http://www.loratracker.uk/

tasmod

Thanks, yes just found the OS description of the number as letter relationship.  Quite logical.

Never used that format always used the two letters but does make it easier in programming.

OV00 is the only 'odd' square in the UK and only exposed at low tide.

alanj

I have this working.

I use it with an ESP 8266 very accurate.

Happy to share.

Go Up