Trying to push the limits of App Inventor! Snippets and Tutorials from Pura Vida Apps by Taifun.
The code was in a different programming language, C if I recall, a little beyond my comprehension.
Although there is a lot of involved maths, this is not a difficult task. There are some problems, however, firstly due to the different ways latitude and longitude can be defined. All mapping systems have to use a model of the shape of planet Earth, which is not a perfect sphere. The distance between the poles is less that the distance across the equator, making a shape called an ellipsoid. But the ellipsoid is not a perfect fit either.
Most accurate maps show only part of the surface of the Earth, usually a single country or part of a country. So the usual way of getting round the problem is to define an ellipsoid whose surface has a good alignment with reality over the area the maps is to cover. The Ordnance Survey maps of the United Kingdom use three such reference ellipsoids, one for the Island of Great Britain (England, Scotland and Wales), another for the Island of Ireland, and finally a third ellipsoid to deal with the Channel Islands. This article primarily describes Great Britain, but the same maths can be applied to most of the other two areas.
The traditional ellipsoid adopted for Great Britain is called the Airy Spheroid. This was defined in 1830 by Britain's Astronomer Royal, Sir George Airy. This has a semi-major axis of 6377563.396 metres, and a semi-minor axis of 6356256.910 metres. In the equations that follow we use a constant called the eccentricity, this is simply the diffence between the squares of the axes, divided by the square of the major axis.
function ecc(major, minor)
{
return (major*major - minor*minor) / (major*major);
}
Using the formula above with the Airy axes we get an eccentricity of 0.006670539761597337
Let us assume that that our measurements of latitude and longitude are based on the Airy spheroid. This may not be true, but I'll cover that later on. Firstly we need to convert our measurements from degrees to radians. There is accomplished by multiplying our latitude or longitude by PI / 180. I shall call this conversion factor deg2rad. A few other constants are also defined:
function LLtoNE(lat, lon)
{
var deg2rad = Math.PI / 180;
var rad2deg = 180.0 / Math.PI;
var phi = lat * deg2rad; // convert latitude to radians
var lam = lon * deg2rad; // convert longitude to radians
a = 6377563.396; // OSGB semi-major axis
b = 6356256.91; // OSGB semi-minor axis
e0 = 400000; // OSGB easting of false origin
n0 = -100000; // OSGB northing of false origin
f0 = 0.9996012717; // OSGB scale factor on central meridian
e2 = 0.0066705397616; // OSGB eccentricity squared
lam0 = -0.034906585039886591; // OSGB false east
phi0 = 0.85521133347722145; // OSGB false north
var af0 = a * f0;
var bf0 = b * f0;
// easting
var slat2 = Math.sin(phi) * Math.sin(phi);
var nu = af0 / (Math.sqrt(1 - (e2 * (slat2))));
var rho = (nu * (1 - e2)) / (1 - (e2 * slat2));
var eta2 = (nu / rho) - 1;
var p = lam - lam0;
var IV = nu * Math.cos(phi);
var clat3 = Math.pow(Math.cos(phi),3);
var tlat2 = Math.tan(phi) * Math.tan(phi);
var V = (nu / 6) * clat3 * ((nu / rho) - tlat2);
var clat5 = Math.pow(Math.cos(phi), 5);
var tlat4 = Math.pow(Math.tan(phi), 4);
var VI = (nu / 120) * clat5 * ((5 - (18 * tlat2)) + tlat4 + (14 * eta2) - (58 * tlat2 * eta2));
east = e0 + (p * IV) + (Math.pow(p, 3) * V) + (Math.pow(p, 5) * VI);
// northing
var n = (af0 - bf0) / (af0 + bf0);
var M = Marc(bf0, n, phi0, phi);
var I = M + (n0);
var II = (nu / 2) * Math.sin(phi) * Math.cos(phi);
var III = ((nu / 24) * Math.sin(phi) * Math.pow(Math.cos(phi), 3)) * (5 - Math.pow(Math.tan(phi), 2) + (9 * eta2));
var IIIA = ((nu / 720) * Math.sin(phi) * clat5) * (61 - (58 * tlat2) + tlat4);
north = I + ((p * p) * II) + (Math.pow(p, 4) * III) + (Math.pow(p, 6) * IIIA);
east = Math.round(east); // round to whole number
north = Math.round(north); // round to whole number
nstr = String(north); // convert to string
estr = String(east); // ditto
}
function Marc(bf0, n, phi0, phi)
{
var 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))) * (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)))));
return(Marc);
}
This is quite involved maths, and I am not going to try to explain it! This 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. These projections are often used for paper (flat) maps. You can find some explanation in the OS document.
The Javascript maths functions are fairly obvious, I hope! You will have to scroll far to the right to see some of the code. Some of the expressions are quite long-winded, so I have split them up a little to avoid too many confusing brackets and Math functions.
We can now try this with some test figures, latitude 52° 39’ 27.2531" north, 1° 43’ 4.5177" east. These figures are taken from a worked example in the Ordnance Survey document.
lat = 52.65757; // decimal version of latitude
lon = 1.71791; // decimal version of longitude
LLtoNE(lat, lon);
document.write("And the results are:");
document.write("Northings: ",nstr, " metres.");
document.write("Eastings: ", estr, " metres.");
And the results are:
Northings: 313177 metres.
Eastings: 651409 metres.
The NGR northings and eastings are shown along the axes of any Ordnance Survey map. Normally the NGR is shown as letters and numbers. The letters replace the 100km digit(s) of the eastings and northings. This could easily be determined from a simple look-up table, which would need to handle the ninety one 100km 'squares'. So here is a weird formula, source unknown, which does the job in less space:
function NE2NGR(east, north)
{
var eX = east / 500000;
var nX = north / 500000;
var 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;
var eing = String(east);
var ning = String(north);
var lnth = eing.length;
eing = eing.substring(lnth - 5, lnth);
lnth = ning.length;
ning = ning.substring(lnth - 5, lnth);
ngr = String.fromCharCode(tmp + 65) + String.fromCharCode(eX + 65) + " " + eing + " " + ning;
return ngr;
}
So for a worked example using the same latitude and longitude, we run the code:
lat = 52.65757;
lon = 1.71791;
LLtoNE(lat, lon);
document.write("And the results are:");
document.write("Northings: ",nstr, " metres.");
document.write("Eastings: ", estr, " metres.");
ngr = NE2NGR(east, north);
document.write("NGR:", ngr);
And the results are:
Northings: 313177 metres.
Eastings: 651409 metres.
NGR: TG 51409 13177
This is a one metre grid reference, it can be reduced to the standard 100m reference by omitting the last two figures of each five figure group: TG 514 131.
This seems to be the best algorithm and instruction I could find