VB.Net Convert GPS to OS Grid

Convert Latitude/Longitude to OS Easting/Northing

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 the OSTN15 file to your project


Public Function LatLonToOsXy(ll As LatLon) As EastNorth
 Dim a, b, eSq, f0, lon0, lat0, n0, e0, n, v, p, netaSq, m, i, ii, iii, iiia, iv, vv, vi, north, east As Double
 Dim ret As EastNorth
 ll.Lat = DegToRad(ll.Lat)
 ll.lon = DegToRad(ll.lon)
 a = 6378137
 b = 6356752.3141
 f0 = 0.9996012717 ' this is a scale factor
 lat0 = DegToRad(49)
 lon0 = DegToRad(-2)
 e0 = 400000
 n0 = -100000
 eSq = (Pow2(a) - Pow2(b)) / Pow2(a)
 n = (a - b) / (a + b)
 v = (a * f0) * Pow_Min0_5((1 - eSq * Pow2(Math.Sin(ll.Lat))))
 p = (a * f0) * (1 - eSq) * Pow_Min1_5((1 - (eSq * (Pow2(Math.Sin(ll.Lat))))))
 netaSq = (v / p) - 1
 m = (b * f0) * ((1 + n + ((5 / 4) * Pow2(n)) + ((5 / 4) * Pow3(n))) * (ll.Lat - lat0) - (((3 * n) + (3 * Pow2(n)) + ((21 / 8) * Pow3(n))) * Math.Sin(ll.Lat - lat0) * Math.Cos(ll.Lat + lat0)) + (((15 / 8) * Pow2(n)) + ((15 / 8) * Pow3(n))) * Math.Sin(2 * (ll.Lat - lat0)) * Math.Cos(2 * (ll.Lat + lat0)) - (35 / 24) * Pow3(n) * Math.Sin(3 * (ll.Lat - lat0)) * Math.Cos(3 * (ll.Lat + lat0)))
 i = m + n0
 ii = (v / 2) * Math.Sin(ll.Lat) * Math.Cos(ll.Lat)
 iii = (v / 24) * Math.Sin(ll.Lat) * Pow3(Math.Cos(ll.Lat)) * (5 - Pow2(Math.Tan(ll.Lat)) + (9 * netaSq))
 iiia = (v / 720) * (Math.Sin(ll.Lat) * Pow5(Math.Cos(ll.Lat))) * (61 - (58 * Pow2(Math.Tan(ll.Lat))) + Pow4(Math.Tan(ll.Lat)))
 iv = v * Math.Cos(ll.Lat)
 vv = (v / 6) * Pow3(Math.Cos(ll.Lat)) * ((v / p) - Pow2(Math.Tan(ll.Lat)))
 vi = ((v / 120) * Pow5(Math.Cos(ll.Lat))) * (5 - (18 * Pow2(Math.Tan(ll.Lat))) + Pow4(Math.Tan(ll.Lat)) + (14 * netaSq) - (58 * Pow2(Math.Tan(ll.Lat)) * netaSq))
 Dim ll0 As Double
 ll0 = ll.lon - lon0
 north = i
 north += (ii * Pow2(ll0))
 north += (iii * Pow4(ll0))
 north += (iiia * Pow6(ll0))
 east = e0
 east += (iv * ll0)
 east += (vv * Pow3(ll0))
 east += (vi * Pow5(ll0))
 ret.East = east
 ret.North = north
 Dim adj As OSTN = OstnXY(east, north)
 ret.East += adj.X
 ret.North += adj.Y
 Return ret
End Function
'
Private Function DegToRad(Deg As Double) As Double
 Return Deg * Math.PI / 180
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 Pow4(n As Double) As Double
 Return n * n * n * n
End Function
'
Private Function Pow5(n As Double) As Double
 Return n * n * n * n * n
End Function
'
Private Function Pow6(n As Double) As Double
 Return n * n * n * 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