More on Converting British National Grid to Latitude and Longitude

So it turns out the internet is awesome.

After I posted a python script that can be used to convert a .csv of Northings and Eastings in British National Grid to a Latitude and Longitude, several people followed suit and sent through their codes in a variety of other languages. 

I've just changed over my website and didn't want to lose the codes, so I've put them all here in one handy place. Hope you find them useful.

 

Objective C

Written by David Gibson


// Convert OSGB36 easting/northing to WSG84 latitude and longitude
– (CLLocation*) latLon_WSG84
{
double E = self.easting;
double N = self.northing;

// E, N are the British national grid coordinates – eastings and northings
double a = 6377563.396, b = 6356256.909; // The Airy 180 semi-major and semi-minor axes used for OSGB36 (m)
double F0 = 0.9996012717; // scale factor on the central meridian
double lat0 = 49*M_PI/180; // Latitude of true origin (radians)
double lon0 = -2*M_PI/180; // Longtitude of true origin and central meridian (radians)
double N0 = -100000, E0 = 400000; // Northing & easting of true origin (m)
double e2 = 1 – (b*b)/(a*a); // eccentricity squared
double n = (a-b)/(a+b);

// Initialise the iterative variables
double lat = lat0, M = 0;

while(N-N0-M >= 0.00001) // Accurate to 0.01mm
{
lat = (N-N0-M)/(a*F0) + lat;
double M1 = (1. + n + (5/4)*pow(n,2) + (5/4)*pow(n, 3)) * (lat-lat0);
double M2 = (3.*n + 3.*pow(n,2) + (21/8)*pow(n,3)) * sin(lat-lat0) * cos(lat+lat0);
double M3 = ((15/8)*pow(n,2) + (15/8)*pow(n,3)) * sin(2*(lat-lat0)) * cos(2*(lat+lat0));
double M4 = (35/24)*pow(n,3) * sin(3*(lat-lat0)) * cos(3*(lat+lat0));

// meridional arc
M = b * F0 * (M1 – M2 + M3 – M4);
}

// transverse radius of curvature
double nu = a*F0/sqrt(1-e2*pow(sin(lat),2));

// meridional radius of curvature
double rho = a*F0*(1-e2)*pow((1-e2*pow(sin(lat),2)),(-1.5));
double eta2 = nu/rho-1.;

double secLat = 1./cos(lat);
double VII = tan(lat)/(2*rho*nu);
double VIII = tan(lat)/(24*rho*pow(nu,3))*(5+3*pow(tan(lat),2)+eta2-9*pow(tan(lat),2)*eta2);
double IX = tan(lat)/(720*rho*pow(nu,5))*(61+90*pow(tan(lat),2)+45*pow(tan(lat),4));
double X = secLat/nu;
double XI = secLat/(6*pow(nu,3))*(nu/rho+2*pow(tan(lat),2));
double XII = secLat/(120*pow(nu,5))*(5+28*pow(tan(lat),2)+24*pow(tan(lat),4));
double XIIA = secLat/(5040*pow(nu,7))*(61+662*pow(tan(lat),2)+1320*pow(tan(lat),4)+720*pow(tan(lat),6));
double dE = E-E0;

// These are on the wrong ellipsoid currently: Airy1830. (Denoted by _1)
double lat_1 = lat – VII*pow(dE,2) + VIII*pow(dE,4) – IX*pow(dE,6);
double lon_1 = lon0 + X*dE – XI*pow(dE,3) + XII*pow(dE,5) – XIIA*pow(dE,7);

// Obj-C log
NSLog(@”Firstbash %f, %f”, lat_1*180/M_PI, lon_1*180/M_PI);

// Want to convert to the GRS80 ellipsoid.
// First convert to cartesian from spherical polar coordinates
double H = 0; // Third spherical coord.
double x_1 = (nu/F0 + H)*cos(lat_1)*cos(lon_1);
double y_1 = (nu/F0+ H)*cos(lat_1)*sin(lon_1);
double z_1 = ((1-e2)*nu/F0 +H)*sin(lat_1);

// Perform Helmut transform (to go between Airy 1830 (_1) and GRS80 (_2))
double s = -20.4894*pow(10,-6); // The scale factor -1
double tx = 446.448, ty = -125.157, tz = 542.060; // The translations along x,y,z axes respectively
double rxs = 0.1502, rys = 0.2470, rzs = 0.8421; // The rotations along x,y,z respectively, in seconds
double rx = rxs*M_PI/(180*3600.), ry = rys*M_PI/(180*3600.), rz = rzs*M_PI/(180*3600.); // In radians
double x_2 = tx + (1+s)*x_1 + (-rz)*y_1 + (ry)*z_1;
double y_2 = ty + (rz)*x_1 + (1+s)*y_1 + (-rx)*z_1;
double z_2 = tz + (-ry)*x_1 + (rx)*y_1 + (1+s)*z_1;

// Back to spherical polar coordinates from cartesian
// Need some of the characteristics of the new ellipsoid
double a_2 =6378137.000, b_2 = 6356752.3141; // The GSR80 semi-major and semi-minor axes used for WGS84(m)
double e2_2 = 1- (b_2*b_2)/(a_2*a_2); // The eccentricity of the GRS80 ellipsoid
double p = sqrt(pow(x_2,2) + pow(y_2,2));

// Lat is obtained by an iterative proceedure:
lat = atan2(z_2,(p*(1-e2_2))); // Initial value
double latold = 2*M_PI;

double nu_2 = 0.;

while(abs(lat – latold)>pow(10,-16))
{
double latTemp = lat;
lat = latold;
latold = latTemp;
nu_2 = a_2/sqrt(1-e2_2*pow(sin(latold),2));
lat = atan2(z_2+e2_2*nu_2*sin(latold), p);
}

// Lon and height are then pretty easy
double lon = atan2(y_2,x_2);
H = p/cos(lat) – nu_2;

// Obj-C log
NSLog(@”%f, %f”, (lat-lat_1)*180/M_PI, (lon – lon_1)*180/M_PI);

// Convert to degrees
lat = lat*180/M_PI;
lon = lon*180/M_PI;

// create Obj-C location object – alternatively, output the ‘lat’ and ‘lon’ variables above
CLLocation* loc = [[CLLocation alloc] initWithLatitude:lat longitude:lon];

return loc;
}

