OS grid - latitude longitude conversion

[RTFC]

Two small spherical trig UDFs using Airy's ellipsoid to convert between the British Ordnance Survey grid (easting and northing, in metres) and latitude-longitude coordinates (in decimal degrees), with alternative ellipsoids for Ireland and Channel Islands (see annotations at the end). All four parameters need to be pre-declared (as the new values are parsed ByRef). Source was partly gleaned from here. I doubt whether anyone outside of Britain will find this of use, but you never know, so here goes:

expandcollapse popup

; Example OS grid to lat/lon
$easting=651409 ; in meters
$northing=313177    ; in meters
$lat=0  ; in decimal degrees
$lon=0  ; in decimal degrees

; returns results in pre-declared vars $lat, $lon
_OSgrid_ToLatLon($easting, $northing, $lat, $lon)

MsgBox(0,"OSgrid_TOLatLon", "OS coordinates Easting: " & $easting & ", Northing: " & $northing & @CR & _
        "Latitude: " & $lat & ", Longitude: " & $lon & @CR & "NGR: " & _OSNE2BNG($easting, $northing))


; Example lat/lon to OS grid
$lat = 52.65757
$lon = 1.71791
$easting=0
$northing=0
_OSgrid_FromLatLon($lat, $lon, $easting, $northing)

MsgBox(0,"OSgrid_FROMLatLon", "OS coordinates Easting: " & $easting & ", Northing: " & $northing & @CR & _
        "Latitude: " & $lat & ", Longitude: " & $lon & @CR & "NGR: " & _OSNE2BNG($easting, $northing))



Func _OSgrid_ToLatLon($E, $N, ByRef $lat, ByRef $lon)
; converts Ordnance Survey grid easting and norting, in meters
; into latitude-longitude coordinates, in decimal degrees
; source not recorded; adapted from a code I wrote many years ago.

    Local $deg2rad=3.141592653589793238462643/180
    Local $rad2deg=180/3.141592653589793238462643

    ; Airy 1830, national grid parameters
    Local   $a = 6377563.396        ; semi-major ellipsoidal axis
    Local $b = 6356256.910      ; semi-minor ellipsoidal axis
    Local $N0= -100000          ; northing of true origin
    Local $E0=  400000          ; easting of true origin
    Local $F0= 0.9996012717     ; scale factor on central meridian
    Local $phi0= 49 * $deg2rad  ; latitude of true origin
    Local $lam0= -2 * $deg2rad  ; longitude of true origin and central meridian
    Local $e2=($a^2-$b^2)/$a^2  ; squared ellipticity

    Local $phi_prime=(($N-$N0)/($a*$F0))+$phi0
    Local $M=_GetM($phi_prime)
    Local $iter=0
    While abs($N-$N0-$M)>=0.01 and $iter<=100   ; max iterations to converge
        $iter+=1
        $phi_prime=(($N-$N0-$M)/($a*$F0))+$phi_prime
        $M=_GetM($phi_prime)
    WEnd

    Local $nu   = $a*$F0*(1-$e2*sin($phi_prime)^2)^-0.5
    Local $rho  = $a*$F0*(1-$e2)*(1-$e2*sin($phi_prime)^2)^-1.5
    Local $eta2 = ($nu/$rho)-1
    Local $VII  = tan($phi_prime)/(2*$rho*$nu)
    Local $VIII = (tan($phi_prime)/(24*$rho*$nu^3))*(5+3*tan($phi_prime)^2+$eta2-$eta2*9*(tan($phi_prime)^2))
    Local $IX   = (tan($phi_prime)/(720*$rho*$nu^5))*(61+(90*tan($phi_prime)^2)+45*tan($phi_prime)^4)
    Local $X    = _sec($phi_prime)/$nu
    Local $XI   = (_sec($phi_prime)/(6*$nu^3))*(($nu/$rho)+2*tan($phi_prime)^2)
    Local $XII  = (_sec($phi_prime)/(120*$nu^5))*(5+(28*tan($phi_prime)^2)+24*tan($phi_prime)^4)
    Local $XIIA = (_sec($phi_prime)/(5040*$nu^7))*(61+(662*tan($phi_prime)^2)+(1320*tan($phi_prime)^4)+(720*tan($phi_prime)^6))
    Local $phi  = $phi_prime-($VII*($E-$E0)^2)+($VIII*($E-$E0)^4)-($IX*($E-$E0)^6)
    Local $lamb = $lam0+(($X*($E-$E0))-($XI*($E-$E0)^3)+($XII*($E-$E0)^5)-($XIIA*($E-$E0)^7))

    $lat=$phi*$rad2deg
    $lon=$lamb*$rad2deg

