## OSGB36 Easting and Northings - ETRS89 Latitude Longitude Conversion

I had thought that the script I created previously to convert ETRS89 latitude and longitude to OSGB36 eastings and northings would be sufficient for my purposes. I wanted latlng -> east/north primarily because I was building a mapping app which would allow a user to draw on it and thus generate a list of eastings and northings for import in AutoCAD Civil 3D - the idea being the user could draw tram lines etc. within a city, in the browser, using Open Street Map as a background layer; then if they wanted they could import this into C3D as polylines.

The problem came about when I wanted to incorporate some data from Ordnance Survey for postcode regions. These regions are defined as a polyline, but these polylines are specified as OSGB36 eastings and northings, this is a problem for my mapping app as it uses Leaflet which in turn uses ETRS89 lat lngs! So I needed a script to convert back again, from OSGB36 eastings and northings back into ETRS89 latitude / longitude.

This was a little more involved than the previous script - it requires a couple of iterations to get back to ETRS89 latitude and longitude from OSGB36 eastings and northings.

Firstly, we need to reverse the OSTN15 shift. We do this by passing in the OSGB36 eastings and northings to get our first guess of what shift was applied to the ETRS89 eastings and northings. With this new adjusted value we try again to find the shift, improving each time until we're satisfied we're as close as we're going to get - we know we're at this point when the difference between the previous shift and the new shift is less than 0.0001m in both the eastings and northings. With the correct shifts identified, we can reverse these. Simply subtract the shift from the OSGB36 eastings and northings to get the ETRS89 eastings and northings.

OSGB36 Easting | 484228.591 |

OSGB36 Northing | 155542.948 |

Iteration | Easting Shift | Northing Shift |
---|---|---|

1 | 97.6817754548 | -78.8695975497 |

2 | 97.6793378976 | -78.8727917185 |

3 | 97.6793377760 | -78.8727917606 |

ETRS89 Easting | 484130.912 |

ETRS89 Northing | 155621.821 |

We now have the ETRS89 eastings and northings. Next, we need to convert these to latitude and longitude. Again, this is an iterative procedure. Firstly, we need the ellipsoid constants we used when converting ETRS89 to OSGB36; and the projection constants, the projection we're using is the National Grid projection:

Ellipsoid Constants | ||
---|---|---|

Semi-major axis | a | 6,378,137.0000m |

Semi-minor axis | b | 6,356,752.3141m |

Eccentricity squared | e^2 | 0.00669438003551 |

Projection Constants | ||

Scale factor on central meridian | F0 | 0.9996012717 |

True origin, latitude: 49 degrees N | lat0 | 0.855211333477 (radians) |

True origin, longitude; 2 degrees W | lon0 | -0.0349065850399 (radians) |

Map coordinates of true origin | E0 | 400,000m |

Map coordinates of true origin | N0 | -100,000m |

n | 0.0016792203978 |

Now we can begin the second iterative procedure. We begin by calculating lat' and M.

lat' | 0.8953051391 |

M | 255420.6896490931 |

Then we check to see if the absolute value of
*(N - N0 - M) => 0.01mm*. If it is, we obtain a new value for lat', again recomputing M. We keep doing this until
we have the final lat' which satisfies the above constraint.

lat' | 0.8953367047 |

Now the iterative calculations are done, it's just a case of converting the ETRS89 easting and northings to latitude and longitude. Again, this is a fairly lengthy process so I've just copied in the results of each stage below. For the actual calculations see the code on my GitHub page.

p | 6371842.800815761 |

v | 6388631.164616689 |

n2 | 0.0026347737 |

VII | 1.5330891453e-14 |

VIII | 3.0173150202e-28 |

IX | 7.9358993746e-42 |

X | 2.5034231282e-07 |

XI | 4.2101765047e-21 |

XII | 1.3383844198e-34 |

XIIA | 5.1274240608e-48 |

ETRS89 lat | 51.292798 |

ETRS89 lng | -0.793407 |