Javascript

Written by Chris Dack

function OSGB36toWGS84(E,N){

//E, N are the British national grid coordinates - eastings and northings
var a = 6377563.396;
var b = 6356256.909; //The Airy 180 semi-major and semi-minor axes used for OSGB36 (m)
var F0 = 0.9996012717; //scale factor on the central meridian
var lat0 = 49*Math.PI/180; //Latitude of true origin (radians)
var lon0 = -2*Math.PI/180; //Longtitude of true origin and central meridian (radians)
var N0 = -100000;
var E0 = 400000; //Northing & easting of true origin (m)
var e2 = 1 - (b*b)/(a*a); //eccentricity squared
var n = (a-b)/(a+b);

//Initialise the iterative variables
lat = lat0;
M = 0;

while (N-N0-M >= 0.00001){ //Accurate to 0.01mm

lat = (N-N0-M)/(a*F0) + lat;
M1 = (1 + n + (5./4)*Math.pow(n,2) + (5./4)*Math.pow(n,3)) * (lat-lat0);
M2 = (3*n + 3*Math.pow(n,2) + (21./8)*Math.pow(n,3)) * Math.sin(lat-lat0) * Math.cos(lat+lat0);
M3 = ((15./8)*Math.pow(n,2) + (15./8)*Math.pow(n,3)) * Math.sin(2*(lat-lat0)) * Math.cos(2*(lat+lat0));
M4 = (35./24)*Math.pow(n,3) * Math.sin(3*(lat-lat0)) * Math.cos(3*(lat+lat0));
//meridional arc
M = b * F0 * (M1 - M2 + M3 - M4) ;
}
//transverse radius of curvature
var nu = a*F0/Math.sqrt(1-e2*Math.pow(Math.sin(lat),2));

//meridional radius of curvature

var rho = a*F0*(1-e2)*Math.pow((1-e2*Math.pow(Math.sin(lat),2)),(-1.5));
var eta2 = nu/rho-1
var secLat = 1./Math.cos(lat);
var VII = Math.tan(lat)/(2*rho*nu);
var VIII = Math.tan(lat)/(24*rho*Math.pow(nu,3))*(5+3*Math.pow(Math.tan(lat),2)+eta2-9*Math.pow(Math.tan(lat),2)*eta2);
var IX = Math.tan(lat)/(720*rho*Math.pow(nu,5))*(61+90*Math.pow(Math.tan(lat),2)+45*Math.pow(Math.tan(lat),4));
var X = secLat/nu;
var XI = secLat/(6*Math.pow(nu,3))*(nu/rho+2*Math.pow(Math.tan(lat),2));
var XII = secLat/(120*Math.pow(nu,5))*(5+28*Math.pow(Math.tan(lat),2)+24*Math.pow(Math.tan(lat),4));
var XIIA = secLat/(5040*Math.pow(nu,7))*(61+662*Math.pow(Math.tan(lat),2)+1320*Math.pow(Math.tan(lat),4)+720*Math.pow(Math.tan(lat),6));
var dE = E-E0;

//These are on the wrong ellipsoid currently: Airy1830. (Denoted by _1)
var lat_1 = lat - VII*Math.pow(dE,2) + VIII*Math.pow(dE,4) - IX*Math.pow(dE,6);
var lon_1 = lon0 + X*dE - XI*Math.pow(dE,3) + XII*Math.pow(dE,5) - XIIA*Math.pow(dE,7);

//Want to convert to the GRS80 ellipsoid.
//First convert to cartesian from spherical polar coordinates
var H = 0 //Third spherical coord.
var x_1 = (nu/F0 + H)*Math.cos(lat_1)*Math.cos(lon_1);
var y_1 = (nu/F0+ H)*Math.cos(lat_1)*Math.sin(lon_1);
var z_1 = ((1-e2)*nu/F0 +H)*Math.sin(lat_1);

//Perform Helmut transform (to go between Airy 1830 (_1) and GRS80 (_2))
var s = -20.4894*Math.pow(10,-6); //The scale factor -1
var tx = 446.448; //The translations along x,y,z axes respectively
var ty = -125.157;
var tz = 542.060;
var rxs = 0.1502;
var rys = 0.2470;
var rzs = 0.8421; //The rotations along x,y,z respectively, in seconds
var rx = rxs*Math.PI/(180*3600);
var ry= rys*Math.PI/(180*3600);
var rz = rzs*Math.PI/(180*3600); //In radians
var x_2 = tx + (1+s)*x_1 + (-rz)*y_1 + (ry)*z_1;
var y_2 = ty + (rz)*x_1 + (1+s)*y_1 + (-rx)*z_1;
var z_2 = tz + (-ry)*x_1 + (rx)*y_1 + (1+s)*z_1;

//Back to spherical polar coordinates from cartesian
//Need some of the characteristics of the new ellipsoid
var a_2 = 6378137.000;
var b_2 = 6356752.3141; //The GSR80 semi-major and semi-minor axes used for WGS84(m)
var e2_2 = 1- (b_2*b_2)/(a_2*a_2); //The eccentricity of the GRS80 ellipsoid
var p = Math.sqrt(Math.pow(x_2,2) + Math.pow(y_2,2));

//Lat is obtained by an iterative proceedure:
var lat = Math.atan2(z_2,(p*(1-e2_2))); //Initial value
var latold = 2*Math.PI;

while (Math.abs(lat - latold)>Math.pow(10,-16)){
//console.log(Math.abs(lat - latold));
var temp;
var temp2;
var nu_2;
temp1 = lat;
temp2 = latold;
latold = temp1;
lat =temp2;

lat = latold;
latold = lat;
nu_2 = a_2/Math.sqrt(1-e2_2*Math.pow(Math.sin(latold),2));
lat = Math.atan2(z_2+e2_2*nu_2*Math.sin(latold), p);
}
//Lon and height are then pretty easy
var lon = Math.atan2(y_2,x_2);
var H = p/Math.cos(lat) - nu_2;

//Convert to degrees
lat = lat*180/Math.PI;
lon = lon*180/Math.PI;

//Job's a good'n.
return [lat,lon];
}

