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