All posts by Hannah Fry

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

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

Converting Latitude and Longitude to British National grid

ORIGINALLY POSTED FEBRUARY 1, 2012

This code reads in a .csv file called LatLon, expecting two columns with headers – Latitude and Longitude (in WGS84, decimal form). If the script is run in the same directory as LatLon.csv, it will spit out a second file called LatLonandBNG.csv, with two additional columns: OSGB36 Eastings and Northings respectively.

For the inverse transform – OSGB36 to WGS84 – please refer to this post, where a you can find the relevant python script, and more details on the algorithms involved.

#This code converts lat lon (WGS84) to british national grid (OSBG36)
from scipy import *

import csv

def WGS84toOSGB36(lat, lon):

#First convert to radians
#These are on the wrong ellipsoid currently: GRS80. (Denoted by _1)
lat_1 = lat*pi/180
lon_1 = lon*pi/180

#Want to convert to the Airy 1830 ellipsoid, which has the following:
a_1, b_1 =6378137.000, 6356752.3141 #The GSR80 semi-major and semi-minor axes used for WGS84(m)
e2_1 = 1- (b_1*b_1)/(a_1*a_1) #The eccentricity of the GRS80 ellipsoid
nu_1 = a_1/sqrt(1-e2_1*sin(lat_1)**2)

#First convert to cartesian from spherical polar coordinates
H = 0 #Third spherical coord.
x_1 = (nu_1 + H)*cos(lat_1)*cos(lon_1)
y_1 = (nu_1+ H)*cos(lat_1)*sin(lon_1)
z_1 = ((1-e2_1)*nu_1 +H)*sin(lat_1)

#Perform Helmut transform (to go between GRS80 (_1) and Airy 1830 (_2))
s = 20.4894*10**-6 #The scale factor -1
tx, ty, tz = -446.448, 125.157, -542.060 #The translations along x,y,z axes respectively
rxs,rys,rzs = -0.1502, -0.2470, -0.8421#The rotations along x,y,z respectively, in seconds
rx, ry, rz = rxs*pi/(180*3600.), rys*pi/(180*3600.), rzs*pi/(180*3600.) #In radians
x_2 = tx + (1+s)*x_1 + (-rz)*y_1 + (ry)*z_1
y_2 = ty + (rz)*x_1+ (1+s)*y_1 + (-rx)*z_1
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
a, b = 6377563.396, 6356256.909 #The GSR80 semi-major and semi-minor axes used for WGS84(m)
e2 = 1- (b*b)/(a*a) #The eccentricity of the Airy 1830 ellipsoid
p = sqrt(x_2**2 + y_2**2)

#Lat is obtained by an iterative proceedure:
lat = arctan2(z_2,(p*(1-e2))) #Initial value
latold = 2*pi
while abs(lat - latold)>10**-16:
lat, latold = latold, lat
nu = a/sqrt(1-e2*sin(latold)**2)
lat = arctan2(z_2+e2*nu*sin(latold), p)

#Lon and height are then pretty easy
lon = arctan2(y_2,x_2)
H = p/cos(lat) - nu

#E, N are the British national grid coordinates - eastings and northings
F0 = 0.9996012717 #scale factor on the central meridian
lat0 = 49*pi/180#Latitude of true origin (radians)
lon0 = -2*pi/180#Longtitude of true origin and central meridian (radians)
N0, E0 = -100000, 400000#Northing & easting of true origin (m)
n = (a-b)/(a+b)

#meridional radius of curvature
rho = a*F0*(1-e2)*(1-e2*sin(lat)**2)**(-1.5)
eta2 = nu*F0/rho-1

M1 = (1 + n + (5/4)*n**2 + (5/4)*n**3) * (lat-lat0)
M2 = (3*n + 3*n**2 + (21/8)*n**3) * sin(lat-lat0) * cos(lat+lat0)
M3 = ((15/8)*n**2 + (15/8)*n**3) * sin(2*(lat-lat0)) * cos(2*(lat+lat0))
M4 = (35/24)*n**3 * sin(3*(lat-lat0)) * cos(3*(lat+lat0))

#meridional arc
M = b * F0 * (M1 - M2 + M3 - M4)

I = M + N0
II = nu*F0*sin(lat)*cos(lat)/2
III = nu*F0*sin(lat)*cos(lat)**3*(5- tan(lat)**2 + 9*eta2)/24
IIIA = nu*F0*sin(lat)*cos(lat)**5*(61- 58*tan(lat)**2 + tan(lat)**4)/720
IV = nu*F0*cos(lat)
V = nu*F0*cos(lat)**3*(nu/rho - tan(lat)**2)/6
VI = nu*F0*cos(lat)**5*(5 - 18* tan(lat)**2 + tan(lat)**4 + 14*eta2 - 58*eta2*tan(lat)**2)/120

N = I + II*(lon-lon0)**2 + III*(lon- lon0)**4 + IIIA*(lon-lon0)**6
E = E0 + IV*(lon-lon0) + V*(lon- lon0)**3 + VI*(lon- lon0)**5 

#Job's a good'n.
return E,N

#Read in from a file
LatLon = csv.reader(open('LatLon.csv', 'rU'), delimiter = ',')
LatLon.next()

#Get the output file ready
outputFile = open('LatLonandBNG.csv', 'wb')
output=csv.writer(outputFile,delimiter=',')
output.writerow(['Lat', 'Lon', 'E', 'N'])

#Loop through the data
for line in LatLon:
lat = line[0]
lon = line[1]
E, N = WGS84toOSGB36(float(lat), float(lon))
output.writerow([str(lat), str(lon), str(E), str(N)])

#Close the output file
outputFile.close()

Converting Latitude and Longitude to British National grid

ORIGINALLY POSTED FEBRUARY 1, 2012

This code reads in a .csv file called LatLon, expecting two columns with headers – Latitude and Longitude (in WGS84, decimal form). If the script is run in the same directory as LatLon.csv, it will spit out a second file called LatLonandBNG.csv, with two additional columns: OSGB36 Eastings and Northings respectively.

For the inverse transform – OSGB36 to WGS84 – please refer to this post, where a you can find the relevant python script, and more details on the algorithms involved.

#This code converts lat lon (WGS84) to british national grid (OSBG36)
from scipy import *

import csv

def WGS84toOSGB36(lat, lon):

#First convert to radians
#These are on the wrong ellipsoid currently: GRS80. (Denoted by _1)
lat_1 = lat*pi/180
lon_1 = lon*pi/180

#Want to convert to the Airy 1830 ellipsoid, which has the following:
a_1, b_1 =6378137.000, 6356752.3141 #The GSR80 semi-major and semi-minor axes used for WGS84(m)
e2_1 = 1- (b_1*b_1)/(a_1*a_1) #The eccentricity of the GRS80 ellipsoid
nu_1 = a_1/sqrt(1-e2_1*sin(lat_1)**2)

#First convert to cartesian from spherical polar coordinates
H = 0 #Third spherical coord.
x_1 = (nu_1 + H)*cos(lat_1)*cos(lon_1)
y_1 = (nu_1+ H)*cos(lat_1)*sin(lon_1)
z_1 = ((1-e2_1)*nu_1 +H)*sin(lat_1)

#Perform Helmut transform (to go between GRS80 (_1) and Airy 1830 (_2))
s = 20.4894*10**-6 #The scale factor -1
tx, ty, tz = -446.448, 125.157, -542.060 #The translations along x,y,z axes respectively
rxs,rys,rzs = -0.1502, -0.2470, -0.8421#The rotations along x,y,z respectively, in seconds
rx, ry, rz = rxs*pi/(180*3600.), rys*pi/(180*3600.), rzs*pi/(180*3600.) #In radians
x_2 = tx + (1+s)*x_1 + (-rz)*y_1 + (ry)*z_1
y_2 = ty + (rz)*x_1+ (1+s)*y_1 + (-rx)*z_1
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
a, b = 6377563.396, 6356256.909 #The GSR80 semi-major and semi-minor axes used for WGS84(m)
e2 = 1- (b*b)/(a*a) #The eccentricity of the Airy 1830 ellipsoid
p = sqrt(x_2**2 + y_2**2)