EndFunc


Func _OSgrid_FromLatLon($lat, $lon, ByRef $easting, ByRef $northing)
; converts latitude-longitude coordinates, in decimal degrees
; into Ordnance Survey grid easting and norting, in meters
; http://www.dorcus.co.uk/carabus/ll_ngr.html

    Local $deg2rad=3.141592653589793238462643/180
    Local $rad2deg=180/3.141592653589793238462643
    Local $phi = $lat * $deg2rad; convert latitude to radians
    Local $lam = $lon * $deg2rad; convert longitude to radians
    Local   $a = 6377563.396        ; semi-major ellipsoidal axis
    Local $b = 6356256.910      ; semi-minor ellipsoidal axis
    Local $N0= -100000          ; northing of true origin
    Local $E0=  400000          ; easting of true origin
    Local $F0= 0.9996012717     ; scale factor on central meridian
    Local $e2=($a^2-$b^2)/$a^2  ; squared ellipticity
    Local $phi0= 49 * $deg2rad  ; latitude of true origin
    Local $lam0= -2 * $deg2rad  ; longitude of true origin and central meridian
    Local $af0 = $a * $F0
    Local $bf0 = $b * $F0

    ; easting
    Local $slat2 = sin($phi) * sin($phi)
    Local $nu = $af0 / (sqrt(1 - ($e2 * ($slat2))))
    Local $rho = ($nu * (1 - $e2)) / (1 - ($e2 * $slat2))
    Local $eta2 = ($nu / $rho) - 1
    Local $p = $lam - $lam0
    Local $IV = $nu * cos($phi)
    Local $clat3 = (cos($phi))^3
    Local $tlat2 = tan($phi) * tan($phi)
    Local $V = ($nu / 6) * $clat3 * (($nu / $rho) - $tlat2)
    Local $clat5 = (cos($phi))^5
    Local $tlat4 = (tan($phi))^4
    Local $VI = ($nu / 120) * $clat5 * ((5 - (18 * $tlat2)) + $tlat4 + (14 * $eta2) - (58 * $tlat2 * $eta2))
    $easting = Round($e0 + ($p * $IV) + (($p^3) * $V) + (($p^5) * $VI),0)

    ; northing
    Local $n = ($af0 - $bf0) / ($af0 + $bf0)
    Local $M = _GetM($phi)
    Local $I = $M + ($n0)
    Local $II = ($nu / 2) * sin($phi) * cos($phi)
    Local $III = (($nu / 24) * sin($phi) * (cos($phi)^3)) * (5 - (tan($phi)^2) + (9 * $eta2))
    Local $IIIA = (($nu / 720) * sin($phi) * $clat5) * (61 - (58 * $tlat2) + $tlat4)
    $northing = Round($I + (($p * $p) * $II) + (($p^4) * $III) + (($p^6) * $IIIA),0)

EndFunc


