×
INTELLIGENT WORK FORUMS
FOR COMPUTER PROFESSIONALS

Are you a
Computer / IT professional?
Join Tek-Tips Forums!
• Talk With Other Members
• Be Notified Of Responses
• Keyword Search
Favorite Forums
• Automated Signatures
• Best Of All, It's Free!

*Tek-Tips's functionality depends on members receiving e-mail. By joining you are opting in to receive e-mail.

#### Posting Guidelines

Promoting, selling, recruiting, coursework and thesis posting is forbidden.

# OS Grid Ref to Lat/Long

## OS Grid Ref to Lat/Long

(OP)
I'm trying to convert from Ordnance Survey Grid Ref to Longtitude/Latitude using VBA.

The forumla for this is freely available (http://www.ordnancesurvey.co.uk/oswebsite/gps/docs/convertingcoordinatesEN.pdf)

I have posted my code below. It seems to be running ok.. Except it is miles out to the east/north! Has anyone ever attempted anything like this before or is there anyone who can spot the error. Many Thanks

#### CODE

Option Compare Database
Option Explicit

Const a As Double = 6377563.396
Const b As Double = 6356256.91
Const e2 As Double = (a - b) / a
Const N0 As Double = -100000
Const E0 As Double = 400000
Const f0 As Double = 0.999601272
Const phi0 As Double = 0.855211333
Const lambda0 As Double = -0.034906585
Const n As Double = (a - b) / (a + b)

Public Function ConvertOSToLatLon(E__1 As Double, n__2 As Double, lonorlat As Variant) As Variant
Dim phi As Double
Dim M As Double
Dim v As Double
Dim p As Double
Dim vii As Double
Dim viii As Double
Dim ix As Double
Dim x As Double
Dim xi As Double
Dim xii As Double
Dim xiia As Double
Dim n2 As Double
Dim e__3 As Double
Dim lon As Double
Dim lat As Double

phi = (n__2 - N0) / (a * f0) + phi0
M = b * f0 * ((1 + n + (5 \ 4) * n * n + (5 \ 4) * n * n * n) * (phi - phi0) - (3 * n + 3 * n * n + (21 \ 8) * n * n * n) * Sin(phi - phi0) * Cos(phi + phi0) + ((15 \ 8) * n * n + (15 \ 8) * n * n * n) * Sin(2 * (phi - phi0)) * Cos(2 * (phi + phi0)) - (35 \ 24) * n * n * n * Sin(3 * (phi - phi0)) * Cos(3 * (phi + phi0)))

Do While n__2 - N0 - M >= 0.01
phi = (n__2 - N0 - M) / (a * f0) + phi
M = b * f0 * ((1 + n + (5 \ 4) * n * n + (5 \ 4) * n * n * n) * (phi - phi0) - (3 * n + 3 * n * n + (21 \ 8) * n * n * n) * Sin(phi - phi0) * Cos(phi + phi0) + ((15 \ 8) * n * n + (15 \ 8) * n * n * n) * Sin(2 * (phi - phi0)) * Cos(2 * (phi + phi0)) - (35 \ 24) * n * n * n * Sin(3 * (phi - phi0)) * Cos(3 * (phi + phi0)))
Loop

v = a * f0 * ((1 - e2 * Sin(phi) * Sin(phi)) ^ -0.5)

p = a * f0 * ((1 - e2 * Sin(phi) * Sin(phi) ^ -0.5)) * (1 - e2)
n2 = v / p - 1

vii = Tan(phi) / (2 * p * v)
viii = Tan(phi) / (24 * p * v * v * v) * (5 + 3 * Tan(phi) * Tan(phi) + n2 - 9 * Tan(phi) * Tan(phi) * n2)
ix = Tan(phi) / (720 * p * ((v) ^ 5)) * (61 + 90 * Tan(phi) * Tan(phi) + 45 * ((Tan(phi)) ^ 4))
x = (1 / Cos(phi)) / v
xi = (1 / Cos(phi)) / (6 * ((v) ^ 3)) * (v / p + 2 * ((Tan(phi) ^ 2)))
xii = (1 / Cos(phi)) / (120 * ((v) ^ 5)) * (5 + 28 * (Tan(phi) ^ 2) + 24 * (Tan(phi) ^ 4))
xiia = (1 / Cos(phi)) / (5040 * ((v) ^ 7)) * (61 + 662 * (Tan(phi) ^ 2) + 1320 * (Tan(phi) ^ 4)) + 720 * (Tan(phi) ^ 6)

e__3 = (E__1 - E0)

lon = (phi - vii * e__3 * e__3 + viii * e__3 * e__3 * e__3 * e__3 - ix * e__3 * e__3 * e__3 * e__3 * e__3 * e__3) * 180 / (4 * Atn(1))
lat = (lambda0 + x * e__3 - xi * e__3 * e__3 * e__3 + xii * e__3 * e__3 * e__3 * e__3 * e__3 - xiia * e__3 * e__3 * e__3 * e__3 * e__3 * e__3 * e__3) * 180 / (4 * Atn(1))

Select Case lonorlat
Case Is = "lon"
ConvertOSToLatLon = lon
Case Is = "lat"
ConvertOSToLatLon = lat
Case Else
ConvertOSToLatLon = "Choose Lon or Lat"
End Select

End Function

### RE: OS Grid Ref to Lat/Long

I believe at a minimum your formula for e2 should be.

Const e2 As Double = 1 - (b ^ 2) / (a ^ 2)

### RE: OS Grid Ref to Lat/Long

I believe this works, but it took about a thousand debug statements to verify.

#### CODE

Option Explicit

Public Type EastingNorthing
Easting As Double
Northing As Double
End Type

Public Type LatLong
Lattitude As Double
Longitude As Double
End Type

Const pi As Double = 3.14159265
Const a As Double = 6377563.396
Const b As Double = 6356256.91
Const e2 As Double = 1 - (b ^ 2) / (a ^ 2)
Const N0 As Double = -100000
Const E0 As Double = 400000
Const F0 As Double = 0.999601272
Const phi0 As Double = 0.855211333
Const lambda0 As Double = -0.034906585

Public Sub testGetEastingNorthing()
'latt = 52 39 27.2351 N
'long = 1 43 4.5177 E
Dim LL As LatLong
Dim EN As EastingNorthing

'Test lat long to Easting Northing
EN = getEastingNorthing(LL)
Debug.Print "Easting " & EN.Easting
Debug.Print "Northing " & EN.Northing

'Test Easting Northing to Lat Long
EN.Easting = 651409.903
EN.Northing = 313177.27
LL = getLatLong(EN)
End Sub

Public Function getEastingNorthing(LL As LatLong) As EastingNorthing
Dim phi As Double
Dim lambda As Double
Dim v As Double
Dim sinPhi As Double
Dim phiMinusPhi0 As Double
Dim rho As Double
Dim eta2 As Double
Dim M As Double
Dim lambdaMinusLambda0 As Double
Dim cosPhi As Double
Dim tanPhi As Double
Dim I As Double
Dim II As Double
Dim III As Double
Dim IIIA As Double
Dim IV As Double
Dim V_ As Double
Dim VI As Double
Dim East As Double
Dim North As Double

phi = LL.Lattitude
lambda = LL.Longitude
sinPhi = Sin(phi)
cosPhi = Cos(phi)
tanPhi = Tan(phi)
lambdaMinusLambda0 = lambda - lambda0

v = getV(phi)
rho = getRho(phi)
eta2 = getEta2(phi)
M = getM(phi)
I = M + N0
Debug.Print "I " & I
II = v / 2 * sinPhi * cosPhi
Debug.Print "II " & II
III = (v / 24) * sinPhi * cosPhi ^ 3 * (5 - tanPhi ^ 2 + 9 * eta2)
Debug.Print "III " & III
IIIA = (v / 720) * sinPhi * cosPhi ^ 5 * (61 - 58 * tanPhi ^ 2 + tanPhi ^ 4)
Debug.Print "IIIA" & IIIA
IV = v * cosPhi
Debug.Print "IV " & IV
V_ = v / 6 * cosPhi ^ 3 * (v / rho - tanPhi ^ 2)
Debug.Print "V " & V_
VI = v / 120 * cosPhi ^ 5 * (5 - 18 * tanPhi ^ 2 + tanPhi ^ 4 + 14 * eta2 - 58 * tanPhi ^ 2 * eta2)
Debug.Print "VI " & VI
North = I + II * lambdaMinusLambda0 ^ 2 + III * lambdaMinusLambda0 ^ 4 + IIIA * lambdaMinusLambda0 ^ 6
East = E0 + IV * (lambdaMinusLambda0) + V_ * lambdaMinusLambda0 ^ 3 + VI * lambdaMinusLambda0 ^ 5
getEastingNorthing.Easting = East
getEastingNorthing.Northing = North
End Function

Public Function getLatLong(EN As EastingNorthing) As LatLong
Dim phiPrime As Double
Dim phiPrimeNew As Double
Dim phi As Double
Dim tanPhi As Double
Dim secPhi As Double
Dim rho As Double
Dim v As Double
Dim eta2 As Double
Dim VII As Double
Dim VIII As Double
Dim IX As Double
Dim X As Double
Dim XI As Double
Dim XII As Double
Dim XIIA As Double
Dim lat As Double
Dim lon As Double
Dim EminusE0
phiPrime = getPhiPrime(EN.Northing)
phiPrimeNew = getPhiPrimeNew(EN.Northing, phiPrime)
phi = phiPrimeNew
tanPhi = Tan(phi)
secPhi = 1 / Cos(phi)
v = getV(phi)
rho = getRho(phi)
eta2 = getEta2(phi)
VII = tanPhi / (2 * rho * v)
Debug.Print "VII " & VII
VIII = (tanPhi / (24 * rho * v ^ 3)) * (5 + 3 * tanPhi ^ 2 + eta2 - 9 * (tanPhi ^ 2) * eta2)
Debug.Print "VIII " & VIII
IX = (tanPhi / (720 * rho * v ^ 5)) * (61 + 90 * tanPhi ^ 2 + 45 * tanPhi ^ 4)
Debug.Print "IX " & IX
X = secPhi / v
Debug.Print "X " & X
XI = (secPhi / (6 * v ^ 3)) * (v / rho + 2 * tanPhi ^ 2)
Debug.Print "XI " & XI
XII = (secPhi / (120 * v ^ 5)) * (5 + 28 * tanPhi ^ 2 + 24 * tanPhi ^ 4)
Debug.Print "XII " & XII
XIIA = (secPhi / (5040 * v ^ 7)) * (61 + 662 * tanPhi ^ 2 + 1320 * tanPhi ^ 4 + 720 * tanPhi ^ 6)
Debug.Print "XIIA " & XIIA
EminusE0 = EN.Easting - E0
getLatLong.Lattitude = phi - VII * EminusE0 ^ 2 + VIII * EminusE0 ^ 4 - IX * EminusE0 ^ 6
Debug.Print "Lat " & getLatLong.Lattitude
getLatLong.Longitude = lambda0 + X * EminusE0 - XI * EminusE0 ^ 3 + XII * EminusE0 ^ 5 - XIIA * EminusE0 ^ 7
Debug.Print "long " & getLatLong.Longitude
End Function
Public Function getRho(phi As Double) As Double
getRho = a * F0 * (1 - e2) * (1 - e2 * Sin(phi) ^ 2) ^ (-1.5)
Debug.Print "rho " & getRho
End Function
Public Function getV(phi As Double) As Double
getV = a * F0 * (1 - e2 * Sin(phi) ^ 2) ^ (-0.5)
Debug.Print "V " & getV
End Function

Public Function getM(phi As Double) As Double
Dim M As Double
Dim N As Double
Dim n2 As Double
Dim n3 As Double
Dim phiMinusPhi0 As Double
Dim phiPlusPhi0 As Double
phiMinusPhi0 = phi - phi0
phiPlusPhi0 = phi + phi0
N = (a - b) / (a + b)
Debug.Print "n " & N
n2 = N ^ 2
n3 = N ^ 3
M = (1 + N + (5 / 4) * n2 + (5 / 4) * n3) * phiMinusPhi0
'Debug.Print M
M = M - (3 * N + 3 * n2 + (21 / 8) * n3) * Sin(phiMinusPhi0) * Cos(phiPlusPhi0)
'Debug.Print M
M = M + ((15 / 8) * n2 + (15 / 8) * n3) * Sin(2 * phiMinusPhi0) * Cos(2 * phiPlusPhi0)
'Debug.Print M
M = M - (35 / 24) * n3 * Sin(3 * phiMinusPhi0) * Cos(3 * phiPlusPhi0)
'Debug.Print M
M = b * F0 * M
Debug.Print "M " & M
getM = M
End Function
Public Function getEta2(phi As Double) As Double
getEta2 = getV(phi) / getRho(phi) - 1
Debug.Print "eta2 " & getEta2
End Function

Public Function getPhiPrime(N As Double) As Double
getPhiPrime = ((N - N0) / (a * F0)) + phi0
Debug.Print "Phi Prime " & getPhiPrime
End Function

Public Function getPhiPrimeNew(N As Double, phiPrimeOld As Double) As Double
Dim tempPhi As Double
Dim M As Double
M = getM(phiPrimeOld)
tempPhi = ((N - N0 - M) / (a * F0)) + phiPrimeOld
If (N - N0 - M) > 0.01 Then
getPhiPrimeNew = getPhiPrimeNew(N, tempPhi)
Else
getPhiPrimeNew = phiPrimeOld
Debug.Print "Phi Prime New " & getPhiPrimeNew
Exit Function
End If
End Function

Public Function getRadians(deg As Double) As Double
getRadians = deg * (pi / 180)
End Function
Public Function getDecimalDegrees(deg As Long, min As Long, sec As Double)
getDecimalDegrees = deg + min / 60 + sec / 3600
End Function
Dim str As String
Dim decdegs As Double
Dim fullDegs As Integer
Dim decMinutes As Double
Dim fullMinutes As Integer
Dim decSeconds As Double
decdegs = (rads * 180) / pi
fullDegs = Fix(decdegs)
str = fullDegs & "deg "
decMinutes = (decdegs - fullDegs) * 60
fullMinutes = Fix(decMinutes)
str = str & fullMinutes & "' "
decSeconds = (decMinutes - fullMinutes) * 60
str = str & Round(decSeconds, 4)
End Function

Here are the debug results

#### CODE

'Get EN
V 6388502.33516818
rho 6372756.44173155
eta2 2.4708136236813E-03
n 1.67322025032505E-03
M 406688.292429058
I 306688.292429058
II 1540407.91057736
III 156068.755159861
IIIA-20671.1229062775
IV 3875120.58133976
V -170000.78105085
VI -101344.70456831
Easting 651409.903072026
Northing 313177.266802808

'Get Lat Long
Phi Prime 0.920023245543362
n 1.67322025032505E-03
Phi Prime New 0.920066209037565
V 6388523.34338502
rho 6372819.31125047
eta2 2.46422052274897E-03
VII 1.61305624639685E-14
VIII 3.33955473318196E-28
IX 9.41985612153949E-42
X 2.58400624846302E-07
XI 4.69859698192467E-21
XII 1.6124316535575E-34
XIIA 6.65773158313229E-48
Lat 0.919047977367157
long 2.99833878575216E-02
Lattitude 52deg 39' 27.2532
Longitue 1deg 43' 4.5177

### RE: OS Grid Ref to Lat/Long

(OP)
That is far above and beyond what I was expecting MajP and you've certainly taught me a thing or three!

I was just popping on to post a code that works but isn't fantastically accurate. Thank you loads for your help.

#### CODE

Option Explicit
Const a As Double = 6377563.396
Const b As Double = 6356256.91
Const E0 As Double = 400000
Const N0 As Double = -100000
Const F0 As Double = 0.999601272
Const phi0 As Double = 49
Const LAM0 As Double = -2
Const N As Double = (a - b) / (a + b)
Const lambda0 As Double = -0.034906585
Const e2 As Double = (a - b) / a

Function InitialLat(North, N0, afo, phi0, N, bfo)
Dim PHI1 As Double
Dim PHI2 As Double
Dim M As Double

'First PHI value (PHI1)
PHI1 = ((North - N0) / afo) + phi0

'Calculate M
M = Marc(bfo, N, phi0, PHI1)

'Calculate new PHI value (PHI2)
PHI2 = ((North - N0 - M) / afo) + PHI1

'Iterate to get final value for InitialLat
Do While Abs(North - N0 - M) > 0.00001
PHI2 = ((North - N0 - M) / afo) + PHI1
M = Marc(bfo, N, phi0, PHI2)
PHI1 = PHI2
Loop

InitialLat = PHI2

End Function

Function Marc(bf0, N, phi0, phi)

Marc = bf0 * (((1 + N + ((5 / 4) * (N ^ 2)) + ((5 / 4) * (N ^ 3))) * (phi - phi0)) _
- (((3 * N) + (3 * (N ^ 2)) + ((21 / 8) * (N ^ 3))) * (Sin(phi - phi0)) * (Cos(phi + phi0))) _
+ ((((15 / 8) * (N ^ 2)) + ((15 / 8) * (N ^ 3))) * (Sin(2 * (phi - phi0))) * (Cos(2 * (phi + phi0)))) _
- (((35 / 24) * (N ^ 3)) * (Sin(3 * (phi - phi0))) * (Cos(3 * (phi + phi0)))))

End Function

Function E_N_to_Lat(East, North)
Dim pi As Double
Dim af0 As Double
Dim bf0 As Double
Dim e2 As Double
Dim N As Double
Dim Et As Double
Dim PHId As Double
Dim nu As Double
Dim rho As Double
Dim eta2 As Double
Dim VII As Double
Dim VIII As Double
Dim IX As Double

pi = (4 * Atn(1))
pi = 3.14159265358979
RadPHI0 = phi0 * (pi / 180)
RadLAM0 = LAM0 * (pi / 180)

'Compute af0, bf0, e squared (e2), n and Et
af0 = a * F0
bf0 = b * F0
e2 = ((af0 ^ 2) - (bf0 ^ 2)) / (af0 ^ 2)
N = (af0 - bf0) / (af0 + bf0)
Et = East - E0

'Compute initial value for latitude (PHI) in radians
PHId = InitialLat(North, N0, af0, RadPHI0, N, bf0)

'Compute nu, rho and eta2 using value for PHId
nu = af0 / (Sqr(1 - (e2 * ((Sin(PHId)) ^ 2))))
rho = (nu * (1 - e2)) / (1 - (e2 * (Sin(PHId)) ^ 2))
eta2 = (nu / rho) - 1

'Compute Latitude
VII = (Tan(PHId)) / (2 * rho * nu)
VIII = ((Tan(PHId)) / (24 * rho * (nu ^ 3))) * (5 + (3 * ((Tan(PHId)) ^ 2)) + eta2 - (9 * eta2 * ((Tan(PHId)) ^ 2)))
IX = ((Tan(PHId)) / (720 * rho * (nu ^ 5))) * (61 + (90 * ((Tan(PHId)) ^ 2)) + (45 * ((Tan(PHId)) ^ 4)))

E_N_to_Lat = (180 / pi) * (PHId - ((Et ^ 2) * VII) + ((Et ^ 4) * VIII) - ((Et ^ 6) * IX))

End Function

Function E_N_to_Long(East, North) ', a, b, e0, n0, f0, PHI0, LAM0)
Dim pi As Double
Dim af0 As Double
Dim bf0 As Double
Dim e2 As Double
Dim N As Double
Dim Et As Double
Dim PHId As Double
Dim nu As Double
Dim rho As Double
Dim eta2 As Double
Dim X As Double
Dim XI As Double
Dim XII As Double
Dim XIIA As Double

pi = 3.14159265358979
RadPHI0 = phi0 * (pi / 180)
RadLAM0 = LAM0 * (pi / 180)

'Compute af0, bf0, e squared (e2), n and Et
af0 = a * F0
bf0 = b * F0
e2 = ((af0 ^ 2) - (bf0 ^ 2)) / (af0 ^ 2)
N = (af0 - bf0) / (af0 + bf0)
Et = East - E0

'Compute initial value for latitude (PHI) in radians
PHId = InitialLat(North, N0, af0, RadPHI0, N, bf0)

'Compute nu, rho and eta2 using value for PHId
nu = af0 / (Sqr(1 - (e2 * ((Sin(PHId)) ^ 2))))
rho = (nu * (1 - e2)) / (1 - (e2 * (Sin(PHId)) ^ 2))
eta2 = (nu / rho) - 1

'Compute Longitude
X = ((Cos(PHId)) ^ -1) / nu
XI = (((Cos(PHId)) ^ -1) / (6 * (nu ^ 3))) * ((nu / rho) + (2 * ((Tan(PHId)) ^ 2)))
XII = (((Cos(PHId)) ^ -1) / (120 * (nu ^ 5))) * (5 + (28 * ((Tan(PHId)) ^ 2)) + (24 * ((Tan(PHId)) ^ 4)))
XIIA = (((Cos(PHId)) ^ -1) / (5040 * (nu ^ 7))) * (61 + (662 * ((Tan(PHId)) ^ 2)) + (1320 * ((Tan(PHId)) ^ 4)) + (720 * ((Tan(PHId)) ^ 6)))

E_N_to_Long = (180 / pi) * (RadLAM0 + (Et * X) - ((Et ^ 3) * XI) + ((Et ^ 5) * XII) - ((Et ^ 7) * XIIA))

End Function

#### Red Flag This Post

Please let us know here why this post is inappropriate. Reasons such as off-topic, duplicates, flames, illegal, vulgar, or students posting their homework.

#### Red Flag Submitted

Thank you for helping keep Tek-Tips Forums free from inappropriate posts.
The Tek-Tips staff will check this out and take appropriate action.

Close Box

# Join Tek-Tips® Today!

Join your peers on the Internet's largest technical computer professional community.
It's easy to join and it's free.

Here's Why Members Love Tek-Tips Forums:

• Talk To Other Members
• Notification Of Responses To Questions
• Favorite Forums One Click Access
• Keyword Search Of All Posts, And More...

Register now while it's still free!