#Lat is obtained by an iterative proceedure:
lat = arctan2(z_2,(p*(1-e2))) #Initial value
latold = 2*pi
while abs(lat - latold)>10**-16:
lat, latold = latold, lat
nu = a/sqrt(1-e2*sin(latold)**2)
lat = arctan2(z_2+e2*nu*sin(latold), p)

#Lon and height are then pretty easy
lon = arctan2(y_2,x_2)
H = p/cos(lat) - nu

#E, N are the British national grid coordinates - eastings and northings
F0 = 0.9996012717 #scale factor on the central meridian
lat0 = 49*pi/180#Latitude of true origin (radians)
lon0 = -2*pi/180#Longtitude of true origin and central meridian (radians)
N0, E0 = -100000, 400000#Northing & easting of true origin (m)
n = (a-b)/(a+b)

#meridional radius of curvature
rho = a*F0*(1-e2)*(1-e2*sin(lat)**2)**(-1.5)
eta2 = nu*F0/rho-1

M1 = (1 + n + (5/4)*n**2 + (5/4)*n**3) * (lat-lat0)
M2 = (3*n + 3*n**2 + (21/8)*n**3) * sin(lat-lat0) * cos(lat+lat0)
M3 = ((15/8)*n**2 + (15/8)*n**3) * sin(2*(lat-lat0)) * cos(2*(lat+lat0))
M4 = (35/24)*n**3 * sin(3*(lat-lat0)) * cos(3*(lat+lat0))

#meridional arc
M = b * F0 * (M1 - M2 + M3 - M4)

I = M + N0
II = nu*F0*sin(lat)*cos(lat)/2
III = nu*F0*sin(lat)*cos(lat)**3*(5- tan(lat)**2 + 9*eta2)/24
IIIA = nu*F0*sin(lat)*cos(lat)**5*(61- 58*tan(lat)**2 + tan(lat)**4)/720
IV = nu*F0*cos(lat)
V = nu*F0*cos(lat)**3*(nu/rho - tan(lat)**2)/6
VI = nu*F0*cos(lat)**5*(5 - 18* tan(lat)**2 + tan(lat)**4 + 14*eta2 - 58*eta2*tan(lat)**2)/120

N = I + II*(lon-lon0)**2 + III*(lon- lon0)**4 + IIIA*(lon-lon0)**6
E = E0 + IV*(lon-lon0) + V*(lon- lon0)**3 + VI*(lon- lon0)**5 

#Job's a good'n.
return E,N

#Read in from a file
LatLon = csv.reader(open('LatLon.csv', 'rU'), delimiter = ',')
LatLon.next()

#Get the output file ready
outputFile = open('LatLonandBNG.csv', 'wb')
output=csv.writer(outputFile,delimiter=',')
output.writerow(['Lat', 'Lon', 'E', 'N'])

#Loop through the data
for line in LatLon:
lat = line[0]
lon = line[1]
E, N = WGS84toOSGB36(float(lat), float(lon))
output.writerow([str(lat), str(lon), str(E), str(N)])

#Close the output file
outputFile.close()

Converting British National Grid to Latitude and Longitude II

ORIGINALLY POSTED FEBRUARY 1, 2012

A few months ago, I wrote a python script to convert British National grid coordinates (OSGB36) to latitude and longitude (WGS84).  A fellow blogger Andrzej Bieniek very kindly pointed out that the algorithm was only accurate to around 100m because of an additional subtlety which I hadn’t taken into account.

I’ll have a quick bash at explaining the reasoning behind this difference,  but if you want to skip to the new version of the code (now accurate to 5m) I’ve posted it at the bottom of the page in all its delicious glory. Alternatively, if you are looking for a similarly accurate script applying the inverse transform – ie. WGS84 lat, lon to OSGB36 Eastings, Northings – you can find it here.

Because the Earth’s surface is not smooth, any attempt to represent a geographical area in two dimensions  – on paper, on screen etc – will naturally result in some distortions and discontinuities. One way to get round this, is to approximate the Earth using a ‘3D ellipse’ (or ‘ellipsoid’, if we’re being proper) try and fit it as best you can to the Earth’s squashed sphere shape, cross your fingers & not get too worried about ‘land’, ‘oceans’ and other such trivial inconveniences.

If you’re looking to consider large distances, or perhaps the whole globe at once – as is the case with WGS84 –  you’ll need to fit the ellipsoid to the entire Earth’s surface at once.  If, however, you’re only interested in mapping a small area of the globe –  as is the case with British National Grid  – it may be that a much smaller ellipse would work much better:

 

Of course, the picture above – taken from the Ordinance Survey website – is a greatly exaggerated illustration of this idea, but I think it demonstrates the concept pretty well.  This is exactly what I hadn’t realised when doing my previous post: the ellipsoids that are used for the two systems of British National grid and WGS84 latitude and longitude are ever so slightly different – hence why the results could be up to 100m out. British National grid is based on the Airy* 1830 ellipsoid, while latitude and longitude use the GSR80 (or WGS 1984).

The origins of these two, A1830 and GSR80 are in slightly different positions, their axes are a little rotated, and one is a tiny bit larger than the other – as in the picture below.

It’s not sufficient to just find  ‘lat/lon style’ coordinates from Eastings and Northings. You’ve also got to take the differences between the two ellipsoids before you can reproduce the latitude and longitude coordinates.

 

So, on top of the original post, to convert between the two ellipsoids you must scale, translate and rotate every point on the surface.

Actually, once you get down to it, the maths isn’t that bad – although the transform is a great deal simpler if you first convert to Cartesian coordinates, and that’s the only annoying bit. The full process is actually all contained on theOrdinance Survey guide  (just annoyingly not all in one place) and can be summarised as follows:

1) Use your Eastings and Northings to find lat/lon style coordinates on the Airy 1830 ellipsoid. (As in the previous post, or on pages 41&42 of the guide).

2) Convert those lat/lon coords to 3D Cartesian coordinates. (As on page 38 of the guide). Take the third spherical polar coordinate – height –   as zero in your calculations. (As though the point sits on the surface of the ellipsoid).

3) Use Helmert’s datum transform to convert between the two ellipsoids, and anchor your coordinates on the GSR80 system. (As on page 29 of the guide).

4) Your answer will be in Cartesian coordinates, so the final step is to convert back to spherical polars to get the final lat/lon. (As on page 39 of the guide).

All the constants you will need are also there in the guide – although you might have to hunt for them a little. Or just nick ‘em from my script below.

As before, the code reads in a .csv file called BNG, expecting two columns with headers – East and North. If the script is run in the same directory as BNG.csv, it will spit out a second file called BNGandLatLon.csv, with both sets of coordinate (lat and lon will be in decimal form). If you don’t have python, you can find instructions for the bits you need on another of my previous posts. Aren’t I nice to you?

*The same chap as Airy’s function, which I virtually have tattooed on the insides of my eyelids from my fluid dynamics days. Not that this is remotely interesting, I just think he’s a top fella is all.

#This code converts british national grid to lat lon

from math import sqrt, pi, sin, cos, tan, atan2 as arctan2
import csv

def OSGB36toWGS84(E,N):

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

    #Initialise the iterative variables
    lat,M = lat0, 0

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

    #transverse radius of curvature
    nu = a*F0/sqrt(1-e2*sin(lat)**2)

    #meridional radius of curvature
    rho = a*F0*(1-e2)*(1-e2*sin(lat)**2)**(-1.5)
    eta2 = nu/rho-1

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

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

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

    #Perform Helmut transform (to go between Airy 1830 (_1) and GRS80 (_2))
    s = -20.4894*10**-6 #The scale factor -1
    tx, ty, tz = 446.448, -125.157, + 542.060 #The translations along x,y,z axes respectively
    rxs,rys,rzs = 0.1502,  0.2470,  0.8421  #The rotations along x,y,z respectively, in seconds
    rx, ry, rz = rxs*pi/(180*3600.), rys*pi/(180*3600.), rzs*pi/(180*3600.) #In radians
    x_2 = tx + (1+s)*x_1 + (-rz)*y_1 + (ry)*z_1
    y_2 = ty + (rz)*x_1  + (1+s)*y_1 + (-rx)*z_1
    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    
    a_2, b_2 =6378137.000, 6356752.3141 #The GSR80 semi-major and semi-minor axes used for WGS84(m)
    e2_2 = 1- (b_2*b_2)/(a_2*a_2)   #The eccentricity of the GRS80 ellipsoid
    p = sqrt(x_2**2 + y_2**2)

    #Lat is obtained by an iterative proceedure:   
    lat = arctan2(z_2,(p*(1-e2_2))) #Initial value
    latold = 2*pi
    while abs(lat - latold)>10**-16: 
        lat, latold = latold, lat
        nu_2 = a_2/sqrt(1-e2_2*sin(latold)**2)
        lat = arctan2(z_2+e2_2*nu_2*sin(latold), p)

    #Lon and height are then pretty easy
    lon = arctan2(y_2,x_2)
    H = p/cos(lat) - nu_2

