UKC

The questions get tougher...

Assuming that I now know two 10 figure grid references on the same OSGB sheet how do I calculate the bearing from one to the other (without a map and a compass)?

Some function of sine cosine or tangent? But what happens for due N S E or W where there's no right triangle? And how do you account for each quarter of the compass?

Any calculating climbers help out with this one?

In reply to The Crow: Why do you want to do it? is it for a spreadsheet formulae or are you just testing us? You could workout which quadrant and eliminate the zero's and infinities at due north etc with some "if" statements, cotan or inv-tan of opposite over adjacent(difference in long. and lat.) to work out the basic angle and add the appropriate amount for which quadrant from the "if's". There's sure to be some neat way with complex numbers, but I can't remember any of that stuff from school, sorry!

In reply to The Crow:

Well for any triangle the following hold true where H = Hypotenuse, A = Adjacent side and O = opposite side. Call the angle Theta.

Sin(Theta) = Opposite/ Hypotenuse
Cos(Theta) = Adjacent / Hypotenuse
Tan(Theta) = Opposite / Adjacent

In the latter case if the adjacent length is zero, i.e. then Tan(Theta) is effectively infinity.

How do you deal with the cases where the two points are dead North etc? Simple investigation of the values:

p1 = 1234567890
p1 = 1234509876

The eastings are the same so all that you need to deal with are the northings to determine if P1 is north or south of P2.

