Dan0
January 22, 2019, 10:16am
1
1 Like
Not a task for the standard AVR-based Arduino, because it does not support the 64 bit double data type (required for accuracy).
srnet
January 22, 2019, 6:23pm
3
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.
Dan0
January 25, 2019, 8:48am
4
thanks for the replies
jremington:
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
January 25, 2019, 8:52am
5
jremington:
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
jonp
January 28, 2019, 10:14am
6
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);
}
system
January 28, 2019, 10:18am
7
...and that's why we use code tags
jonp
January 28, 2019, 10:36am
8
apologies for not putting in code tags - first post to any forum!
system
January 28, 2019, 1:17pm
9
And yet you managed to find the TTY tags.
srnet
January 28, 2019, 2:12pm
10
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 !
Dan0
February 11, 2019, 1:29pm
11
interesting, will have a little go of that later
tasmod
February 11, 2019, 5:12pm
12
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
February 11, 2019, 5:21pm
13
Have a look at the grid square index of the UK and why ST is (3) (1) should be obvious.
tasmod
February 11, 2019, 5:35pm
14
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
May 11, 2019, 7:56pm
15
I have this working.
I use it with an ESP 8266 very accurate.
Happy to share.