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 the Ordinance 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()

Nice article – well done Hannah

Thanks for example. Python formatting is messed up – perhaps post download link?

Hi Richard – thanks for pointing that out. Should be sorted now, but otherwise, here is the download link:

http://dl.dropbox.com/u/20570844/OSGB36toWGS84.py

Thank you for this Hannah, this really helped to get the postcode search up and running on FlatmateRooms. Ordnance Survey only offers postcode geocoding in british national grid. It’s much better than nothing, but I was very surprised what they ask developers to figure out. In order to use postcode data with Bing or Google Maps API every web developer needs it in WGS84. My previous attempts at conversion ended up consistently several streets off, which is unacceptable for local navigation.

Hint: email them and ask for a link, it’s well deserved.

Python suggestions and legibility:

* double import of scipy on the top of the script

* import only the functions needed from scipy, not * as it would help me and perhaps others to understand what is stdlib and what is coming from scipy

* uppercase variable names are used for global variables and camelCase variable names are not very PEP8, spaces after a hash in a comment, also variables named N, n and nu in the same function make it a bit scary, but I really apprecciate the conciseness of the solution despite complexity of the task at hand.

Really good things on this blog.

Thanks Frank – all good points. I’ll see if I can tidy it up a bit. Appreciate the comments!

If anyone else is having problems with scipy dependency (Fortran compiler required on OS X). You can just use stdlib math module. Instead of

“from scipy import *” use

“from math import sqrt, pi, sin, cos, tan, atan2 as arctan2″ and all good.

Hello,

First of all thanks for writing this.

… I’m having troubles getting the program to work.

I have the BNG reference TL98942058, what do I need to put into the BNG.csv file for it to run and give me the Long/Lat co-ordinates?

Josh

Hi Josh,

I think that the TL is a grid reference to a particular 100km square in the British isles. After that, I would assume the first 4 digits are Eastings and the last four are Northings. See this link on the ordinance survey website, which appears to have been written for children:

http://www.ordnancesurvey.co.uk/oswebsite/gi/nationalgrid/nationalgrid.pdf

Unfortunately the BNG.csv file above is based on a grid with origin somewhere near the Isles of Scilly, whereas your reference appears to have its own origin within the TL grid square. You’d need to first convert between the two before the above script can be used I think..

Sorry I can’t be more help!

Hi

Many thanks Hannah for that script. I was looking for something like that since some time back.

Josh: to convert those letters at the beginning of the BNG coordinates into numbers so as to using Hannah’s function you have to do the following:

That couple of letter identifies a certain 100km square within the British national grid. They refer respectively to eastings and northings. Just count the number 100km squares you pass over from the bottom left corner to the one you are interested in. E.g. TL is 5 left (east) and 2 up (north). So the function’s input should be 59894 (eastings) and 22058 (northings).

I know it’s too late but actually it should be 598940 (eastings) and 220580 (northings).

Hannah,

Thank you very much for the code. Great stuff – it works brilliantly and has saved me loads of time coding and testing the conversion.

I noticed that the closed loop differences (errors?) between the conversion BNG->WGS84 and back again depend on latitude. I see upto a metre difference in the back and forward test locations at the most southery latitudes. Is this due to the iteration in BNG-> WGS or some other effect that is latitide dependent? It’s not an issue (for my purposes) but I wondered if you noticed it?

Thanks!

Ciaran

Hannah,

This is too good. I was using pyproj bindings but I was struggling to include them in compiled exe (using py2exe).

This saves my life

Thanks!!

I don’t fully understand all the above but I’m impressed! Will this convert BNG to lat/long in bulk or one at a time? I need to convert loads of points and doing them individually is making my wrist ache. Many thanks in advance for any advice.

Yes! It will do a whole list all at once. Just make an excel spreadsheet with your BNG coordinates in it: East in the first column and North in the second. Include “East” and “North” as headers at the top of those columns and then save the sheet as “BNG.csv”. Download this python script and save it in the same folder as BNG.csv. Run the python script and kaboom! You should have a second file appear called BNGandLatLon.csv which will have both sets of coordinates inside.

Hope that makes some sense, but feel free to come back if you need more help!

Many thanks Hannah, that seems easy enough. The nearest I got to a python script previously was quoting The Life of Brian so the simple instructions are appreciated!

Wow! You’re too clever!!!!

Thank you, this code works well. But you need to change

lat, latold = latold, lat

to

latold = lat

otherwise it is and endless loop.

Hi Kamel! If I remember rightly:

lat, latold = latold, lat

just swaps the pointers over for the two variables. Shouldn’t cause any problems with infinite loops..