C#

Written by Andy Nichols


using System;

namespace Whatever.Namespace.You.Want
{
//adapted from the python script at
//http://hannahfry.co.uk/2012/02/01/converting-british-national-grid-to-latitude-and-longitude-ii/
//where variables have non-descriptive names it's because I couldn't see what they represent
public static class GeoCoordinatesConverter
{
private const double Airy180SemiMajorAxis = 6377563.396;
private const double Airy180SemiMinorAxis = 6356256.909;
private const double ScaleFactorOnCentralMeridian = 0.9996012717;
private const int NorthingOfTrueOrigin = -100000;
private const int EastingOfTrueOrigin = 400000;
private const double LatitudeOfTrueOrigin = 49 * Math.PI / 180;
private const double LongtitudeOfTrueOrigin = -2 * Math.PI / 180;
private const double N = (Airy180SemiMajorAxis - Airy180SemiMinorAxis) / (Airy180SemiMajorAxis + Airy180SemiMinorAxis);
private const double Grs80SemiMajorAxis = 6378137.000;
private const double Grs80SemiMinorAxis = 6356752.3141;
private const double ScaleFactor = -0.0000204894;
private const double XRadians = 0.1502 * Math.PI / 648000;
private const double YRadians = 0.2470 * Math.PI / 648000;
private const double ZRadians = 0.8421 * Math.PI / 648000;

public static WGS84 Convert(OSGB36 coordinates)
{
var latitude = InitializeLatitude(coordinates);
var grs80 = CalculateGrs80(coordinates, latitude);
latitude = RecalculateLatitude(grs80);

return new WGS84
{
Latitude = latitude * 180 / Math.PI,
Longtitude = Math.Atan2(grs80.Y, grs80.X) * 180 / Math.PI
};
}

private static double InitializeLatitude(OSGB36 coordinates)
{
var latitude = LatitudeOfTrueOrigin;
var meridionalArc = 0.0;
while (coordinates.Northing - NorthingOfTrueOrigin - meridionalArc >= 0.00001)
{
latitude += (coordinates.Northing - NorthingOfTrueOrigin - meridionalArc) /
(Airy180SemiMajorAxis * ScaleFactorOnCentralMeridian);
var meridionalArc1 = (1.0 + N + (1.25 * Math.Pow(N, 2)) + (1.25 * Math.Pow(N, 3))) *
(latitude - LatitudeOfTrueOrigin);
var meridionalArc2 = ((3 * N) + (3 * Math.Pow(N, 2)) + ((21.0 / 8) * Math.Pow(N, 3))) *
Math.Sin(latitude - LatitudeOfTrueOrigin) *
Math.Cos(latitude + LatitudeOfTrueOrigin);
var meridionalArc3 = ((15.0 / 8) * Math.Pow(N, 2) + (15.0 / 8) * Math.Pow(N, 3)) *
Math.Sin(2 * (latitude - LatitudeOfTrueOrigin)) *
Math.Cos(2 * (latitude + LatitudeOfTrueOrigin));
var meridionalArc4 = (35.0 / 24) * Math.Pow(N, 3) * Math.Sin(3 * (latitude - LatitudeOfTrueOrigin)) *
Math.Cos(3 * (latitude + LatitudeOfTrueOrigin));
meridionalArc = Airy180SemiMinorAxis * ScaleFactorOnCentralMeridian *
(meridionalArc1 - meridionalArc2 + meridionalArc3 - meridionalArc4);
}
return latitude;
}

private static Cartesian CalculateGrs80(OSGB36 coordinates, double latitude)
{
var eccentricitySquared = 1 - (Math.Pow(Airy180SemiMinorAxis, 2) / Math.Pow(Airy180SemiMajorAxis, 2));
var transverseRadiusOfCurvature = Airy180SemiMajorAxis * ScaleFactorOnCentralMeridian /
Math.Sqrt(1 - (eccentricitySquared * latitude.SinToPower(2)));

var airy1830 = CalculateAiry1830(coordinates, latitude, transverseRadiusOfCurvature, eccentricitySquared);
var cartesian = new Cartesian
{
X = (transverseRadiusOfCurvature / ScaleFactorOnCentralMeridian) * Math.Cos(airy1830.Latitude) * Math.Cos(airy1830.Longtitude),
Y = (transverseRadiusOfCurvature / ScaleFactorOnCentralMeridian) * Math.Cos(airy1830.Latitude) * Math.Sin(airy1830.Longtitude),
Z = ((1 - eccentricitySquared) * transverseRadiusOfCurvature / ScaleFactorOnCentralMeridian) * Math.Sin(airy1830.Latitude)
};
var grs80 = new Cartesian
{
X = 446.448 + (1 + ScaleFactor) * cartesian.X - ZRadians * cartesian.Y + YRadians * cartesian.Z,
Y = -125.157 + ZRadians * cartesian.X + (1 + ScaleFactor) * cartesian.Y - XRadians * cartesian.Z,
Z = 542.060 - YRadians * cartesian.X + XRadians * cartesian.Y + (1 + ScaleFactor) * cartesian.Z
};
return grs80;
}

private static LatitudeLongtitude CalculateAiry1830(OSGB36 coordinates, double latitude,
double transverseRadiusOfCurvature, double eccentricitySquared)
{
var meridionalRadiusOfCurvature = Airy180SemiMajorAxis * ScaleFactorOnCentralMeridian * (1 - eccentricitySquared) * Math.Pow((1 - (eccentricitySquared * latitude.SinToPower(2))), -1.5);
var eta2 = (transverseRadiusOfCurvature / meridionalRadiusOfCurvature) - 1;
var secLat = 1.0 / Math.Cos(latitude);
var vii = Math.Tan(latitude) / (2 * meridionalRadiusOfCurvature * transverseRadiusOfCurvature);
var viii = Math.Tan(latitude) / (24 * meridionalRadiusOfCurvature * Math.Pow(transverseRadiusOfCurvature, 3)) *
(5 + 3 * latitude.TanToPower(2) + eta2 - 9 * latitude.TanToPower(2) * eta2);
var ix = Math.Tan(latitude) / (720 * meridionalRadiusOfCurvature * Math.Pow(transverseRadiusOfCurvature, 5)) *
(61 + 90 * latitude.TanToPower(2) + 45 * latitude.TanToPower(4));
var x = secLat / transverseRadiusOfCurvature;
var xi = secLat / (6 * Math.Pow(transverseRadiusOfCurvature, 3)) *
(transverseRadiusOfCurvature / meridionalRadiusOfCurvature + 2 * latitude.TanToPower(2));
var xii = secLat / (120 * Math.Pow(transverseRadiusOfCurvature, 5)) *
(5 + 28 * latitude.TanToPower(2) + 24 * latitude.TanToPower(4));
var xiia = secLat / (5040 * Math.Pow(transverseRadiusOfCurvature, 7)) *
(61 + 662 * latitude.TanToPower(2) + 1320 * latitude.TanToPower(4) + 720 * latitude.TanToPower(6));

var eastingDifference = coordinates.Easting - EastingOfTrueOrigin;
return new LatitudeLongtitude
{
Latitude = latitude - vii * Math.Pow(eastingDifference, 2) + viii * Math.Pow(eastingDifference, 4) -
ix * Math.Pow(eastingDifference, 6),
Longtitude = LongtitudeOfTrueOrigin + x * eastingDifference - xi * Math.Pow(eastingDifference, 3) +
xii * Math.Pow(eastingDifference, 5) - xiia * Math.Pow(eastingDifference, 7)
};
}

private static double RecalculateLatitude(Cartesian grs80)
{
var eccentricityOfGrs80Ellipsoid = 1 - Math.Pow(Grs80SemiMinorAxis, 2) / Math.Pow(Grs80SemiMajorAxis, 2);
var sphericalPolar = Math.Sqrt(Math.Pow(grs80.X, 2) + Math.Pow(grs80.Y, 2));

var latitude = Math.Atan2(grs80.Z, (sphericalPolar * (1 - eccentricityOfGrs80Ellipsoid)));
var latitudeOld = 2 * Math.PI;
while (Math.Abs(latitude - latitudeOld) > 0.0000000000000001)
{
latitudeOld = latitude;
var transverseRadiusOfCurvature2 = Grs80SemiMajorAxis /
Math.Sqrt(1 - eccentricityOfGrs80Ellipsoid * latitudeOld.SinToPower(2));
latitude = Math.Atan2(grs80.Z + eccentricityOfGrs80Ellipsoid * transverseRadiusOfCurvature2 * Math.Sin(latitudeOld),
sphericalPolar);
}
return latitude;
}

private static double SinToPower(this double number, double power)
{
return Math.Pow(Math.Sin(number), power);
}

private static double TanToPower(this double number, double power)
{
return Math.Pow(Math.Tan(number), power);
}

private class LatitudeLongtitude
{
public double Longtitude { get; set; }
public double Latitude { get; set; }
}

private class Cartesian
{
public double X { get; set; }
public double Y { get; set; }
public double Z { get; set; }
}

public class WGS84
{
public double Longtitude { get; set; }
public double Latitude { get; set; }
}

public class OSGB36
{
public double Easting { get; set; }
public double Northing { get; set; }
}
}
}