OS Grid Ref to Lat/Long 1

Felix18807 · Sep 6, 2011

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

MajP · Sep 7, 2011

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

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

MajP · Sep 7, 2011

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
  
  LL.Lattitude = getRadians(getDecimalDegrees(52, 39, 27.2531))
  LL.Longitude = getRadians(getDecimalDegrees(1, 43, 4.5177))
  '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)
  Debug.Print "Lattitude " & degMinSecFromRadians(LL.Lattitude)
  Debug.Print "Longitue " & degMinSecFromRadians(LL.Longitude)
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
Public Function degMinSecFromRadians(rads As Double) As String
  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)
  degMinSecFromRadians = str
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

Felix18807 · Sep 8, 2011

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 RadPHI0 As Double
Dim RadLAM0 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))
'Convert angle measures to radians
    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 RadPHI0 As Double
Dim RadLAM0 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

'Convert angle measures to radians
    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

Tek-Tips is the largest IT community on the Internet today!

Members share and learn making Tek-Tips Forums the best source of peer-reviewed technical information on the Internet!

OS Grid Ref to Lat/Long 1

Felix18807

Programmer

MajP

Technical User

MajP

Technical User

Felix18807

Programmer

Similar threads

Part and Inventory Search

Sponsor