This is excellent! Thanks!

Hannah, many thanks for publishing this.

Just one problem that I can see: in the first while loop you calculate fractions (5/4), (21/8), (15/8), (35/24). Without decimal points anywhere, these will be evaluated using integer division, eg

$ python

>>> print (5/4)

1

(versus)

>>> print (5/4.0)

1.25

and so on.

Thanks John – and very well spotted. Have updated the script so should be fine now.

Hi Hannah,

I’d like to re-use this code within the DART project (www.dartproject.info). I would like to credit you: have you released this under a licence or is it public domain?

Many thanks for this.

Ant

Hi Ant,

Thank you for the message. You’re more than welcome to use the code – although a credit would be lovely. Best of luck with the project!

Hannah

Many thanks

Many thanks for your implementation Hannah! I required the translation for a project and this got me out of a pickle!

I rewrote the implementation as an Objective-C method, which I’ve pasted below in case anybody else might find it useful. Easting and northing are set at the top of the method and lat/lon returned.

Thanks again!

__________________________

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

}

Amazing! Thanks Dave – I was busy trying to convert the script but having trouble – you’ve saved me!

Hello Hannah,

thanks for putting this script together. I’m having trouble running it in python though. I am a newbie at python so don’t have a good understanding of what might be wrong but wonder if anyone can help.

Using python 2.7.3 the script falls over at the iteration over the BNG.csv file with the error ‘File “OSGB36toWGS84.py”, line 115, in for E,N in BNG: ValueError: too many values to unpack’.

I’ve checked out this error on Stackoverflow and think it’s to do with line 150 ‘for E, N in BNG:’ finding one item (E and N) instead of two separate E and N values. But I’m having real problems understanding how to proceed. Some help would be v welcome.

Cheers, Mark

Found the error to be with the .csv file I was reading into python. I was running OSGB36toWGS84.py from R and the .csv R was writing had too many items in each row.

Thanks again

Thank you very much for this code, you and Python are amazing!

NodeJS ImplementationMany thanks for this code, got me out of a big hole. I’ve converted this to Javascript to use under Node. However, this will work just as well in the browser.

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

}

vans authentic 通販 メンズ トレーナー 着こなし http://meirixinp.breitlingcarrier2013.com/

Thank you so much for this script Hannah. It’s saved me a bundle of time whilst trying to learn Python.

One quick question the for loop at the end [105-107] – if I run this multiple times, will it append data to the BNGandLatLon.csv or simply overwrite?

seem to have answered my own question here:

the line “outputFile = open(‘BNGandLatLon.csv’, ‘wb’)” [Line 100]

Opens a file for writing only in binary format. Overwrites the file if the file exists. If the file does not exist, creates a new file for writing.

However if changed to “outputFile = open(‘BNGandLatLon.csv’, ‘ab’)”

The result is:

Opens a file for appending in binary format. The file pointer is at the end of the file if the file exists. That is, the file is in the append mode. If the file does not exist, it creates a new file for writing.

Info from <a href="http://www.tutorialspoint.com/python/python_files_io.htm" title="Here"?

Hi Hannah,

thanks a lot for the code! Just wanted you to know that it’s still being used. In my case, to do work on flood risk for Scotland.

Cheers!

Andrew

This was massively helpful for me. Many thanks for posting this. I converted it to C# which I’ve copied below so others can use it if they need to. I tried to give variables descriptive names. If anyone spots places where I used any misleading ones due to misunderstanding the transformation please let me know.

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

}

}

}

WEP isn’t that secure anymore, but you may not have a real choice with older technology.

You have a couple of choices here as you will not be capable to add Wi-Fi networks that are locked for your i – Phone with no

password.

Since these electromagnetic waves, like RFID, Wi – Fi, and Bluetooth, issue inside

a certain volume of radiation, there exists cause for concern if these SAR levels

are too high. A Wi-Fi connection is established with all the destination of the pictures and all sorts of the data is saved on your hard disk drive, with no extra work.

Wifi extender options If you’ve got a wireless card or built-in WLAN

in your notebook and the other laptops all have Wi – Fi (represents Wireless Fidelity), the solution is both basic and free.

Naturally, each of the tips and advice given above might

seem hard operate in itself but following them does

pay off in the long run. The first attribute from the article is good talking

about writing concise, in addition to knitwear. home setting up apple airport

express as wifi extender Linked – In headlines are

searchable fields using the ‘People Search’ function when someone is looking for particular skills, interests, qualifications, or

credentials. One thing that you want to start out thinking

is always that all salesletters are is salesmanship in publications form.