#Uncomment this line if you want to print the results
    #print [(lat-lat_1)*180/pi, (lon - lon_1)*180/pi]

    #Convert to degrees
    lat = lat*180/pi
    lon = lon*180/pi

    #Job's a good'n. 
    return lat, lon

#Read in from a file
BNG = csv.reader(open('BNG.csv', 'rU'), delimiter = ',')
BNG.next()

#Get the output file ready
outputFile = open('BNGandLatLon.csv', 'wb')
output=csv.writer(outputFile,delimiter=',')
output.writerow(['Lat', 'Lon', 'E', 'N'])

#Loop through the data
for E,N in BNG:
    lat, lon = OSGB36toWGS84(float(E), float(N))
    output.writerow([str(lat), str(lon), str(E), str(N)])
#Close the output file
outputFile.close()

Converting British National Grid to Latitude and Longitude II

ORIGINALLY POSTED FEBRUARY 1, 2012

A few months ago, I wrote a python script to convert British National grid coordinates (OSGB36) to latitude and longitude (WGS84).  A fellow blogger Andrzej Bieniek very kindly pointed out that the algorithm was only accurate to around 100m because of an additional subtlety which I hadn’t taken into account.

I’ll have a quick bash at explaining the reasoning behind this difference,  but if you want to skip to the new version of the code (now accurate to 5m) I’ve posted it at the bottom of the page in all its delicious glory. Alternatively, if you are looking for a similarly accurate script applying the inverse transform – ie. WGS84 lat, lon to OSGB36 Eastings, Northings – you can find it here.

Because the Earth’s surface is not smooth, any attempt to represent a geographical area in two dimensions  – on paper, on screen etc – will naturally result in some distortions and discontinuities. One way to get round this, is to approximate the Earth using a ‘3D ellipse’ (or ‘ellipsoid’, if we’re being proper) try and fit it as best you can to the Earth’s squashed sphere shape, cross your fingers & not get too worried about ‘land’, ‘oceans’ and other such trivial inconveniences.

If you’re looking to consider large distances, or perhaps the whole globe at once – as is the case with WGS84 –  you’ll need to fit the ellipsoid to the entire Earth’s surface at once.  If, however, you’re only interested in mapping a small area of the globe –  as is the case with British National Grid  – it may be that a much smaller ellipse would work much better:

 

Of course, the picture above – taken from the Ordinance Survey website – is a greatly exaggerated illustration of this idea, but I think it demonstrates the concept pretty well.  This is exactly what I hadn’t realised when doing my previous post: the ellipsoids that are used for the two systems of British National grid and WGS84 latitude and longitude are ever so slightly different – hence why the results could be up to 100m out. British National grid is based on the Airy* 1830 ellipsoid, while latitude and longitude use the GSR80 (or WGS 1984).

The origins of these two, A1830 and GSR80 are in slightly different positions, their axes are a little rotated, and one is a tiny bit larger than the other – as in the picture below.

It’s not sufficient to just find  ‘lat/lon style’ coordinates from Eastings and Northings. You’ve also got to take the differences between the two ellipsoids before you can reproduce the latitude and longitude coordinates.

 

So, on top of the original post, to convert between the two ellipsoids you must scale, translate and rotate every point on the surface.

Actually, once you get down to it, the maths isn’t that bad – although the transform is a great deal simpler if you first convert to Cartesian coordinates, and that’s the only annoying bit. The full process is actually all contained on theOrdinance Survey guide  (just annoyingly not all in one place) and can be summarised as follows:

1) Use your Eastings and Northings to find lat/lon style coordinates on the Airy 1830 ellipsoid. (As in the previous post, or on pages 41&42 of the guide).

2) Convert those lat/lon coords to 3D Cartesian coordinates. (As on page 38 of the guide). Take the third spherical polar coordinate – height –   as zero in your calculations. (As though the point sits on the surface of the ellipsoid).

3) Use Helmert’s datum transform to convert between the two ellipsoids, and anchor your coordinates on the GSR80 system. (As on page 29 of the guide).

4) Your answer will be in Cartesian coordinates, so the final step is to convert back to spherical polars to get the final lat/lon. (As on page 39 of the guide).

All the constants you will need are also there in the guide – although you might have to hunt for them a little. Or just nick ‘em from my script below.

As before, the code reads in a .csv file called BNG, expecting two columns with headers – East and North. If the script is run in the same directory as BNG.csv, it will spit out a second file called BNGandLatLon.csv, with both sets of coordinate (lat and lon will be in decimal form). If you don’t have python, you can find instructions for the bits you need on another of my previous posts. Aren’t I nice to you?

*The same chap as Airy’s function, which I virtually have tattooed on the insides of my eyelids from my fluid dynamics days. Not that this is remotely interesting, I just think he’s a top fella is all.

#This code converts british national grid to lat lon

from math import sqrt, pi, sin, cos, tan, atan2 as arctan2
import csv

def OSGB36toWGS84(E,N):

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

    #Initialise the iterative variables
    lat,M = lat0, 0

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

    #transverse radius of curvature
    nu = a*F0/sqrt(1-e2*sin(lat)**2)

    #meridional radius of curvature
    rho = a*F0*(1-e2)*(1-e2*sin(lat)**2)**(-1.5)
    eta2 = nu/rho-1

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

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

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

    #Perform Helmut transform (to go between Airy 1830 (_1) and GRS80 (_2))
    s = -20.4894*10**-6 #The scale factor -1
    tx, ty, tz = 446.448, -125.157, + 542.060 #The translations along x,y,z axes respectively
    rxs,rys,rzs = 0.1502,  0.2470,  0.8421  #The rotations along x,y,z respectively, in seconds
    rx, ry, rz = rxs*pi/(180*3600.), rys*pi/(180*3600.), rzs*pi/(180*3600.) #In radians
    x_2 = tx + (1+s)*x_1 + (-rz)*y_1 + (ry)*z_1
    y_2 = ty + (rz)*x_1  + (1+s)*y_1 + (-rx)*z_1
    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    
    a_2, b_2 =6378137.000, 6356752.3141 #The GSR80 semi-major and semi-minor axes used for WGS84(m)
    e2_2 = 1- (b_2*b_2)/(a_2*a_2)   #The eccentricity of the GRS80 ellipsoid
    p = sqrt(x_2**2 + y_2**2)

    #Lat is obtained by an iterative proceedure:   
    lat = arctan2(z_2,(p*(1-e2_2))) #Initial value
    latold = 2*pi
    while abs(lat - latold)>10**-16: 
        lat, latold = latold, lat
        nu_2 = a_2/sqrt(1-e2_2*sin(latold)**2)
        lat = arctan2(z_2+e2_2*nu_2*sin(latold), p)

    #Lon and height are then pretty easy
    lon = arctan2(y_2,x_2)
    H = p/cos(lat) - nu_2

#Uncomment this line if you want to print the results
    #print [(lat-lat_1)*180/pi, (lon - lon_1)*180/pi]

    #Convert to degrees
    lat = lat*180/pi
    lon = lon*180/pi

    #Job's a good'n. 
    return lat, lon

#Read in from a file
BNG = csv.reader(open('BNG.csv', 'rU'), delimiter = ',')
BNG.next()

