Skip to content

Instantly share code, notes, and snippets.

@lfigueira
Last active March 24, 2021 11:14
Show Gist options
  • Save lfigueira/58ad4a5cea1a5d8a92d7 to your computer and use it in GitHub Desktop.
Save lfigueira/58ad4a5cea1a5d8a92d7 to your computer and use it in GitHub Desktop.
Python script to convert British National grid coordinates (OSGB36 Eastings, Northings) to WGS84 latitude and longitude. The code from this script was copied/adapted from Hannah Fry's blog; the original code and post can be found here: http://www.hannahfry.co.uk/blog/2012/02/01/converting-british-national-grid-to-latitude-and-longitude-ii, toghe…
#
# Converts eastings and northings (British national grid coordinates) to Lat/Long
#
# Original code author: Hannah Fry; see code/comments here:
# http://www.hannahfry.co.uk/blog/2012/02/01/converting-british-national-grid-to-latitude-and-longitude-ii
#
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
# The Airy 180 semi-major and semi-minor axes used for OSGB36 (m)
a, b = 6377563.396, 6356256.909
F0 = 0.9996012717 # scale factor on the central meridian
lat0 = 49 * pi / 180 # Latitude of true origin (radians)
# Longtitude of true origin and central meridian (radians)
lon0 = -2 * pi / 180
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
# The translations along x,y,z axes respectively
tx, ty, tz = 446.448, -125.157, + 542.060
# The rotations along x,y,z respectively, in seconds
rxs, rys, rzs = 0.1502, 0.2470, 0.8421
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
# The GSR80 semi-major and semi-minor axes used for WGS84(m)
a_2, b_2 = 6378137.000, 6356752.3141
# The eccentricity of the GRS80 ellipsoid
e2_2 = 1 - (b_2 * b_2) / (a_2 * a_2)
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
@Benblob688
Copy link

I believe there is an error with this code.
using the E/N:
(-405000, -625000)
I get lat/lon of:
(48.4818508703034, -12.902610178120215)
the correct answer is:
(43.835377, -12.011950)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment