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