#Get the output file ready
outputFile = open('BNGandLatLon.csv', 'wb')
output=csv.writer(outputFile,delimiter=',')
output.writerow(['Lat', 'Lon', 'E', 'N'])

#Loop through the data
for E,N in BNG:
    lat, lon = OSGB36toWGS84(float(E), float(N))
    output.writerow([str(lat), str(lon), str(E), str(N)])
#Close the output file
outputFile.close()

Converting British National Grid to Latitude and Longitude I

ORIGINALLY POSTED OCTOBER 10, 2011

 

EDIT: Andrzej Bieniek brought to my attention that this version is correct to 100m. For a more accurate script (accurate to 5m) see my new post.

I have recently started to deal with a lot of geographical data in my research and and begun to realise it is difficult to go to long without stumbling across the sticky world of map projections – something I knew almost nothing about a year ago.

There is a nice blog post explaining the background by James Cheshire which I shan’t attempt to reproduce, but explanations aside, I found myself today trying to convert a long list of British National grid coordinates into Latitude and Longitude.

Of course, if you’ve got a week spare, you can do these conversions using QGis or ArcGis, but all I was looking for was just a quick .csv to .csv converter  of one coordinate system to the other. After a fruitless half hour of googling, I decided it would just be quicker to write a python script myself to get the job done. And because I know I’ll use it again, here it is.

The code reads in a .csv file called BNG, expecting two columns with headers – East and North. If the script is run in the same directory as BNG.csv, it will spit out a second file called BNGandLatLon.csv, with both sets of coordinate. All the maths behind the conversion can be found on page 41 of a pdf written by the chaps at ordnance survey. Sorry about the lack of highlighting – I’ll sort it out when I have a spare moment.

#This code converts british national grid to lat lon

from scipy import *
import csv

def OSGB36toWGS84(E,N):

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

#Initialise the iterative variables
lat,M = lat0, 0

while N-N0-M >= 0.00001: #Accurate to 0.01mm
lat = (N-N0-M)/(a*F0) + lat;
M1 = (1 + n + (5/4)*n**2 + (5/4)*n**3) * (lat-lat0)
M2 = (3*n + 3*n**2 + (21/8)*n**3) * sin(lat-lat0) * cos(lat+lat0)
M3 = ((15/8)*n**2 + (15/8)*n**3) * sin(2*(lat-lat0)) * cos(2*(lat+lat0))
M4 = (35/24)*n**3 * sin(3*(lat-lat0)) * cos(3*(lat+lat0))
#meridional arc
M = b * F0 * (M1 - M2 + M3 - M4)

#transverse radius of curvature
nu = a*F0/sqrt(1-e2*sin(lat)**2)
#meridional radius of curvature
rho = a*F0*(1-e2)*(1-e2*sin(lat)**2)**(-1.5)
eta2 = nu/rho-1

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

lat = lat - VII*dE**2 + VIII*dE**4 - IX*dE**6
lon = lon0 + X*dE - XI*dE**3 + XII*dE**5 - XIIA*dE**7

#Convert to degrees
lat = lat*180/pi
lon = lon*180/pi

return lat, lon

#Read in from a file
BNG = csv.reader(open('BNG.csv', 'rU'), delimiter = ',')
BNG.next()

#Get the output file ready
outputFile = open('BNGandLatLon.csv', 'wb')
output=csv.writer(outputFile,delimiter=',')
output.writerow(['Lat', 'Lon', 'E', 'N'])

#Loop through the data
for E,N in BNG:
lat, lon = OSGB36toWGS84(float(E), float(N))
output.writerow([str(lat), str(lon), str(E), str(N)])
#Close the output file
outputFile.close()

Converting British National Grid to Latitude and Longitude I

ORIGINALLY POSTED OCTOBER 10, 2011

 

EDIT: Andrzej Bieniek brought to my attention that this version is correct to 100m. For a more accurate script (accurate to 5m) see my new post.

I have recently started to deal with a lot of geographical data in my research and and begun to realise it is difficult to go to long without stumbling across the sticky world of map projections – something I knew almost nothing about a year ago.

There is a nice blog post explaining the background by James Cheshire which I shan’t attempt to reproduce, but explanations aside, I found myself today trying to convert a long list of British National grid coordinates into Latitude and Longitude.

Of course, if you’ve got a week spare, you can do these conversions using QGis or ArcGis, but all I was looking for was just a quick .csv to .csv converter  of one coordinate system to the other. After a fruitless half hour of googling, I decided it would just be quicker to write a python script myself to get the job done. And because I know I’ll use it again, here it is.

The code reads in a .csv file called BNG, expecting two columns with headers – East and North. If the script is run in the same directory as BNG.csv, it will spit out a second file called BNGandLatLon.csv, with both sets of coordinate. All the maths behind the conversion can be found on page 41 of a pdf written by the chaps at ordnance survey. Sorry about the lack of highlighting – I’ll sort it out when I have a spare moment.

#This code converts british national grid to lat lon

from scipy import *
import csv

def OSGB36toWGS84(E,N):

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

#Initialise the iterative variables
lat,M = lat0, 0

while N-N0-M >= 0.00001: #Accurate to 0.01mm
lat = (N-N0-M)/(a*F0) + lat;
M1 = (1 + n + (5/4)*n**2 + (5/4)*n**3) * (lat-lat0)
M2 = (3*n + 3*n**2 + (21/8)*n**3) * sin(lat-lat0) * cos(lat+lat0)
M3 = ((15/8)*n**2 + (15/8)*n**3) * sin(2*(lat-lat0)) * cos(2*(lat+lat0))
M4 = (35/24)*n**3 * sin(3*(lat-lat0)) * cos(3*(lat+lat0))
#meridional arc
M = b * F0 * (M1 - M2 + M3 - M4)

#transverse radius of curvature
nu = a*F0/sqrt(1-e2*sin(lat)**2)
#meridional radius of curvature
rho = a*F0*(1-e2)*(1-e2*sin(lat)**2)**(-1.5)
eta2 = nu/rho-1

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

lat = lat - VII*dE**2 + VIII*dE**4 - IX*dE**6
lon = lon0 + X*dE - XI*dE**3 + XII*dE**5 - XIIA*dE**7

#Convert to degrees
lat = lat*180/pi
lon = lon*180/pi

return lat, lon

#Read in from a file
BNG = csv.reader(open('BNG.csv', 'rU'), delimiter = ',')
BNG.next()

#Get the output file ready
outputFile = open('BNGandLatLon.csv', 'wb')
output=csv.writer(outputFile,delimiter=',')
output.writerow(['Lat', 'Lon', 'E', 'N'])

#Loop through the data
for E,N in BNG:
lat, lon = OSGB36toWGS84(float(E), float(N))
output.writerow([str(lat), str(lon), str(E), str(N)])
#Close the output file
outputFile.close()

Bezier Curves

ORIGINALLY POSTED NOVEMBER 16, 2011

This week: Bezier curves, how to draw them in Python and particularly, how to decide where to put the control points.

A Bezier curve, is a special type of parametric curve used frequently in computer graphics. You may have seen them in Powerpoint, or using the pen tool in illustrator, but for my purpose, they are a lovely way to visualise flow paths on a map, like this:

This particular visualisation explores the London public transport network and was made by my rather awesome colleague Dr Ed Manley. It shows London's most popular destinations by origin, derived from a large dataset of morning peak Oyster Card trips.

The curves themselves are anchored at both ends (origin P0, destination P3) and get their shape from two control points (blue P1 and red P2), as in the following picture:

 

The curve itself is created from the fairly simple expression (everything in bold is a vector):

latex.png

The parameterisation t starts from zero, so that the curve begins at P0, and as it increases bends the line towards the two control points, first P1 and then P2. Finally, when t  is 1, the line finishes at P3.

It’s all very clever stuff, and there is lots and lots of literature available on the web about them. But – as far as I can see – no one ever tells you how to decide where to put your control points, and certainly doesn’t offer an algorithm to automatically do so.

But fear not people! I have written one. And ugly and fudgy as it is, it does the job.

The code below takes 6 inputs:

-origin points,

-destination points,

-a length for the blue line,

-a length for the red line (both as a fraction of the length of the green line)

-the angle between the blue and green lines

