VB.NET - Convert OS Grid to Lat/Long

Convert OS Grid References to Latitude/Longitude

Notes

The page uses algorhythmns explained in "A Guide to Coordinate Systems in Great Britain" published by Ordnance Survey. Where you copy, publish or distribute the contents of this document to third parties, you must acknowledge Ordnance Survey as the source of the information by including the attribution statement Copyright Ordnance Survey 2018
This conversion only applies to Grid References in England, Wales and Scotland. Northern Ireland uses the Irish Grid system
The "Introduction" page in this section explains how to add teh OSTN15 file to your project


Public Function OsXyToLatLon(EN As EastNorth) As LatLon
 Dim ret As LatLon
 Dim a, b, n, eSq, f0, lon0, lat0, n0, e0, v, p, netaSq, m, no, ea, latBar, vii, viii, ix, x, xi, xii, xiia, lat, lon As Double
 ea = EN.East
 no = EN.North
 Dim adj As OSTN = ArrOSTN(ea, no)
 ea = EN.East - adj.X
 no = EN.North - adj.Y
 ' ea and no within a few mm at this point - OSTN15 flexi sheet accuracy will improve if we adjrect at adjusted point
 adj = ArrOSTN(ea, no)
 ea = EN.East - adj.X
 no = EN.North - adj.Y
 ' ea and no now within about 1mm - very close to limit accuracy on OSTN15
 a = 6378137.0
 b = 6356752.31414036
 f0 = 0.9996012717
 lat0 = DegToRad(49)
 lon0 = DegToRad(-2)
 e0 = 400000
 n0 = -100000
 eSq = (Pow2(a) - Pow2(b)) / Pow2(a)
 n = (a - b) / (a + b)
 latBar = ((no - n0) / (a * f0)) + lat0
 m = (b * f0) * ((1 + n + ((5 / 4) * Pow2(n)) + ((5 / 4) * Pow3(n))) * (latBar - lat0) - (((3 * n) + (3 * Pow2(n)) + ((21 / 8) * Pow3(n))) * Math.Sin(latBar - lat0) * Math.Cos(latBar + lat0)) + (((15 / 8) * Pow2(n)) + ((15 / 8) * Pow3(n))) * Math.Sin(2 * (latBar - lat0)) * Math.Cos(2 * (latBar + lat0)) - (35 / 24) * Pow3(n) * Math.Sin(3 * (latBar - lat0)) * Math.Cos(3 * (latBar + lat0)))
 While (Math.Abs(no - n0 - m) >= 0.000001)
  latBar = ((no - n0 - m) / (a * f0)) + latBar
  m = (b * f0) * ((1 + n + ((5 / 4) * Pow2(n)) + ((5 / 4) * Pow3(n))) * (latBar - lat0) - (((3 * n) + (3 * Pow2(n)) + ((21 / 8) * Pow3(n))) * Math.Sin(latBar - lat0) * Math.Cos(latBar + lat0)) + (((15 / 8) * Pow2(n)) + ((15 / 8) * Pow3(n))) * Math.Sin(2 * (latBar - lat0)) * Math.Cos(2 * (latBar + lat0)) - (35 / 24) * Pow3(n) * Math.Sin(3 * (latBar - lat0)) * Math.Cos(3 * (latBar + lat0)))
 End While
 latBar = ((no - n0 - m) / (a * f0)) + latBar
 v = (a * f0) * Pow_Min0_5((1 - eSq * Pow2(Math.Sin(latBar))))
 p = (a * f0) * (1 - eSq) * Pow_Min1_5((1 - (eSq * Pow2(Math.Sin(latBar)))))
 netaSq = (v / p) - 1
 vii = (Math.Tan(latBar) / (2 * p * v))
 viii = (Math.Tan(latBar) / (24 * p * Pow3(v))) * (5 + (3 * Pow2(Math.Tan(latBar))) + netaSq - (9 * Pow2(Math.Tan(latBar)) * netaSq))
 ix = (Math.Tan(latBar) / (720 * p * (v ^ 5))) * (61 + (90 * (Math.Tan(latBar) ^ 2)) + (45 * (Math.Tan(latBar) ^ 4)))
 x = Sec(latBar) / v
 xi = (Sec(latBar) / (6 * (v ^ 3))) * ((v / p) + (2 * (Math.Tan(latBar) ^ 2)))
 xii = (Sec(latBar) / (120 * (v ^ 5))) * (5 + (28 * (Math.Tan(latBar) ^ 2)) + (24 * (Math.Tan(latBar) ^ 4)))
 xiia = (Sec(latBar) / (5040 * (v ^ 7))) * (61 + (662 * (Math.Tan(latBar) ^ 2)) + (1320 * (Math.Tan(latBar) ^ 4)) + (720 * (Math.Tan(latBar) ^ 6)))
 lat = latBar - (vii * ((ea - e0) ^ 2)) + (viii * ((ea - e0) ^ 4)) - (ix * ((ea - e0) ^ 6))
 lon = lon0 + (x * (ea - e0)) - (xi * ((ea - e0) ^ 3)) + (xii * ((ea - e0) ^ 5)) - (xiia * ((ea - e0) ^ 7))
 lat = RadToDeg(lat)
 lon = RadToDeg(lon)
 ret.Lat = lat
 ret.lon = lon
 Return ret