boB
(It's to the North!)

In reply to Tony Little:

I want to do it since in a spreadsheet formula it'll save me days of work the hard way...

I'm not just testing.

In reply to The Crow:

Ok well to do it from lon / lat in javascript its:

/*
* Calculate distance (in km) between two points specified by latitude/longitude with Haversine formula
*
*/
LatLong.distHaversine = function(p1, p2) {
var R = 6371; // earth's mean radius in km
var dLat = p2.lat - p1.lat;
var dLong = p2.lon - p1.lon;

var a = Math.sin(dLat/2) * Math.sin(dLat/2) +
Math.cos(p1.lat) * Math.cos(p2.lat) * Math.sin(dLong/2) * Math.sin(dLong/2);
var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
var d = R * c;

return d;
}


/*
* ditto using law of cosines.
*/
LatLong.distCosineLaw = function(p1, p2) {
var R = 6371; // earth's mean radius in km
var d = Math.acos(Math.sin(p1.lat)*Math.sin(p2.lat) +
Math.cos(p1.lat)*Math.cos(p2.lat)*Math.cos(p2.lon-p1.lon)) * R;
return d;
}


/*
* LatLong constructor:
*/
function LatLong(degLat, degLong) {
this.lat = LatLong.llToRad(degLat);
this.lon = LatLong.llToRad(degLong);
}


/*
* convert lat/long in degrees to radians, for handling input values
*
* this is very flexible on formats, allowing signed decimal degrees (numeric or text), or
* deg-min-sec suffixed by compass direction (NSEW). A variety of separators are accepted
* (eg 3º 37' 09"W) or fixed-width format without separators (eg 0033709W). Seconds and minutes
* may be omitted. Minimal validation is done.
*/
LatLong.llToRad = function(llDeg) {
if (!isNaN(llDeg)) return llDeg * Math.PI / 180; // signed decimal degrees without NSEW

llDeg = llDeg.replace(/[s]*$/,''); // strip trailing whitespace
var dir = llDeg.slice(-1).toUpperCase(); // compass dir'n
if (!/[NSEW]/.test(dir)) return NaN; // check for correct compass direction
llDeg = llDeg.slice(0,-1); // and lose it off the end
var dms = llDeg.split(/[s:,°º′'″"]/); // check for separators indicating d/m/s
if (dms[dms.length-1] == '') dms.length--; // trailing separator? (see note below)
switch (dms.length) { // convert to decimal degrees...
case 3: // interpret 3-part result as d/m/s
var deg = dms[0]/1 + dms[1]/60 + dms[2]/3600; break;
case 2: // interpret 2-part result as d/m
var deg = dms[0]/1 + dms[1]/60; break;
case 1: // non-separated format dddmmss
if (/[NS]/.test(dir)) llDeg = '0' + llDeg; // - normalise N/S to 3-digit degrees
var deg = llDeg.slice(0,3)/1 + llDeg.slice(3,5)/60 + llDeg.slice(5)/3600; break;
default: return NaN;
}
if (/[WS]/.test(dir)) deg = -deg; // take west and south as -ve
return deg * Math.PI / 180; // then convert to radians
}
// note: 'x-'.split(/-/) should give ['x',''] but in IE just gives ['x']


<p>Lat 1: <input name="lat1" value="53 09 02N"> Long 1: <input name="long1" value="001 50 40W"></p>
<p>Lat 2: <input name="lat2" value="52 12 17N"> Long 2: <input name="long2" value="000 08 26E"></p>

<input type="button" value="calculate distance"
onClick="result.value = LatLong.distHaversine(new LatLong(lat1.value, long1.value),
new LatLong(lat2.value, long2.value)) + ' km'">
<input name="result" value="">


and.... (split coz its longer than the max message size)

...and to get to and from NGR <--> lon/lat its

/*
* convert geodesic co-ordinates to OS grid reference
*/
LatLong.prototype.LatLongToOSGrid = function() {
var a = 6377563.396, b = 6356256.910; // Airy 1830 major & minor semi-axes
var F0 = 0.9996012717; // NatGrid scale factor on central meridian
var lat0 = 49*Math.PI/180, lon0 = -2*Math.PI/180; // NatGrid true origin
var N0 = -100000, E0 = 400000; // northing & easting of true origin, metres

var e2 = 1 - (b*b)/(a*a); // eccentricity squared
var n = (a-b)/(a+b), n2 = n*n, n3 = n*n*n;

var lat = this.lat;
var cosLat = Math.cos(lat), sinLat = Math.sin(lat);
var nu = a*F0/Math.sqrt(1-e2*sinLat*sinLat); // transverse radius of curvature
var rho = a*F0*(1-e2)/Math.pow(1-e2*sinLat*sinLat, 1.5); // meridional radius of curvature
var eta2 = nu/rho-1;


var Ma = (1 + n + (5/4)*n2 + (5/4)*n3) * (lat-lat0);
var Mb = (3*n + 3*n*n + (21/8)*n3) * Math.sin(lat-lat0) * Math.cos(lat+lat0);
var Mc = ((15/8)*n2 + (15/8)*n3) * Math.sin(2*(lat-lat0)) * Math.cos(2*(lat+lat0));
var Md = (35/24)*n3 * Math.sin(3*(lat-lat0)) * Math.cos(3*(lat+lat0));
var M = b * F0 * (Ma - Mb + Mc - Md); // meridional arc

var cos3lat = cosLat*cosLat*cosLat;
var cos5lat = cos3lat*cosLat*cosLat;
var tan2lat = Math.tan(lat)*Math.tan(lat);
var tan4lat = tan2lat*tan2lat;

var I = M + N0;
var II = (nu/2)*sinLat*cosLat;
var III = (nu/24)*sinLat*cos3lat*(5-tan2lat+9*eta2);
var IIIA = (nu/720)*sinLat*cos5lat*(61-58*tan2lat+tan4lat);
var IV = nu*cosLat;
var V = (nu/6)*cos3lat*(nu/rho-tan2lat);
var VI = (nu/120) * cos5lat * (5 - 18*tan2lat + tan4lat + 14*eta2 - 58*tan2lat*eta2);

var dLon = this.lon-lon0;
var dLon2 = dLon*dLon, dLon3 = dLon2*dLon, dLon4 = dLon3*dLon, dLon5 = dLon4*dLon, dLon6 = dLon5*dLon;

var grid = new Object();
grid.N = I + II*dLon2 + III*dLon4 + IIIA*dLon6;
grid.E = E0 + IV*dLon + V*dLon3 + VI*dLon5;

return LatLong.gridrefNumToLet(grid.E, grid.N, 8);
}


/*
* convert OS grid reference to geodesic co-ordinates
*/
LatLong.OSGridToLatLong = function(gridRef) {
var gr = LatLong.gridrefLetToNum(gridRef);
var E = gr[0], N = gr[1];

var a = 6377563.396, b = 6356256.910; // Airy 1830 major & minor semi-axes
var F0 = 0.9996012717; // NatGrid scale factor on central meridian
var lat0 = 49*Math.PI/180, lon0 = -2*Math.PI/180; // NatGrid true origin
var N0 = -100000, E0 = 400000; // northing & easting of true origin, metres

var e2 = 1 - (b*b)/(a*a); // eccentricity squared
var n = (a-b)/(a+b), n2 = n*n, n3 = n*n*n;

var lat=lat0, M=0;
do {
lat = (N-N0-M)/(a*F0) + lat;

var Ma = (1 + n + (5/4)*n2 + (5/4)*n3) * (lat-lat0);
var Mb = (3*n + 3*n*n + (21/8)*n3) * Math.sin(lat-lat0) * Math.cos(lat+lat0);
var Mc = ((15/8)*n2 + (15/8)*n3) * Math.sin(2*(lat-lat0)) * Math.cos(2*(lat+lat0));
var Md = (35/24)*n3 * Math.sin(3*(lat-lat0)) * Math.cos(3*(lat+lat0));
var M = b * F0 * (Ma - Mb + Mc - Md); // meridional arc

} while (N-N0-M >= 0.00001); // ie until < 0.01mm

var cosLat = Math.cos(lat), sinLat = Math.sin(lat);
var nu = a*F0/Math.sqrt(1-e2*sinLat*sinLat); // transverse radius of curvature
var rho = a*F0*(1-e2)/Math.pow(1-e2*sinLat*sinLat, 1.5); // meridional radius of curvature
var eta2 = nu/rho-1;

var tanLat = Math.tan(lat);
var tan2lat = tanLat*tanLat, tan4lat = tan2lat*tan2lat, tan6lat = tan4lat*tan2lat;
var secLat = 1/cosLat;
var nu3 = nu*nu*nu, nu5 = nu3*nu*nu, nu7 = nu5*nu*nu;
var VII = tanLat/(2*rho*nu);
var VIII = tanLat/(24*rho*nu3)*(5+3*tan2lat+eta2-9*tan2lat*eta2);
var IX = tanLat/(720*rho*nu5)*(61+90*tan2lat+45*tan4lat);
var X = secLat/nu;
var XI = secLat/(6*nu3)*(nu/rho+2*tan2lat);
var XII = secLat/(120*nu5)*(5+28*tan2lat+24*tan4lat);
var XIIA = secLat/(5040*nu7)*(61+662*tan2lat+1320*tan4lat+720*tan6lat);

var dE = (E-E0), dE2 = dE*dE, dE3 = dE2*dE, dE4 = dE2*dE2, dE5 = dE3*dE2, dE6 = dE4*dE2, dE7 = dE5*dE2;
lat = lat - VII*dE2 + VIII*dE4 - IX*dE6;
var lon = lon0 + X*dE - XI*dE3 + XII*dE5 - XIIA*dE7;

return new LatLong(lat*180/Math.PI, lon*180/Math.PI);
}


/*
* convert standard grid reference ('SU387148') to fully numeric ref ([438700,114800])
* returned co-ordinates are in metres, centred on grid square for conversion to lat/long
*
* note that northern-most grid squares will give 7-digit northings
* no error-checking is done on gridref (bad input will give bad results or NaN)
*/
LatLong.gridrefLetToNum = function(gridref) {
// get numeric values of letter references, mapping A->0, B->1, C->2, etc:
var l1 = gridref.toUpperCase().charCodeAt(0) - 'A'.charCodeAt(0);
var l2 = gridref.toUpperCase().charCodeAt(1) - 'A'.charCodeAt(0);
// shuffle down letters after 'I' since 'I' is not used in grid:
if (l1 > 7) l1--;
if (l2 > 7) l2--;

// convert grid letters into 100km-square indexes from false origin (grid square SV):
var e = ((l1-2)%5)*5 + (l2%5);
var n = (19-Math.floor(l1/5)*5) - Math.floor(l2/5);

// skip grid letters to get numeric part of ref, stripping any spaces:
gridref = gridref.slice(2).replace(/ /g,'');

// append numeric part of references to grid index:
e += gridref.slice(0, gridref.length/2);
n += gridref.slice(gridref.length/2);

// normalise to 1m grid, rounding up to centre of grid square:
switch (gridref.length) {
case 6: e += '50'; n += '50'; break;
case 8: e += '5'; n += '5'; break;
// 10-digit refs are already 1m
}

return [e, n];
}


/*
* convert numeric grid reference (in metres) to standard-form grid ref
*/
LatLong.gridrefNumToLet = function(e, n, digits) {
// get the 100km-grid indices
var e100k = Math.floor(e/100000), n100k = Math.floor(n/100000);

if (e100k<0 || e100k>6 || n100k<0 || n100k>12) return '';

// translate those into numeric equivalents of the grid letters
var l1 = (19-n100k) - (19-n100k)%5 + Math.floor((e100k+10)/5);
var l2 = (19-n100k)*5%25 + e100k%5;

// compensate for skipped 'I' and build grid letter-pairs
if (l1 > 7) l1++;
if (l2 > 7) l2++;
var letPair = String.fromCharCode(l1+'A'.charCodeAt(0), l2+'A'.charCodeAt(0));

// strip 100km-grid indices from easting & northing, and reduce precision
e = Math.floor((e%100000)/Math.pow(10,5-digits/2));
n = Math.floor((n%100000)/Math.pow(10,5-digits/2));

var gridRef = letPair + e.padLZ(digits/2) + n.padLZ(digits/2);

return gridRef;
}


/*
* pad a number with sufficient leading zeros to make it w chars wide
*/
Number.prototype.padLZ = function(w) {
var n = this.toString();
for (var i=0; i<w-n.length; i++) n = '0' + n;
return n;
}

In reply to Bob:

I understand the old SOHCAHTOA bit from schooldays but what is simple investigation of values and can I do this in a spreadsheet with a function?

Arse I wish I was smarter!

In reply to chrisw:

Bloody hell!

I'm not going to manage that in Excel am I? Even with macros...

Thanks v' much though, I might put that through to someone in IT.

In reply to The Crow:

no, its not quite as simple as it sounds. If you had the lon lat values in an excel spreadsheet I can probably come up with a simplified version that will give fairly accurate results for points upto about 200km apart (after that cumulative error gets too bad). Its the NGR to lon/lat that takes most of the work.

Have you not thought about using the google earth API?? OK not a spreadsheet directly but could be used to scrape data into CSV and into a spreadsheet from there perhaps??

In reply to The Crow: Assume the grids are abcdefghij & qrstuvwxyz.

You need E = qrstu - abcde and N = vwxyz - fghij

Firstly if N=0:
- if E is +ve Bearing = 90 degrees
- if E is -ve Bearing = 270 degrees

For all other cases:
Bearing = inverse tan (E/N) + correction.

if E +ve and N +ve - no correction
if E +ve and N -ve - add 180 degrees
if E -ve and N -ve - add 180 degrees
if E -ve and N +ve - add 360 degrees

I think I've got that all correct but can't promise anything.

In reply to chrisw:

Points over 50km apart don't matter and I can provide the OSGB grids as a lat/log easily (I've a spreadsheet that does that - not that I wrote the functions myself).

I don't suppose? ;oP

In reply to Mark Stevenson:

Thanks Mark I'll try that an see how it pans out...

Mail me an email addr that accepts attachments... I will see what I can do. Might not be tonight tho.

Hint:: When not chargin ppl hundreds of pounds an hour for IT consultancy I enjoy beer, wine, free BMC membership for a year etc

In reply to The Crow:

I think what Bob meant was inspection of the values

ie testing whether either the Easting or the Northing is identical - if so then the bearing will be one of the Cardinal points of the compass (N S E or W) depending on whether the remaining part (easting or northing as the case may be) of P1 is less than or greater than P2

In reply to Anonymous:

Indeed yes

To The Crow:

Are you intending doing something along the lines of an electronic route card? If so I have one somewhere that I did in Excel - Allows entry of points, determines if the length of each individual length is perhaps too long for accurate navigation (configurable); calculates bearings; calculates timings according to Naismith with Tranter's variations if you need them.

boB

In reply to The Crow:

I assume you are just trying to find the Grid bearing from OS Co-ords.

dE / dN gives the TAN where dE is the difference in Eastings between the two points etc.

Simple logic test to sort the quadrant out based upon the sign of dE and dN - if both +ve then Bearing = TAN. If both -ve then Brg = 180 - TAN. If dE +ve and dN -ve then 180 + TAN and if dN +ve and dE -ve then 360+ TAN.

If you want to calc true bearing then you will need to allow for t-T correction but no point unless a Geodetic issue. Also MSL correction element. If you have lat and long instead of OSGB Co-ords then do the calc on the Spheroid since its a Spheroidal Co-ordinate system.

But I guess it is as I noted at the begining, a simple bearing from co-ords on the OS grid and ignoring all the errors because its a projection.

In reply to ian2u: Correction - if both dE and dN -ve then Brg = 180 + TAN
TAN is the ATAN of dE/dN and the signs are rigurous.

In reply to The Crow: Ok, I must be bored this evening because I've written a spreadsheet formula that works and handles due N,S,E and W with logic. Unfortunately it takes up about 15 cells to do it...I can e-mail it if you like.

In reply to The Crow: OK, not quite as big a task as expected. 1 spreadsheet ready when you are.

In reply to chrisw: snap!

In reply to The Crow:

=IF((LEFT(C6,5)-LEFT(B6,5))=0,IF(RIGHT(C6,5)>RIGHT(B6,5),0,IF(RIGHT(C6,5)=RIGHT(B6,5),"nowhere",180)),90-DEGREES(ATAN((RIGHT(C6,5)-RIGHT(B6,5))/(LEFT(C6,5)-LEFT(B6,5))))+IF((LEFT(C6,5)-LEFT(B6,5))<0,180,0))

where:

B6 = 10 digit from coords
C6 = 10 digit to coords

In reply to satori:

ps. B6 and C6 need to be formatted as text to preserve any leasding zeros.

eg. 0012367890 rather than 12367890

In reply to satori:

Yup thats pretty much what I came back with - being 10 fig means you dont have to fudge the letters back to numbers and as distances are < 50km there is no need to take projection variation into account - ie no real need to convert back to lon/lat for complex calcs. String slicing and the basic tan calcs are enough.

In reply to chrisw:

Mine was =CONCATENATE(IF((LEFT(B4;5)-LEFT(A4;5))=0;IF(RIGHT(B4;5)>RIGHT(A4;5);0;IF(RIGHT(B4;5)=RIGHT(A4;5);"ERR";180));90-DEGREES(ATAN((RIGHT(B4;5)-RIGHT(A4;5))/(LEFT(B4;5)-LEFT(A4;5))))+IF((LEFT(B4;5)-LEFT(A4;5))<0;180;0));CHAR(186))

Based on the example in the Advanced Algorithms book.

Checked on a 1km square I get:

From To Bearing
3350030000 3360030000 90°
3350030000 3350029900 180°
3360030000 3350030000 270°
3350029900 3350030000 0°
3350029900 3360030000 45°

Sorry for poor cut n paste

In reply to all of you guys:

Wow. Thanks very much all of you. It's not for aroute card as it happens. I need the distances and bearings from a protected habitat to various agricultural air pollution sources to put into an atmospheric deposition model. You've all saved me lots and lots of repetitive work.

:oD

In reply to chrisw:

can remember years and years ago writing something for an air traffic radar simulator (only for educational purposes)

guts of it was rectangular to polar coordinate conversion and back again

in the beginning we had a strange effect where an aircraft would "split" into 2 equal and opposite tracks on crossing a quadrant boundary. turned out that we were using variables of insufficient precision and the conversion error alternated between quadrants