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; } } } }