Jump to content

OS grid - latitude longitude conversion


RTFC
 Share

Recommended Posts

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:

; 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 by RTFC
Link to comment
Share on other sites

You're welcome, KLM. :)

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

Link to comment
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now
 Share

×
×
  • Create New...