Converting UK Ordnance Survey Grid Ref to Lat/Long

[mike_eddie]

I am writing an application where I want to plot a UK grid ref on a google map. To do this I need to convert from the grid ref to lat/long, and then display this using the google maps API in a web viewer.

 

I have found the following javascript http://www.movable-type.co.uk/scripts/latlong-gridref.html which appears to do the trick. Particulaurly, this part of the code:

/**
 * Converts Ordnance Survey grid reference easting/northing coordinate to (OSGB36) latitude/longitude
 *
 * @param   {OsGridRef} gridref - Easting/northing to be converted to latitude/longitude.
 * @returns {LatLonE}   Latitude/longitude (in OSGB36) of supplied grid reference.
 *
 * @example
 *   var grid = new OsGridRef(651409, 313177);
 *   var p = OsGridRef.osGridToLatLon(grid); // p.toString(): 52°39′27″N, 001°43′04″E
 */
OsGridRef.osGridToLatLon = function(gridref) {
    var E = gridref.easting + 0.5;  // easting of centre of 1m grid square
    var N = gridref.northing + 0.5; // northing of centre of 1m grid square

    var a = 6377563.396, b = 6356256.909;              // Airy 1830 major & minor semi-axes
    var F0 = 0.9996012717;                             // NatGrid scale factor on central meridian
    var φ0 = 49*Math.PI/180, λ0 = -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+, n2 = n*n, n3 = n*n*n;         // n, n², n³

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

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

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

    var cosφ = Math.cos(φ), sinφ = Math.sin(φ);
    var ν = a*F0/Math.sqrt(1-e2*sinφ*sinφ);            // nu = transverse radius of curvature
    var ρ = a*F0*(1-e2)/Math.pow(1-e2*sinφ*sinφ, 1.5); // rho = meridional radius of curvature
    var η2 = ν/ρ-1;                                    // eta = ?

    var tanφ = Math.tan(φ);
    var tan2φ = tanφ*tanφ, tan4φ = tan2φ*tan2φ, tan6φ = tan4φ*tan2φ;
    var secφ = 1/cosφ;
    var ν3 = ν*ν*ν, ν5 = ν3*ν*ν, ν7 = ν5*ν*ν;
    var VII = tanφ/(2*ρ*ν);
    var VIII = tanφ/(24*ρ*ν3)*(5+3*tan2φ+η2-9*tan2φ*η2);
    var IX = tanφ/(720*ρ*ν5)*(61+90*tan2φ+45*tan4φ);
    var X = secφ/ν;
    var XI = secφ/(6*ν3)*(ν/ρ+2*tan2φ);
    var XII = secφ/(120*ν5)*(5+28*tan2φ+24*tan4φ);
    var XIIA = secφ/(5040*ν7)*(61+662*tan2φ+1320*tan4φ+720*tan6φ);

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

    return new LatLonE(φ.toDegrees(), λ.toDegrees(), LatLonE.datum.OSGB36);
};

Ive tried to recreate the above logic in a script in Filemaker. I seem to be getting an answer which is nearly correct, but not quite. Ive checked and double checked my version of this and I cannot see any errors. Could this be a rounding issue, or perhaps Im missing something obvious!