-the angle between the red and green lines (both in radians. I am a mathematician after all).

And it will spit out a totally delicious plot of your curves.

It’s set up to bend the curves clockwise (as in the image at the top of the page), which generally looks nicer when plotting more than one bezier curve.

Have a play with the last four inputs: the longer the length, the more the curve will deviate from the green line. The higher the angle, the steeper the curve will roll in to it’s destination. Indeed, choosing the angles to be on the other side of the green line will result in an s-shaped curve. Anyway – the idea is for you to hack out the bits of the code that are useful.

import scipy as sp
import pylab as plt

#### Inputs

#A list of P0's and P3's. Must be the same length
origins = [[1,0],[-0.5,sp.sqrt(3)/2], [-0.5,-sp.sqrt(3)/2]]
destinations = [[0,0],[0,0],[0,0]]

#The angle the control point will make with the green line
blue_angle = sp.pi/6
red_angle = sp.pi/4

#And the lengths of the lines (as a fraction of the length of the green one)
blue_len = 1./5
red_len = 1./3

### Workings

#Generate the figure
fig = plt.figure()
ax = fig.add_subplot(111)
ax.hold(True)

#Setup the parameterisation
t = sp.linspace(0,1,100)

for i in xrange(len(origins)):
#Read in the origin & destination points
POx,POy = origins[i][0], origins[i][1]
P3x,P3y = destinations[i][0], destinations[i][1]

#Add those to the axes
ax.plot(POx,POy, 'ob')
ax.plot(P3x,P3y, 'or')
ax.plot((POx,P3x),(POy,P3y), 'g')

#Work out r and theta (as if based at P3)
r = ((POx-P3x)**2 + (POy-P3y)**2)**0.5
theta = sp.arctan2((POy-P3y),(POx-P3x))

#Find the relevant angles for the control points
aO =theta + blue_angle+ sp.pi
aD = theta - red_angle

#Work out the control points
P1x, P1y = POx+ blue_len*r*sp.cos(aO), POy + blue_len*r*sp.sin(aO)
P2x, P2y = P3x+ red_len*r*sp.cos(aD), P3y + red_len*r*sp.sin(aD)

#Plot the control points and their vectors
ax.plot((P3x,P2x),(P3y,P2y), 'r')
ax.plot((POx,P1x),(POy,P1y), 'b')
ax.plot(P1x, P1y, 'ob')
ax.plot(P2x, P2y, 'or')

#Use the Bezier formula
Bx = (1-t)**3*POx + 3*(1-t)**2*t*P1x + 3*(1-t)*t**2*P2x + t**3*P3x
By = (1-t)**3*POy + 3*(1-t)**2*t*P1y + 3*(1-t)*t**2*P2y + t**3*P3y

#Plot the Bezier curve
ax.plot(Bx, By, 'k')

#Save it
plt.savefig('totally-awesome-bezier.png')
plt.show()
#Bosch.

Bezier Curves

ORIGINALLY POSTED NOVEMBER 16, 2011

This week: Bezier curves, how to draw them in Python and particularly, how to decide where to put the control points.

A Bezier curve, is a special type of parametric curve used frequently in computer graphics. You may have seen them in Powerpoint, or using the pen tool in illustrator, but for my purpose, they are a lovely way to visualise flow paths on a map, like this:

This particular visualisation explores the London public transport network and was made by my rather awesome colleague Dr Ed Manley. It shows London's most popular destinations by origin, derived from a large dataset of morning peak Oyster Card trips.

The curves themselves are anchored at both ends (origin P0, destination P3) and get their shape from two control points (blue P1 and red P2), as in the following picture:

 

The curve itself is created from the fairly simple expression (everything in bold is a vector):

latex.png

The parameterisation t starts from zero, so that the curve begins at P0, and as it increases bends the line towards the two control points, first P1 and then P2. Finally, when t  is 1, the line finishes at P3.

It’s all very clever stuff, and there is lots and lots of literature available on the web about them. But – as far as I can see – no one ever tells you how to decide where to put your control points, and certainly doesn’t offer an algorithm to automatically do so.

But fear not people! I have written one. And ugly and fudgy as it is, it does the job.

The code below takes 6 inputs:

-origin points,

-destination points,

-a length for the blue line,

-a length for the red line (both as a fraction of the length of the green line)

-the angle between the blue and green lines

-the angle between the red and green lines (both in radians. I am a mathematician after all).

And it will spit out a totally delicious plot of your curves.

It’s set up to bend the curves clockwise (as in the image at the top of the page), which generally looks nicer when plotting more than one bezier curve.

Have a play with the last four inputs: the longer the length, the more the curve will deviate from the green line. The higher the angle, the steeper the curve will roll in to it’s destination. Indeed, choosing the angles to be on the other side of the green line will result in an s-shaped curve. Anyway – the idea is for you to hack out the bits of the code that are useful.

import scipy as sp
import pylab as plt

#### Inputs

#A list of P0's and P3's. Must be the same length
origins = [[1,0],[-0.5,sp.sqrt(3)/2], [-0.5,-sp.sqrt(3)/2]]
destinations = [[0,0],[0,0],[0,0]]

#The angle the control point will make with the green line
blue_angle = sp.pi/6
red_angle = sp.pi/4

#And the lengths of the lines (as a fraction of the length of the green one)
blue_len = 1./5
red_len = 1./3

### Workings

#Generate the figure
fig = plt.figure()
ax = fig.add_subplot(111)
ax.hold(True)

#Setup the parameterisation
t = sp.linspace(0,1,100)

for i in xrange(len(origins)):
#Read in the origin & destination points
POx,POy = origins[i][0], origins[i][1]
P3x,P3y = destinations[i][0], destinations[i][1]

#Add those to the axes
ax.plot(POx,POy, 'ob')
ax.plot(P3x,P3y, 'or')
ax.plot((POx,P3x),(POy,P3y), 'g')

#Work out r and theta (as if based at P3)
r = ((POx-P3x)**2 + (POy-P3y)**2)**0.5
theta = sp.arctan2((POy-P3y),(POx-P3x))

#Find the relevant angles for the control points
aO =theta + blue_angle+ sp.pi
aD = theta - red_angle

#Work out the control points
P1x, P1y = POx+ blue_len*r*sp.cos(aO), POy + blue_len*r*sp.sin(aO)
P2x, P2y = P3x+ red_len*r*sp.cos(aD), P3y + red_len*r*sp.sin(aD)

#Plot the control points and their vectors
ax.plot((P3x,P2x),(P3y,P2y), 'r')
ax.plot((POx,P1x),(POy,P1y), 'b')
ax.plot(P1x, P1y, 'ob')
ax.plot(P2x, P2y, 'or')

#Use the Bezier formula
Bx = (1-t)**3*POx + 3*(1-t)**2*t*P1x + 3*(1-t)*t**2*P2x + t**3*P3x
By = (1-t)**3*POy + 3*(1-t)**2*t*P1y + 3*(1-t)*t**2*P2y + t**3*P3y

#Plot the Bezier curve
ax.plot(Bx, By, 'k')

#Save it
plt.savefig('totally-awesome-bezier.png')
plt.show()
#Bosch.

Counting Crowds

Screen Shot 2013-07-30 at 13.55.23

This month has seen some dramatic developments in Egypt, with widespread protests, Morsi deposed from power – and more recently – outbreaks of violence across the country.

In the early days of the events, a nationwide protest was organised on the 30th of June to call for Morsi’s resignation. Huge numbers of people turned to the streets, with reports that up to 30 million people were involved – leading to the media claiming this was the largest protest in history.

Motivated by this headline, BBC Radio 4′s More or Less looked to investigate these numbers, and find out a bit more about how big crowd numbers are calculated. I was interviewed for a section on the show which you can find on iPlayer: http://www.bbc.co.uk/iplayer/episode/p01c2m3k/More_or_Less_Egypt_Biggest_Protest_in_History/  or download the podcast from the More or Less website:  http://www.bbc.co.uk/podcasts/series/moreorless

In addition, we filmed a short explanation of how to count crowds for the BBC website: : http://www.bbc.co.uk/news/world-23314873 and Ruth Alexander wrote a BBC magazine piece to accompany the work: http://www.bbc.co.uk/news/magazine-23312656

