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…

bng_to_latlon.py

#
# 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

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)