Converting between an OS grid reference and longitude/latitude

Question

The UK Ordnance Survey is the national mapping agency for Great Britain. Its maps refer to locations using a grid of easting and northing coordinates $(E, N)$ (the British National Grid) instead of a graticule of longitude and latitude, but there is an algorithm for converting between the two, detailed below and in more detail in the OS publication A Guide to Coordinate Systems in Great Britain [PDF]. The topic of geodesy is a complex and interesting one and the interested reader is referred to this document for more information.

A short warning: the following algorithm converted between grid references and longitude/latitude coordinates within the same datum only (for the parameters given here, the OSGB36 on the Airy 1830 ellipsoid, a good local fit to the island of Great Britain). It is not a coordinate transformation, and so cannot be used to transform OS map grid references to other datums. In particular, it will not to convert between grid references and longitude/latitude coordinates within the global WGS84 datum used by most GPS devices).

Datum and projection parameters

These parameters apply to the Airy 1830 ellipsoid associated with the OSGB36 datum and the coordinates and conventions of the transverse Mercator projection used on OS maps.

ParameterValueDescription
$a$6377563.396 mSemi-major axis
$b$6356256.909 mSemi-minor axis
$E_0$400000 mMap coordinates of true origin
$N_0$-100000 m
$\phi_0$$49^\circ\,\mathrm{N}$Longitude of true origin
$\lambda_0$$2^\circ\,\mathrm{W}$Latitude of true origin
$F_0$0.9996012717Scale factor on central meridian
Some further quantities and functions

The square of the ellipsoid's eccentricity constant, $e^2 = \frac{a^2-b^2}{a^2}$.

Define $n = \frac{a-b}{a+b}$.

The following function is used in scaling the latitude coordinate:

\begin{align*} f_M(\phi) = bF_0 & \left[ \left(1 + n + \frac{5}{4}n^2 + \frac{5}{4}n^3\right)(\phi-\phi_0) \right.\\ & - \left( 3n + 3n^2 + \frac{21}{8}n^3\right)\sin(\phi-\phi_0)\cos(\phi+\phi_0)\\ & + \left( \frac{15}{8}n^2 + \frac{15}{8}n^3\right)\sin[2(\phi-\phi_0)]\cos[2(\phi+\phi_0)]\\ & \left. - \frac{35}{24}n^3\sin[3(\phi-\phi_0)]\cos[3(\phi+\phi_0)] \right] \end{align*}

Note: In implementing the algorithm below, ensure that you work in double precision (float numerical types) and express all angles in radians. Remember that southerly latitudes and westerly longitudes are expressed as negative angles (e.g. $2^\circ\,\mathrm{W}$ is $-2^\circ$).

Converting eastings and northings to latitude and longitude

Given the grid coordinates $(E, N)$ this algorithm returns the latitude and longitude $(\phi, \lambda)$.

First compute $\phi'$ by the following iterative procedure

  1. Set $M=0$, $\phi'=0$.

  2. While $|N-N_0-M| \ge 10^{-5}\,\mathrm{m}$ compute the following

    $$ \phi' \rightarrow \phi' + \frac{N-N_0-M}{aF_0}$$

    $$M = f_M(\phi')$$

  3. Using this converged value of $\phi'$, calculate the further parameters

    \begin{align*} \nu &= \frac{aF_0}{\sqrt{1-e^2\sin^2\phi'}}\\ \rho &= aF_0\frac{1-e^2}{(1-e^2\sin^2\phi')^{3/2}}\\ \eta^2 &= \frac{\nu}{\rho}-1 \end{align*}

  4. Now compute the following coefficients

    \begin{align*} c_1 &= \frac{\tan\phi'}{2\rho\nu}\\ c_2 &= \frac{\tan\phi'}{24\rho\nu^3}(5 + 3\tan^2\phi' + \eta^2 - 9\eta^2\tan^2\phi')\\ c_3 &= \frac{\tan\phi'}{720\rho\nu^5}(61 + 90\tan^2\phi' + 45\tan^4\phi')\\ d_1 &= \frac{\sec\phi'}{\nu}\\ d_2 &= \frac{\sec\phi'}{6\nu^3}\left(\frac{\nu}{\rho} + 2\tan^2\phi'\right)\\ d_3 &= \frac{\sec\phi'}{120\nu^5}\left(5 + 28\tan^2\phi' + 24\tan^4\phi'\right)\\ d_4 &= \frac{\sec\phi'}{5040\nu^7}\left(61 + 662\tan^2\phi' + 1320\tan^4\phi' + 720\tan^6\phi'\right) \end{align*}

  5. Finally, calculate the latitude, $\phi$, and longitude, $\lambda$, as:

\begin{align*} \phi &= \phi' - c_1(E-E_0)^2 + c_2(E-E_0)^4 - c_3(E-E_0)^6\\ \lambda &= \lambda_0 + d_1(E-E_0) - d_2(E-E_0)^3 + d_3(E-E_0)^5 - d_4(E-E_0)^7 \end{align*}

Converting latitude and longitude to eastings and northings

Given the latitude and longitude $(\phi, \lambda)$, this algorithm returns the grid coordinates $(E, N)$.

  1. Calculate the parameters

    \begin{align*} \nu &= \frac{aF_0}{\sqrt{1-e^2\sin^2\phi}}\\ \rho &= aF_0\frac{1-e^2}{(1-e^2\sin^2\phi)^{3/2}}\\ \eta^2 &= \frac{\nu}{\rho}-1\\ M &= f_M(\phi) \end{align*}

  2. Compute the parameters

    \begin{align*} a_1 &= M + N_0\\ a_2 &= \frac{\nu}{2}\sin\phi\cos\phi\\ a_3 &= \frac{\nu}{24}\sin\phi\cos^3\phi(5 - \tan^2\phi + 9\eta^2)\\ a_4 &= \frac{\nu}{720}\sin\phi\cos^5\phi(61 - 58\tan^2\phi + \tan^4\phi)\\ b_1 &= \nu\cos\phi\\ b_2 &= \frac{\nu}{6}\cos^3\phi\left(\frac{\nu}{\rho} - \tan^2\phi\right)\\ b_3 &= \frac{\nu}{120}\cos^5\phi\left(5 - 18\tan^2\phi + \tan^4\phi + 14\eta^2 - 58\eta^2\tan^2\phi\right) \end{align*}

  3. The northing, $N$, and easting, $E$, coordinates are then given by

    \begin{align*} N &= a_1 + a_2(\lambda-\lambda_0)^2 + a_3(\lambda-\lambda_0)^4 + a_4(\lambda-\lambda_0)^6\\ E &= E_0 + b_1(\lambda-\lambda_0) + b_2(\lambda-\lambda_0)^3 + b_3(\lambda-\lambda_0)^5 \end{align*}

Exercise

Write a pair of Python functions, ll_to_os(phi, lam) and os_to_ll(E, N) to convert between longitude/latitude and OS grid reference. An example test is given below.

E, N = 544735, 258334    # King's College, Cambridge
print('(E, N) =', (E, N))
phi, lam = os_to_ll(E, N)
print('(φ, λ) = ({:.8f}°, {:.8f}°)'.format(phi, lam))
print('       =', deg_to_dms(phi, 'latitude'), deg_to_dms(lam, 'longitude'))
E, N = ll_to_os(phi, lam)
print('(E, N) =', (E, N))

Output:

(E, N) = (544735, 258334)
(φ, λ) = (52.20380073°, 0.11824087°)
       = 52° 12′ 13.6826″ N 0° 7′ 5.6671″ E
(E, N) = (544734.99998566438, 258333.99999784387)

Solution