How to install Networkx on a mac

1) Download the link “networkx-1.6.zip” from here: http://networkx.lanl.gov/download/networkx/

2) Open your downloads folder and click on the file you just downloaded. Hopefully it should create a second, uncompressed file (also in your downloads folder) called “networkx-1.6″
3) Press Command + Space to open spotlight, and search for Terminal. Click on the top result and a white window with some black text should appear.
4) You now need to navigate your terminal to the downloads folder. It might be that all you need to do is type
“cd Downloads” BUT it depends on how your system is set up.
Essentially this terminal is a text version of the folders in finder.
-To see what folders are in your current location you type “ls”.
-To move into a new folder within your current directory (say a folder called “Hannah”) you type “cd Hannah”.
-To go up a folder you type “cd ..”
Move around until you find the downloads folder.
5) From within the downloads folder type “cd networkx-1.6″ into the terminal window
6) Also in the terminal window, type “sudo python setup.py install”
7) Type in your password (don’t worry that nothing appears on the screen)
8) Lots of text should come up. Have a glance through and see if an error comes up. If it does, send me a screen shot, if it doesn’t, go on to step 9.
9) Open up python and type in “import networkx” If no red text pops up then hurrah! You’ve installed everything correctly.
This might be helpful too:
http://networkx.github.io/documentation/latest/tutorial/tutorial.html

 

The return of the Tuesday Teaser

TuesdayTeaserI’ve missed puzzles. I’ve been too busy recently ‘working’ and ‘being a grown up’ to play with any. But this week CASA’s @frogo  has lined me up with a couple of absolute stonkers and they’ve re-ignited my desire to bring back the Tuesday Teasers. I know, right? Brilliant news. <Whoo! Yeah!>

And so.. without further ado. Here’s one from @frogo:

There was once an evil wizard with a weird and slightly unsavoury fetish for jewel encrusted goblins. Let’s call him Harry.

Harry liked to kidnap goblins and lock them up in his underground dungeon. As he did so, he magically implanted a precious gem in each goblin’s forehead – choosing either a ruby or an emerald each time, but refusing to tell the goblin which gem he had implanted. All they knew was that at least one of their fellow prisoners had a ruby, and at least one had an emerald. 

One day, Harry’s friend Germione – who was generally quite a decent soul – made him feel all guilty about being mean to the poor goblins, and persuaded him to give them a chance at freedom.

The next morning before breakfast – and every morning thereafter –  Harry would line up the goblins and ask all those with the rubys in their forehead to step forward – if exactly all of the goblins with a ruby stepped forward, while every goblin with an emerald remained in line, all the goblins would be freed and Harry would promise to never bother them again.

If they wanted, the goblins could all choose to remain in their line and not participate in Harry’s twisted little game & come to no harm.  But, if the incorrect number of rubied goblins stepped forward on any given morning, Harry would massacre every last goblin in a bloody, sickening and frenzied attack. And that Harry has a pretty nasty reputation.

Complete silence is required in the dungeon and the goblins are not allowed to talk to each other nor discuss strategies nor communicate in any way, there are no mirrors and each goblin has no way of knowing which gem they have in their heads without using bad-ass logic to work it out.

So.. how did the goblins win their freedom?

Deprivation, youth services and the London riots

Screen shot 2013-02-21 at 16.58.37

This video came about as a splinter project of some work I’ve been doing at UCL, trying to understand the 2011 London riots from a mathematical perspective.

While we were looking into the data behind the events, it was hard not to notice the striking links between the areas of the city that most rioters came from and the deprivation of those communities. This link was not the focus of our academic project, but I felt it was important not to just leave it as a one-line comment in a journal article. So I tried to investigate further by speaking to charities, youth workers and young people across the capital.

Once the riots had died down, many were very quick to demonise the groups of young people involved, but after taking the time to speak to young people from the poorest areas of the city I heard stories that shocked, appalled and surprised me – even as a Londoner myself. Perhaps the riots were not the result of innate criminal instincts lying dormant in certain groups of society, but rather a manifestation of deeper underlying problems in our cities.

I thought the voices of these young people and the often desperate situations they face deserved to be heard, and so obtained some funding through UCL’s Focus on the Positive to make this film. It turned out to be harder than I thought to get young people to open up on film but I hope that, regardless, this short piece might offer some insight into another way of looking at things.

I should point out that the links to deprivation and government cuts discussed here are not the subject or focus of the Scientific Reports paper ‘A mathematical model of the London Riots and their policing‘. The work described in this video has not been peer reviewed and does not represent the views of UCL or the other authors.

Who says complex networks aren’t funny?

Ok. Well maybe they aren’t. But I recently gave it a good bash at a stand-up comedy routine for UCL’s brightclub.

The night is a sort of academic variety night, where researchers from across the UK are given five minute slots to entertain a 600 seater theatre with a light-hearted look at their work. And here is my attempt – officially: looking at immunisation strategies on scale-free networks; unofficially: mucking about onstage with some good old fashioned ginger jokes.

WARNING. Contains some naughty words. And a hideous powerpoint fail at the start. Action begins from around 0:50.

TEDxUCL


UCL recently hosted their first ever TEDx event, and to my amazement (probably after seeing some seriously dodgy maths themed stand-up I did) they invited me to give a talk.

So if you have a spare ten minutes and fancy watching me ramble on about snooker balls, the solar system, leopards and how they all relate to complexity science then you’ve come to the right place.

Alternatively, use this link.. Or find out more about the event here.

EDIT: I’m totally mind bogglingly  amazed and thrilled-thrilled -thrilled that the lovely chaps at TED decided to promote my video as their talk of the day on the 4th August 2012.

Here is the link to the TED page: http://www.ted.com/talks/hannah_fry_is_life_really_that_complex.html

and to my bio on their site: http://www.ted.com/speakers/hannah_fry.html

Focus on the positive

focus

On Monday, UCL public engagement played host to a new type of event, funded by EPSRC and titled Focus On The Positive.

The idea is beautifully simple – six academics pitch a solution to a real world problem to an audience, who then vote for whichever cause they deem the most worthy. The winner then has £1000 to spend towards funding their solution.

While working on the London riots model with Toby Davies (which you can read more about in this article from Wired magazine) we realised how striking the links were between rioter involvement and deprivation and that the kids involved tended to come from the London boroughs worst hit by recent government cuts.

We thought the focus on the positive event would be the ideal opportunity to raise awareness about these links and the issues faced by young people in London from deprived backgrounds  - in particular given the recent cuts in funding to essential youth services. We decided our best chance at a lasting legacy from the funding would be to make a short film, and pitched our idea to the audience at the event on Monday.

I am incredibly honoured to say that our project won the funding on the night, and we are now in the processes of pre-production on our film, which I will be sure to post here when it’s ready.

In the meantime, the next Focus on the positive event will be held in Bloomsbury theatre on the 19th of June. You can find out more about the night here, buy tickets here, or contact Hilary Jackson more information

Tuesday Teaser 20th March

This week, a politics question:

Since WW2, how many British Prime Ministers first got the job without winning a general election?

As ever, feel free to gloat and post answers and guesses in the comments section. But only after the snake, mind. We don’t want to ruin it for everyone else.

Tuesday Teaser 13th March

Sorry for the lack of Tuesday Teaser last week, it’s because I was too busy living it up here. Back in full swing now, I promise. And to show you that I mean it, this week I’m pulling out one of my absolute favourite questions: 

There is only one country in the world whose name ends in the letter ‘k’. What is it? 

Please feel free to post answers and gloating remarks in the comments section below. And as ever, so as to keep up this delightful atmosphere of fun, I’ve included the silly snake to entertain you while you’re scrolling.

 

And the correct answer is…. Denmark.

The UK doesn’t count, you douches.

Tuesday Teaser 28th Feb

I promised I would try and avoid general knowledge teasers, but I just really like this question – so I hope that you’ll forgive me, just this once.

There are currently 7 Communist countries in the world. Name them.

And for a bonus point.. Which of these 7 are also democracies? 

As ever, feel free to gloat/post your answers after the snake..

 

 

