Convert Easting and Northing Coordinates to Latitude and Longitude using VBA
Many Government bodies and agencies in the UK still publish postcode data using Easting and Northing co-ordinates. Although this may be useful for some, many software systems, such as Tableau, don’t natively support this and prefer the more universally used Latitude and Longitude coordinate system.
Converting to and from these coordinate systems is not a perfect science. The process is an approximation as it attempts to convert 2-dimensional mappings onto a 3-dimensional object, along with taking account of some of the obscurities of the OSGB36 grid system. You can read more about the UK’s grid system at the Ordinance Survey.
There are a few online converters out there, but nothing seemed particularly relevant for batch converting stuff in Excel or Access.
I did, however, find a SQL UDF written by Sandy Mottram at Carra Lucia that, itself, had been converted from a javascript file hosted on a caving club’s website.
Basically, other than the time taken to convert and error check, I can’t really take credit for this.
You can either refer to the function as part of a VBA routine or use it in a worksheet cell. Call in the usual way:
=NEtoLL(521034,169028,”Lat”) for latitude
=NEtoLL(521034,169028,”Lng”) for longitude
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 |
Function NEtoLL(East As Variant, North As Variant, LatOrLng As String) As Variant On Error GoTo errh If IsNull(East) Then GoTo errh If IsNull(North) Then GoTo errh 'Author: partc.co.uk 'Provenance: UDF converted to VBA from SQL, originally written by Sandy Motteram, found at: https://carra-lucia-ltd.co.uk/2014/07/28/how-to-convert-easting-and-northing-points-to-latitude-and-longitude-for-coordinates-calculation/ ' SQL-based UDF originally adapted from javascript at http://www.bdcc.co.uk/LatLngToOSGB.js 'Instructions: 'Latitude and Longitude are calculated based on BOTH the easting and northing values from the OSGB36 'This UDF takes both easting and northing values in OSGB36 projection and you must specify if a latitude or longitude co-ordinate should be returned. 'It first converts E/N values to lat and long in OSGB36 projection, then converts those values to lat/lng in WGS84 projection Dim Pi As Variant Dim K0 As Variant Dim OriginLat As Variant Dim OriginLong As Variant Dim OriginX As Variant Dim OriginY As Variant Dim a As Variant Dim b As Variant Dim e2 As Variant Dim ex As Variant Dim n1 As Variant Dim n2 As Variant Dim n3 As Variant Dim OriginNorthings As Variant Dim lat As Variant Dim lon As Variant Dim Northing As Variant Dim Easting As Variant Pi = 3.14159265358979 K0 = 0.9996012717 ' grid scale factor on central meridean OriginLat = 49# OriginLong = -2# OriginX = 400000 ' 400 kM OriginY = -100000 ' 100 kM a = 6377563.396 ' Airy Spheroid b = 6356256.91 'e2 'ex 'n1 'n2 'n3 'OriginNorthings 'compute interim values a = a * K0 b = b * K0 n1 = (a - b) / (a + b) n2 = n1 * n1 n3 = n2 * n1 lat = OriginLat * Pi / 180# ' to radians e2 = (a * a - b * b) / (a * a) ' first eccentricity ex = (a * a - b * b) / (b * b) ' second eccentricity OriginNorthings = b * lat + b * (n1 * (1# + 5# * n1 * (1# + n1) / 4#) * lat _ - 3# * n1 * (1# + n1 * (1# + 7# * n1 / 8#)) * Sin(lat) * Cos(lat) _ + (15# * n1 * (n1 + n2) / 8#) * Sin(2# * lat) * Cos(2# * lat) _ - (35# * n3 / 24#) * Sin(3# * lat) * Cos(3# * lat)) Northing = North - OriginY Easting = East - OriginX Dim nu As Variant Dim phid As Variant Dim phid2 As Variant Dim t2 As Variant Dim t As Variant Dim q2 As Variant Dim c As Variant Dim s As Variant Dim nphid As Variant Dim dnphid As Variant Dim nu2 As Variant Dim nudivrho As Variant Dim invnurho As Variant Dim rho As Variant Dim eta2 As Variant 'Evaluate M term: latitude of the northing on the centre meridian Northing = Northing + OriginNorthings phid = Northing / (b * (1# + n1 + 5# * (n2 + n3) / 4#)) - 1# phid2 = phid + 1# Do While (Abs(phid2 - phid) > 0.000001) phid = phid2 nphid = b * phid + b * (n1 * (1# + 5# * n1 * (1# + n1) / 4#) * phid _ - 3# * n1 * (1# + n1 * (1# + 7# * n1 / 8#)) * Sin(phid) * Cos(phid) _ + (15# * n1 * (n1 + n2) / 8#) * Sin(2# * phid) * Cos(2# * phid) _ - (35# * n3 / 24#) * Sin(3# * phid) * Cos(3# * phid)) dnphid = b * ((1# + n1 + 5# * (n2 + n3) / 4#) - 3# * (n1 + n2 + 7# * n3 / 8#) * Cos(2# * phid) _ + (15# * (n2 + n3) / 4#) * Cos(4 * phid) - (35# * n3 / 8#) * Cos(6# * phid)) phid2 = phid - (nphid - Northing) / dnphid Loop c = Cos(phid) s = Sin(phid) t = Tan(phid) t2 = t * t q2 = Easting * Easting nu2 = (a * a) / (1# - e2 * s * s) nu = Sqr(nu2) nudivrho = a * a * c * c / (b * b) - c * c + 1# eta2 = nudivrho - 1 rho = nu / nudivrho invnurho = ((1# - e2 * s * s) * (1# - e2 * s * s)) / (a * a * (1# - e2)) lat = phid - t * q2 * invnurho / 2# + (q2 * q2 * (t / (24 * rho * nu2 * nu) * (5 + (3 * t2) + eta2 - (9 * t2 * eta2)))) lon = (Easting / (c * nu)) _ - (Easting * q2 * ((nudivrho + 2# * t2) / (6# * nu2)) / (c * nu)) _ + (q2 * q2 * Easting * (5 + (28 * t2) + (24 * t2 * t2)) / (120 * nu2 * nu2 * nu * c)) lat = lat * 180# / Pi lon = lon * 180# / Pi + OriginLong 'Now convert the lat and long from OSGB36 to WGS84 Dim OGlat As Variant Dim OGlon As Variant Dim height As Variant OGlat = lat OGlon = lon height = 24 'London's mean height above sea level is 24 metres. Adjust for other locations. Dim deg2rad As Variant Dim rad2deg As Variant Dim radOGlat As Variant Dim radOGlon As Variant deg2rad = Pi / 180 rad2deg = 180 / Pi 'first off convert to radians radOGlat = OGlat * deg2rad radOGlon = OGlon * deg2rad 'these are the values for WGS84(GRS80) to OSGB36(Airy) Dim a2 As Variant Dim h As Variant Dim xp As Variant Dim yp As Variant Dim zp As Variant Dim xr As Variant Dim yr As Variant Dim zr As Variant Dim sf As Variant Dim e As Variant Dim v As Variant Dim x As Variant Dim y As Variant Dim z As Variant Dim xrot As Variant Dim yrot As Variant Dim zrot As Variant Dim hx As Variant Dim hy As Variant Dim hz As Variant Dim newLon As Variant Dim newLat As Variant Dim p As Variant Dim errvalue As Variant Dim lat0 As Variant a2 = 6378137 ' WGS84_AXIS e2 = 6.69438037928458E-03 ' WGS84_ECCENTRIC h = height ' height above datum (from $GPGGA sentence) a = 6377563.396 ' OSGB_AXIS e = 0.0066705397616 ' OSGB_ECCENTRIC xp = 446.448 yp = -125.157 zp = 542.06 xr = 0.1502 yr = 0.247 zr = 0.8421 s = -20.4894 'convert to cartesian; lat, lon are in radians sf = s * 0.000001 v = a / (Sqr(1 - (e * (Sin(radOGlat) * Sin(radOGlat))))) x = (v + h) * Cos(radOGlat) * Cos(radOGlon) y = (v + h) * Cos(radOGlat) * Sin(radOGlon) z = ((1 - e) * v + h) * Sin(radOGlat) 'transform cartesian xrot = (xr / 3600) * deg2rad yrot = (yr / 3600) * deg2rad zrot = (zr / 3600) * deg2rad hx = x + (x * sf) - (y * zrot) + (z * yrot) + xp hy = (x * zrot) + y + (y * sf) - (z * xrot) + yp hz = (-1 * x * yrot) + (y * xrot) + z + (z * sf) + zp 'Convert back to lat, lon newLon = Atn(hy / hx) p = Sqr((hx * hx) + (hy * hy)) newLat = Atn(hz / (p * (1 - e2))) v = a2 / (Sqr(1 - e2 * (Sin(newLat) * Sin(newLat)))) errvalue = 1# lat0 = 0 Do While (errvalue > 0.001) lat0 = Atn((hz + e2 * v * Sin(newLat)) / p) errvalue = Abs(lat0 - newLat) newLat = lat0 Loop 'convert back to degrees newLat = newLat * rad2deg newLon = newLon * rad2deg Dim ReturnMe As Variant ReturnMe = 0 If LatOrLng = "Lat" Then ReturnMe = newLat ElseIf LatOrLng = "Lng" Then ReturnMe = newLon Else ReturnMe = 0 End If NEtoLL = ReturnMe 'Debug.Print ReturnMe Exit Function errh: NEtoLL = Null End Function |
Leave a Reply
Want to join the discussion?Feel free to contribute!