Func _GetM($phi_prime)

    Local $deg2rad=3.141592653589793238462643/180
    Local $rad2deg=180/3.141592653589793238462643
    Local   $a = 6377563.396        ; semi-major ellipsoidal axis
    Local $b = 6356256.910      ; semi-minor ellipsoidal axis
    Local $F0= 0.9996012717     ; scale factor on central meridian
    Local $phi0= 49 * $deg2rad  ; latitude of true origin

    Local $phi_diff=$phi_prime-$phi0
    Local $phi_sum =$phi_prime+$phi0
    Local $n_small =($a-$b)/($a+$b)
    Local $part1   =(1+$n_small+(5/4)*$n_small^2+(5/4)*$n_small^3)*$phi_diff
    Local $part2   =(3*$n_small+3*$n_small^2+(21/8)*$n_small^3)*sin($phi_diff)*cos($phi_sum)
    Local $part3   =((15/8)*$n_small^2+(15/8)*$n_small^3)*sin(2*$phi_diff)*cos(2*$phi_sum)
    Local $part4   =((35/24)*$n_small^3)*sin(3*$phi_diff)*cos(3*$phi_sum)

    Return $b*$F0*($part1-$part2+$part3-$part4)
EndFunc


Func _sec($theta)

Return 1/cos($theta)
EndFunc


Func _OSNE2BNG($easting, $northing)
; produces OS grid reference letter codes
; http://www.dorcus.co.uk/carabus/ll_ngr.html

    Local $eX = $easting / 500000
    Local $nX = $northing / 500000
    Local $tmp = floor($eX)-5.0 * floor($nX)+17.0

    $nX = 5 * ($nX - floor($nX))
    $eX = 20 - 5.0 * floor($nX) + floor(5.0 * ($eX - floor($eX)))
    if $eX > 7.5 Then $eX+=1
    if $tmp > 7.5 Then $tmp+=1

    Return Chr($tmp + 65) & Chr($eX + 65) & " " & $easting & " " & $northing
EndFunc

#cs
NB for Ireland, use a Modified Airy ellipsoid:

Local $a = 6377340.189
Local $b = 6356034.447
Local $e0 = 200000
Local $n0 = 250000
Local $f0 = 1.000035
Local $e2 = 0.00667054015
Local $lam0 = -0.13962634015954636615389526147909
Local $phi0 = 0.93375114981696632365417456114141

And the last line of _OSNE2BNG becomes:
Return Chr($tmp + 65) & Chr($eX + 65) & " " & $easting & " " & $northing


For the Channel Islands, use the International 1924 ellipsoid:

Local $a = 6378388.000  ; INT24 ED50 semi-major
Local $b = 6356911.946  ; INT24 ED50 semi-minor
Local $e0 = 500000      ; CI easting of false origin
Local $n0 = 0               ; CI northing of false origin
Local $f0 = 0.9996      ; INT24 ED50 scale factor on central meridian
Local $e2 = 0.0067226700223333  ; INT24 ED50 eccentricity squared
Local $lam0 = -0.0523598775598  ; CI false east
Local $phi0 = 0 * $deg2rad          ; CI false north



Func _OSNE2BNG_ChannelIslands($easting, $northing)
; http://www.dorcus.co.uk/carabus/ll_ngr.html

    Local $eX = $easting / 500000
    Local $nX = $northing / 500000
    Local $tmp = floor($eX)-5.0 * floor($nX)+17.0

    $nX = 5 * ($nX - floor($nX))
    $eX = 20 - 5.0 * floor($nX) + floor(5.0 * ($eX - floor($eX)))
    if $eX > 7.5 Then $eX+=1
    if $tmp > 7.5 Then $tmp+=1

   If $northing>5500000 Then
        Return "WA " & StringLeft(String($easting),6) & " " & StringTrimLeft(String($northing),1)
    Else
      Return "WV "  & StringLeft(String($easting),6) & " " & StringTrimLeft(String($northing),1)
    EndIf

EndFunc
#ce

Enjoy.

Edited April 5, 2015 by RTFC

[KLM]

Thanks

[RTFC]

You're welcome, KLM.

No download count on this one, so it's good to hear it's of use to someone.