Time for the answers, which you can read more about here.

The current communist one-party states are:

The Republic of Cuba (since 1948)

The Peoples Republic of China (since 1949)

The Socialist Republic of Vietnam (Communist since 1954, one-party state since 1976)

Democratic People’s Republic of Korea (since 1948)

Lao People’s Democratic Republic (since 1975)

And the democratic states with communist government:

Cyprus, (Progressive Party of Working People since February 2008),

Nepal (Unified Communist Party of Nepal (Maoist) since April 2008).

 

 

Tuesday Teaser 21st Feb

This week a word puzzle:

Which plural noun in the English language becomes singular when you add an s?

As ever, please feel free to gloat in the comments section below.. Answers after the snake

 

 

Bit of a controversial one this one. So far I have found 4 nouns which fit the above description. They are:

Princes/Princess

Bras/Brass

Millionaires/Millionaress

Billionaire/Billionairess

Any others?

 

Tuesday Teaser 14th Feb

This week’s teaser is about the London tube network:

Which is the only London Underground station whose name does not include any of the letters in the word ‘Mackerel’?

With this question, you either get it or you don’t, so it’s a bit tricky to allow solutions to be posted. Feel free to gloat in the comments section though if you know the answer.

UPDATE: Answer is now posted in the answers section below. Well done to F. Dans for getting it right.

Tuesday Teaser 7th Feb

This week’s teaser is a geography one..

Which state of the USA is the most:

i) Northerly

ii) Easterly

iii) Westerly

iv) Southerly. 

Feel free to post your solutions at the bottom of the page.  I’ve now included a lovely drawing of a little snake to amuse you while you’re scrolling.

UPDATE: The answers are now up, and can be found after the snake.

Time for the answers. First up, F. Dans was right that Attu is the most westerly island of Alaska, but it actually lies seven degrees west of the 180° longitude line, placing it in the Eastern Hemisphere. Hence, it is another member of the Aleutian Islands: Semisopochnoi, which can claim to be the most Easterly point of the US, being 14 minutes West of the 180° longitude line.

You can read more about it here and here on good old Wikipedia.

Regardless, that leaves Alaska as the most Northerly, Easterly, and Westerly state, with Hawaii rounding off the answers as the most Southerly.

Well done to Owen Boswarva for being the first one to post the right answer. You win my unbounded respect.

Welcome to Tuesday Teasers


For the past 8 years or so, I’ve been a member of a pub quiz team with some good friends of mine from my undergraduate days.

We’ve  had differing degrees of success (inversely correlated with the strength of our arch-rivals More Cool Mafia, as it goes) but in that time, I’ve come across some absolutely stonking good quiz questions – the type that most people know the answer to, but are still really tricky to get.

I thought it might be fun to post some of these questions here in a new regular feature (well, we’ll see how long it lasts) called Tuesday Teasers.

And so, to kick us off:

There are 10 parts of the human body, whose names in common usage are spelt with three letters each. What are they?

Slang words – like ‘bum’, ‘tum’ etc –  and medical terms don’t count. Don’t get all smart on me now.

UPDATE: The answers are now posted  at the bottom of the page. So you don’t get bored while you’re scrolling, I’ve included a picture of a lovely little snake.

And the answers are:

Arm, Leg, Toe, Eye, Ear, Lip, Hip, Rib, Gum and Jaw

Well done if you got it right!

Basic python syntax

A friend of mine has just started using Python, and asked me for a quick run down of the basic syntax. Given that, wherever possible, I try to live my life according to the above xkcd comic, I figured it would be worth my while typing out a short tutorial for her. Of course, there are thousands of similar (and probably far better) scripts available on the internet. But this seemed as good a place as any to use as a storage facility.

#Python has a few differerent types of built in objects.

#Integers
itg = 3

#Floats
flt = 3.

#If you divide an integer by another, you will get the floor of the answer
print 'Dividing integers'
print itg/2

#But you can change the object type:
print 'Converting to floats'
print float(itg)/2

#You can also convert to integers
print int(flt)

#And strings
print str(flt)

#Actually, when printing, there is an easier way to do this using the percentage symbol
print 'Using the percentage symbol %f' %flt

#Here %f means - find the float
#Likewise %s means - find the string
print 'Another example %s' %'Here'

#this is a good way to do print statements as you can use more than one simultaneously
print '.2d means to two significant figures like this %.2d or include a float %f' %(flt, flt)

#You can pull out bits of strings using the index labels, which start from zero.
name = 'Minnie'

#Then the command:
print name[0]
#Should print 'M'

#Working backwards, the labels start with -1:
print name[-1]
#Will print 'e'

#This labelling system works for lists too - another important python object
a = [] #An empty list
b = [1, 4, 'hello']#A list with three items.

#Then:
print b[-1]
#Should print 'hello'

#You can add items on to the end of lists like this:
b.append('new item')

#See?
print b

#Which is especially usefull when you're creating lists using loops
for n in name: #Will loop through all the letters in name
    a.append(n)#And add them to the list a as individual items

print a
#See?

#Another thing to note is that no end statements are used in python. Everything is based on indentation

#There are a few other things you can do with lists, because they are objects. Type a. and then press tab in the shell window to
#see the possible commands. they include a.pop, a.index and so on. 

#You can also put for loops inside lists like this - which is the same as above, just written on one line
a2 = [n for n in name]

#Lists are also super useful because they can contain objects themselves.

#To explain this, we need to use an example of a class
class Girl():
    def __init__(self, name, age):
        self.name  = name
        self.age = age

#The init function, within a class, initialises the instance of the class. For example
g = Girl('Minnie', 30)

#Makes a girl object with name minnie and age 30. Because of the self tag, you can now type
print g.name

#to get Minnie or
print g.age

#To get 30.

#Back to lists, you don't even need to give a handle to each instance of the girl class. Instead, you could make a list, 'girls'
girls = []
girls.append(g)
girls.append(Girl('Hannah', 27))

#Now the first item in the girls list takes the place of your g handle above:
print girls[0].name

#Same with the last item
print girls[-1].name

#Infact typing girls[1]. and then tab in the shell will give you the options which we defined in Girl: age and name.

#Geddit?

#There's one thing you need to be careful of with lists/python in general which is when setting things equal to each other.
list1 = [1]
list2 = list1

#Now, list2 isn't *actually* equal to list1, it just points to the same address on the computer.
print list2

#Now if you change list1:
list1[0] = 5

#List2 will automatically update itself
print list2

#You should be careful of this when coding as it can catch you out (it certainly has me in the past). Actually though,
#Although it initally seems stupid, it's actually a really useful (and quick) feature.

#Say you're working out a function f which depends on itself at the previous time step
#f(t+1) = f(f(t))

#You can code it like this:
#fnew = f**2 + 3 (or whatever the dependence is)
#Then when you're done, you can swap the pointers over
#fnew, f = f, fnew

#The other important built in object in python are dictionaries
D = {} #An empty dictionary

#You could use this instead of the girls object above
girlsD = {'Hannah':27, 'Minnie':30}

#Or to add a new item:
girlsD['Nenna'] = 26

#These are especailly useful when doing stuff with networks, but you should bear in mind that dictionaries DO NOT preserve order (lists do).

#Again type girlsD. and then tab in the shell to see the options available to you.

#That pretty much brings us to the end of what you can do without importing any other packages, so we'll move on to the most useful python
#package (as far as mathematical programming is concerned) Scipy.

#It's custom to import it with a shortcut handle
import scipy as sp

#Then anything that's in scipy can be called using the sp handle
sp.exp(5.)

#Actually, you can call the handle anything you like
import scipy as anything
anything.exp(5.)

#Or, not use one
import scipy
scipy.exp(5.)

#Or import everything, so you don't need to use a handle at all
from scipy import *
exp(5.)

#This is not that efficient though tbh, so probs best avoided.

#Now scipy has lots of cool stuff that you can see online, but the most important feature is the arrays. With lists you see, you can't add them
#Or do any arithmetic with them
a = [1,2,3]
b = [3,2,1]

#Now, if you typed a+b, it would produce an error.

#If a and b were arrays though:
a = sp.array([1,2,3])
b = sp.array([3,2,1])