End Function
'
Private Function DegToRad(Deg As Double) As Double
 Return Deg * Math.PI / 180
End Function
'
Private Function RadToDeg(Rad As Double) As Double
 Return Rad * 180 / Math.PI
End Function
'
Private Function Sec(angle As Double) As Double
 Return Math.Sqrt(1 + (Math.Tan(angle) * Math.Tan(angle)))
End Function
'
Private Function Pow2(n As Double) As Double
 Return n * n
End Function
'
Private Function Pow3(n As Double) As Double
 Return n * n * n
End Function
'
Private Function Pow_Min0_5(n As Double) As Double
 Return 1 / (Math.Sqrt(n))
End Function
'
Private Function Pow_Min1_5(n As Double) As Double
 Return 1 / (n * Math.Sqrt(n))
End Function
'
Private Function OstnXY(east As Double, north As Double) As OSTN
 Dim adj As New OSTN
 Dim e1, n1, x0, y0 As Integer
 If east <= 0 Or east >= 700000 Or north <= 0 Or north >= 1250000 Then
  adj.X = 0
  adj.Y = 0
  Return adj
 End If
 If Not OstnIsLoaded() Then
  Dim pat As String = Application.StartupPath
  If Mid(pat, pat.Length, 1) <> "\" Then pat &= "\"
  LoadOSTN(pat & "ostn15.dat")
 End If
 e1 = Math.Floor(east / 1000) : x0 = e1 * 1000
 n1 = Math.Floor(north / 1000) : y0 = n1 * 1000
 Dim se0, se1, se2, se3, sn0, sn1, sn2, sn3, sz0, sz1, sz2, sz3, dx, dy, t, u, se, sn, sz As Double
 se0 = ArrOSTN(e1, n1).X
 sn0 = ArrOSTN(e1, n1).Y
 se1 = ArrOSTN(e1 + 1, n1).X
 sn1 = ArrOSTN(e1 + 1, n1).Y
 se2 = ArrOSTN(e1 + 1, n1 + 1).X
 sn2 = ArrOSTN(e1 + 1, n1 + 1).Y
 se3 = ArrOSTN(e1, n1 + 1).X
 sn3 = ArrOSTN(e1, n1 + 1).Y
 dx = east Mod 1000
 dy = north Mod 1000
 t = dx / 1000
 u = dy / 1000
 se = ((1 - t) * (1 - u) * se0) + (t * (1 - u) * se1) + (t * u * se2) + ((1 - t) * u * se3)
 sn = ((1 - t) * (1 - u) * sn0) + (t * (1 - u) * sn1) + (t * u * sn2) + ((1 - t) * u * sn3)
 sz = ((1 - t) * (1 - u) * sz0) + (t * (1 - u) * sz1) + (t * u * sz2) + ((1 - t) * u * sz3)
 adj.X = se
 adj.Y = sn
 Return adj
End Function
'
Private Function OstnIsLoaded() As Boolean
 If ArrOSTN.Length < 700 Then Return False
 If ArrOSTN(350, 625).X = 0 Then Return False
  Return True
 End Function
'
Public Function LoadOSTN(OstnFile As String) As String
 ' returns empty string if successful or an error message
 ' expects to find 701x1251 entries
 ' each entry is one line containing OSTN-X and OSTN-y Adjustment
 ' entries sorted ascending on Northing_thenBy_Easting
 If OstnIsLoaded() Then Return String.Empty
 Dim ListOSTNpairs As List(Of String) = IO.File.ReadAllLines(OstnFile).ToList
 If ListOSTNpairs.Count <> 701 * 1251 Then Return "Unexpected number of entries in OSTN file"
 Dim parts() As String
 For yy As Integer = 0 To 1250
  For xx As Integer = 0 To 700
   parts = ListOSTNpairs.Item((yy * 701) + xx).Split(","c)
   If parts.Length <> 2 Then Return "Unpaired OSTN entry encountered"
   If Not Double.TryParse(parts(0), ArrOSTN(xx, yy).X) Then Return ("Non numeric data encountered")
   If Not Double.TryParse(parts(1), ArrOSTN(xx, yy).Y) Then Return ("Non numeric data encountered")
  Next
 Next
 ListOSTNpairs.Clear()
 Return String.Empty
End Function
'
Public Structure LatLon
 Dim Lat As Double
 Dim lon As Double
End Structure
'
Public Structure OSTN
 Dim X As Double
 Dim Y As Double
End Structure
'
Public Structure EastNorth
 Dim East As Double
 Dim North As Double
End Structure

Acknowledgement

This page references material Copyright Ordnance Survey 2018 from their publication A Guide to Coordinate Systems in Great Britain

DigitalDan.co.uk