Convert OS Grid Ref to Lat Lon
Set Variable [ $E; Value:Temp::gE + .5 ]
Set Variable [ $N; Value:Temp::gN + .5 ] Set Variable [ $a; Value:6377563.396 ] Set Variable [ $b; Value:6356256.909 ] Set Variable [ $F0; Value:.9996012717 ] Set Variable [ $phi0; Value:49*Pi/180 ] Set Variable [ $delta0; Value:-2*Pi/180 ] Set Variable [ $N0; Value:-100000 ]
Set Variable [ $E0; Value:400000 ]
Set Variable [ $e2; Value:1 - ($b*$b)/($a*$a) ] Set Variable [ $n; Value:($a-$b)/($a+$ ]
Set Variable [ $n2; Value:$n*$n ]
Set Variable [ $n3; Value:$n*$n*$n ]
Set Variable [ $phi; Value:$phi0 ]
Set Variable [ $M; Value:0 ]
Loop
Set Variable [ $phi; Value:($N-$N0-$M)/($a*$F0) + $phi ]
Set Variable [ $Ma; Value:(1 + $n + (5/4)*$n2 + (5/4)*$n3) * ($phi-$phi0) ]
Set Variable [ $Mb; Value:(3*$n + 3*$n*$n + (21/8)*$n3) * Sin ( $phi-$phi0 ) * Cos($phi+$phi0) ] Set Variable [ $Mc; Value:((15/8)*$n2 + (15/8)*$n3) * Sin(2*($phi-$phi0)) * Cos(2*($phi+$phi0)) ] Set Variable [ $Md; Value:(35/24)*$n3 * Sin(3*($phi-$phi0)) * Cos(3*($phi+$phi0)) ]
Set Variable [ $M; Value:$b * $F0 * ($Ma - $Mb + $Mc - $Md) ]
Exit Loop If [ $N-$N0-$M >= .00001 ]
End Loop
Set Variable [ $cos_phi; Value:Cos($phi) ]
Set Variable [ $sin_phi; Value:Sin($phi) ]
Set Variable [ $v; Value:$a*$F0/Sqrt(1-$e2*$sin_phi*$sin_phi) ]
Set Variable [ $rho; Value:$a*$F0*(1-$e2)/((1-$e2*$sin_phi*$sin_phi)^1.5) ]
Set Variable [ $eta2; Value:$v/$rho-1 ]
Set Variable [ $tan_phi; Value:Tan($phi) ]
Set Variable [ $tan2phi; Value:$tan_phi * $tan_phi ]
Set Variable [ $tan4phi; Value:$tan2phi * $tan2phi ]
Set Variable [ $tan6phi; Value:$tan4phi * $tan2phi ]
Set Variable [ $sec_phi; Value:1 / $cos_phi ]
Set Variable [ $v3; Value:$v * $v * $v ]
Set Variable [ $v5; Value:$v3 * $v * $v ]
Set Variable [ $v7; Value:$v5 * $v * $v ]
Set Variable [ $VII; Value:$tan_phi / (2*$rho*$v) ]
Set Variable [ $VIII; Value:$tan_phi/(24*$rho*$v3)*(5+3*$tan2phi+$eta2-9*$tan2phi*$eta2) ]
Set Variable [ $IX; Value:$tan_phi/(720*$rho*$v5)*(61+90*$tan2phi+45*$tan4phi) ]
Set Variable [ $X; Value:$sec_phi/$v ]
Set Variable [ $XI; Value:$sec_phi/(6*$v3)*($v/$rho+2*$tan2phi) ]
Set Variable [ $XII; Value:$sec_phi/(120*$v5)*(5+28*$tan2phi+24*$tan4phi) ]
Set Variable [ $XIIA; Value:$sec_phi/(5040*$v7)*(61+662*$tan2phi+1320*$tan4phi+720*$tan6phi) ] Set Variable [ $dE; Value:($E-$E0) ]
Set Variable [ $dE2; Value:$dE*$dE ]
Set Variable [ $dE3; Value:$dE2 * $dE ]
Set Variable [ $dE4; Value:$dE2 * $dE2 ]
Set Variable [ $dE5; Value:$dE2 * $dE3 ]
Set Variable [ $dE6; Value:$dE4 * $dE2 ]
Set Variable [ $dE7; Value:$dE5 * $dE2 ]
Set Variable [ $phi; Value:$phi - $VII*$dE2 + $VIII*$dE4 - $IX*$dE6 ]
Set Field [ Incident::Latitude; Degrees($phi) ]
Set Field [ Incident::Longitude; Degrees($delta) ]


Alternatively, is there a way to run the original javascript file, pass in my grid ref from the database, convert to lat/long, save this on the record, and then use the google maps API to plot this accordingly.

 

Any pointers would be hugely appreciated!

 

Mike

[jbante]

I haven't looked at your implementation in detail to see if you've translated the math correctly. However, my experience with coordinate systems is that if you're math is nearly correct, your math is probably exactly correct, and there's something more subtle at work throwing things off.

 

First, how do you know that what you're getting isn't quite right? I presume that you're plotting results in a map, and they aren't what you expect?

 

I noticed that the JavaScript you're borrowing from mentions that it's using the OSGB36 datum. Google Maps uses WGS84. (Roughly. Some organizations don't think the projection used by Google Maps is accurate enough.) You'll want to do a transformation between the two before plotting anything in Google Maps.

[mike_eddie]

Thanks for the reply!

 

I think there may be something wrong with the maths somewhere. I had entered a coordinate for north devon, and it was plotting me somewhere in the english channel. Not that far away in the grand scheme of things. However I tried a couple more coordinates up in scotland and the position only moves slightly. 

 

I just can't figure out what is going wrong! Could it be to do with the number of times the loop is being performed? All the other maths seems correct?

[comment]

Are you getting the same results with your script as you do with the original javascript?

[jbante]

Mike, first, follow Comment's suggestion: compare your results with the results of the original JavaScript.

 

After that's done, look at the datum conversion stuff again. What you're describing is consistent with what I was talking about with the different coordinate datums. Being off by different distances in different places is exactly what you should expect if you're trying to plot OSGB36 coordinates on a WGS84 map.

[comment]

What you're describing is consistent with what I was talking about with the different coordinate datums.

 

That too should be very easy to confirm: the linked page has a demo that returns both OSGB36 and WGS84 lat/long.

[mike_eddie]

Thanks guys. 

Finally found the main problem.

 

Java treats N and n as different variables, whereas filemaker treats $n and $N as the same variable, so things were getting over written.

 

There is now a small difference between the two, plotting on opposite sides of main roads etc which I think is due to the OSGB36 vs WGS84 problem as you suggest. Once Ive written the script to convert these Ill add the code here for anyone looking in future!

 

Thanks again for your help!

 

BW,

 

Mike

[calumma]

Mike, did you make any progress with the OSGB36 to WGS84 conversion? I had a go at writing a script to do this some while ago, but my maths is clearly not up to the task! I use a database for a small UK conservation group and plot species distributions in Google Earth. Since I couldn't do the conversion myself I found an online resource that allows me to scrape data into the database. Not an ideal solution though and when I batch import records, the conversions can swamp the server. Implementing your solution would make things much simpler! 

 

Cheers

Lee