#Now a+b would work.
print a+b

#and a*b gives element-wise multiplication
print a*b

#Geddit?

#a and b are 1d arrays, but you can also make vectors and matrices with sp.array:
row= sp.array([[1,2,3]])

#row is a 1x3 array
print sp.shape(row)

#likewise
col = sp.array([[1],[2],[3]])

#Is a 3x1
print sp.shape(col)

#Now, you could do matrix multiplication with these bad boys:
sp.dot(row, col)

#Which isn't the nicest notation, but it'll do.

#There are two other sp things that are quite useful:
xvals = sp.linspace(0, 100, 200)

#Auto generates an array from 0 to 100 with 200 points

#The other:
yvals = sp.arange(0, 10, 0.05)

#Lets you define your stepsize instead.

#Now, lets say we wanted to plot these, we'd need the plotting library
import pylab as plt

#Guess what the command is to plot something?
plt.plot(xvals, yvals)

#And if you want to save the figure?
plt.savefig('test')

#Phwoar. Got to love python.

PhD Thesis

I submitted my PhD thesis in January 2011, and was awarded the degree after my viva in March.

You can see the thesis in in full by clicking on the title page below: and if I do say so myself, it’s a delightful romp through the two-dimensional Navier-Stokes equations. I don’t want to give away the ending, but it should certainly please all you high Reynolds number two-phase flow fans out there.

 

Getting started with Python

This guide should give you a good overview of the packages you’ll need to download and install to start scientific programming using Python.

First, Python itself.  Python 2.7 is the latest version, and it works with the majority of toolkits you’ll need.

So, here it is:

http://python.org/download/releases/2.7.2/

You’ll also need the numerical extensions to python, NumPy (which lets you use arrays and matrices) and SciPy (which has lots of useful functions for integration, linear algebra, Fourier analysis and so on).

Luckily there are ready-made installers of both of these for both mac and windows:

Those for NumPy can be found here at
(For direct download of the 2.7 compatible version, you can use this link for mac or this one for windows)
.
And SciPy is  here (Or you can use the direct downloads of the 2.6 compatible version  for a mac or for windows)
.
If you want to do any plotting or visualisations you’ll need Matplotlib too. It’s very similar to matlab in its syntax and an overall delight to use. Find it here (Or directly download the 2.6 compatible version for windows and mac ).
.
NetworkX is another handy plugin for dealing with everything to do with networks. The easiest possible way to set it up is with another piece of software called easy install (documentation of which can be found here). I have a feeling (don’t quote me) that this comes pre-loaded on the mac,  but for windows, you can access the download page here or just download it directly.
.
Once that’s done, open up a terminal window (start>run>”cmd” for windows) and use the command ”easy_install networkx”. ‘Tis as simple as that.
.
Since my work is in a global dynamics project, I often have to produce plots overlaid on maps. For that matplotlib has a handy little toolbox called Basemap. The sourceforge page includes a lovely windows installer, unfortunately though, there’s no such thing for mac users. I have managed to set up basemap on my mac, but cannot for the life of me remember how. I know it took a long time to work it out, and I know I used this guide, which was about as useful as a chocolate teapot, but other than that you’re on your own I’m afraid.

1 + 1 = 2

A couple of days ago I saw a post on twitter about the proof that 1 + 1 = 2 by Alfred North Whitehead and Bertrand Russell. It appears in their 1910, three volume work, Principia Mathematica.
Despite the obvious importance of their work, I can definitely see the funny side. (Especially that it took them until page 379 to get this far, and they haven’t even defined arithmetic addition yet). But it does raise what I think is an interesting question of when to accept things as obvious, and when some things that should be obvious turn out to be quite counter intuitive.

One brilliant example was given in a talk by UCL’s own Jack Grahl in July. Consider the infinite sum:

 S = 1 - 1 + 1 - 1 + 1 - 1 + \cdots

and ask the question, what does S add up to? Obvious right?

First rewrite S, so that all the +1‘s and -1‘s are grouped together in pairs:
 S = (1 - 1) + (1 - 1) + (1 - 1) + \cdots
Each of these pairs cancel to zero:
 S = (0) + (0) + (0) + \cdots
So, the answer is
 S = 0

Done!. Although, not quite, because the choice of pairing could be done in a different way, leaving a +1 at the front:
 S = 1 + (-1 + 1) +(-1 + 1) + (-1 + 1) + \cdots
The terms in the brackets still cancel, so that:
 S = 1 + (0) + (0) + \cdots
leaving
 S = 1
eh??!?

And you can do this another way, and get a third answer. First writing S out twice:
 S = 1 - 1 + 1 - 1 + 1 - 1 + \cdots
 S = 0 + 1 - 1 + 1 - 1 + 1 - 1 + \cdots
(the second version is just the same, with a zero in the front to act as a kind of spacer).
Adding these two together. Inside each brackets is the top term, plus the bottom term:
 (S+S) = (1+0)+(- 1+1) + (1 - 1) + (-1+1) + \cdots
but again alot of this cancels, leaving:
 2 S = 1 + 0 + 0 + 0 + \cdots
Which can be solved to show
 S = \frac{1}{2}

So one seemingly simple problem has three mathematically legitimate results.

Counter-intuitively perhaps, given that it is a fraction made from the sum of an infinite series of integers, it is the third result of S = \frac{1}{2} which is generally accepted as ‘correct’. However, this answer seems positively sane when compared to the similar problem of determining the sum of the infinite series:
 1 + 2 + 3 + 4 + 5 + \cdots
Each term gets successively larger, so it should be obvious that if you continue adding infinitely many terms, you end up with an infinitely large sum, right? Wrong.

In a letter which the then unknown and relatively uneducated Ramanujan wrote to Hardy in 1913, he demonstrated his proof that the sum of the above series was.. (wait for it) -\frac{1}{12} .

As it turns out, this answer is absolutely correct, and can be verified by the Riemann Zeta function. Anecdote has it that that was when Hardy knew that Ramanujan was the real thing – if he was just barking mad and trying to pull a fast one, he would surely have come up with something much more sensible.

Mixing it up.


 

 

 

 

 

 

 

 

Following on from the storming success of my last puzzle post (I don’t like to brag, but I believe it had over 31 views..) I thought I’d pop up another. This one is potentially a little more involved – it took me a couple of days and some some fairly heavy maths before I was sure of the answer. However, knowing the solution, there is a simple logical argument involving only two physical concepts which I’m sure many people much cleverer than me could spot straight away.

The puzzle is as follows:

A vase is part filled with water, and part with oil, as in the first picture. The two are then mixed together, as in the second.

What happens to the pressure at the bottom of the vase if it is:
a) Cylindrical
b) Shaped as in the picture above.

And for the less fluid-dynamically minded, the two things you need to know to solve the problem are:

1) Pressure is force per unit area. In this case, pressure at a given point can be considered analogous to the weight of the fluid above that point.
2) Oil is lighter than water.

PS. All viscosity, intertia, vorticity, surface tension, compressibility, imiscibility, blah blah blah effects are ignored. Don’t get all smart on me now.

A little puzzle: the solution

I ruddy love this puzzle. (Find the original post here). Have had great fun tormenting friends and colleagues with it over the last few days.

To begin with, we – like Jones – are looking for four distinct numbers which add up to 18.

1,2,3,4
1,2,3,5
1,2,3,6…

and so on.

But Jones has an extra bit of information to us, in that he knows the door number of Prof White’s house.

The key thing to notice however, is that this still isn’t enough for Jones to determine a unique solution to the problem. Therefore, Prof White’s house number must be one which you can factorise into four numbers, which add up to less than 18, in more than one way.

This knocks off almost every possible house number apart from 60, for which there are two possible solutions which fit the bill:

1,2,5,6
1,3,4,5

And 120, for which there are three:

1,4,5,6
1,3,5,8
2,3,4,5

Now, Jones asks how many children are in the Black (smallest) family, and based on Prof. White’s answer he knows the solution to the problem.

In that case, Prof. White must have answered “2″, else Jones would still not have enough information (regardless of whether the house number was 60 or 120).  So the answer must be 2,3,4,5!

Phwoar. Don’t ya just love a good bit of maths?