B4J Question OSGB36 Offsets To Decimal Lat/Long

RichardN · Nov 2, 2023

I see there are several examples of code to convert from decimal Lat/Long to OSGB36 coordinates on the forum.

I have to carry out the reverse process in converting a large number of X+Y offset cordinates based on OSGB36 to decimal Lat + Long. Has anyone been down this path before that might share their experiences?

Example input: X=533716.6, Y=170817.8 Presumably the OS datum used is 49N 02W?

agraham · Nov 2, 2023

It's not clear what you want to do. OSGB36 coordinates are Lat and Long but on a slightly different ellipsoid to WGS84. When you say X and Y coordinates it sounds more like OS National Grid coordinates.

This library, though old, is pure Java and might do what you want. I've used it on on B4J so it should work for you.

GPStoOSGB library

Mapping and location services have a large emphasis in Android. Most phones appear to be equipped with GPS and, as with Basic4ppc, I am sure that people will be writing own GPS applications. Useful only for those in Great Britain this is a port of my Basic4ppc GPStoOSGB library. For those who...

www.b4x.com

RichardN · Nov 2, 2023

Thanks @agraham,

TBH I don't think the client understands the information he has given me and after playing with it a little I am forming the same conclusion about OSNG. Does library you quote convert both ways?

MicroDrie · Nov 2, 2023

I did find a conversion solution that uses a standard library, but I see that the results differ from your example. A problem that was also identified in 2014 with the GPStoOSGB library among various solutions.

Easting: 533716.599990306, Northing: 170817.79999874375
Latitude: 51.42008789194193, Longitude: -0.07669814834550459
-----
Latitude: 51.5074, Longitude: -0.1278
Easting: 529915.825972819, Northing: 180433.99387455918

Latitude: 51.420087891950935, Longitude: -0.07669814820568178
Easting: 533716.599990306, Northing: 170817.79999874375

MicroDrie · Nov 2, 2023

RichardN said:
Thanks @agraham,

TBH I don't think the client understands the information he has given me and after playing with it a little I am forming the same conclusion about OSNG. Does library you quote convert both ways?

My example use a easy-to-use Java library for converting between latitude and longitude coordinates and easting and northing coordinates in both ways.

agraham · Nov 3, 2023

RichardN said:
Does library you quote convert both ways?

It can convert between all of OS Easting and Northings to OS National Grid to OSGB32 to WGS84 in both directions although you generally need to go through OSGB32 as an intermediate step.

Brian Dean · Nov 3, 2023

I did find a conversion solution that uses a standard library, but I see that the results differ from your example.

I have not seen any mention in this thread of corrections for survey errors. After using Transverse Mercator Projection to convert from WGS84 to OSNG it is necessary to correct for survey errors in the grid which can be as much as 100 metres. I do not know if the @agraham library includes this correction; most libraries simply convert between geodetic and cartesian coordinates, as one would expect.

I attach a routine which corrects for such errors in OSNG by bilinear interpolation from a 100km grid of offsets. In 2002 and again in 2015 OSGB published a set of offsets on a one kilometer grid (OSTN02/OSTN15) but I do not have routines for those in B4X.

RichardN · Nov 3, 2023

Thanks @agraham

It can convert between all of OS Easting and Northings to OS National Grid to OSGB32 to WGS84 in both directions although you generally need to go through OSGB32 as an intermediate step.

I believe the following code is what you are driving at, I hope I got this right.....

Polygon showing the London Borough of Lambeth:

Sub gmap_Ready
    
    Dim Points As List
    Points.Initialize
    Dim Lat,Lon As Double
    Dim LL As LatLng
    Dim NatGridRef as String
    
    Do While Borough.NextRow
    
        Lat = Borough.GetDouble("X")                    'Data typically x = 51.420045 y = -0.76722977
        Lon = Borough.GetDouble("Y")
        
        NatGridRef  = OSGB.NEtoNGR(Lat,Lon)
        OSGB.NGRtoOSGB36(NatGridRef)
        
        LL.Initialize(OSGB.NewLatitude,OSGB.NewLongitude)
        Points.Add(LL)
        
    Loop
    
    gmap.AddPolygon(Points,1,fx.Colors.Red,fx.Colors.Transparent,0)
    
    Dim cp As CameraPosition
    cp.Initialize(LL.Latitude,LL.Longitude,15)                'zoom to the last point
    gmap.MoveCamera(cp)
    
End Sub

Let me first say that the polygon data comes direct from government Ordnance Survey sources so I have to assume that it is 100% accurate. Using the metodology above every point seems to be plotted about 110°T/60m from where it is supposed to be but the error appears constant so I should be able to remove it with a simple offset. This will serve my purposes just fine so I will accept the small error. Frankly the mathematics behind this takes me back to Astro-Navigation from first principles. It makes my head spin so i think I will leave it at that!

RichardN · Nov 3, 2023

Using a simple decimal degrees offset the plot of all the city boroughs has worked out quite well. The screen-grab looks slightly jaggy but the screen presentation is quite sharp:

B4X:

        xOffset = -0.0015
        yOffset = 0.00055

Brian Dean · Nov 3, 2023

... every point seems to be plotted about 110°T/60m from where it is supposed to be

Yes - if this fits your purposes it should be fine. Its not as if you are trying to settle some dispute over a the line of a garden fence!

agraham · Nov 3, 2023

I'm surprised about the 60m offset. I've used this library for many years to take GPS coordinates and map them onto OS 25K and 50K digital map images and the errors I have seen are typically a maximum of 10m some of which can be down to the individual phones' GPS receiver. All my phones give slightly different GPS coordinates when sitting next to each other but are all normally within a few metres of each other. But, as long as it works for you that's fine.

By the way the library uses the OS publlished 7 parameter Helmert transformation algorithms to do the transformations to and from OSGB36 and National grid. The WGS84 transformations are also 7 parameter Helmert transformations published a Roger Muggleton. Yes the maths do my head in as well.

The Carabus Web Pages

roumei · Nov 3, 2023

This is my code to convert from EPSG:27700 (OSGB36 / British National Grid) to WGS84 (latitude and longitude) and back. It uses the Helmert Transformation.
You can check the results with BGS's online converter: https://webapps.bgs.ac.uk/data/webservices/convertForm.cfm#bngToLatLng

B4X:

Sub Class_Globals
    Private Root As B4XView
    Private xui As XUI
    
    Private Const Airy1830_A As Double = 6377563.396
    Private Const Airy1830_B As Double = 6356256.909
    Private Const WGS84_a As Double = 6378137
    Private Const WGS84_b As Double = 6356752.3142
    Private Const MathPIdiv180 As Double = cPI / 180
    Private Const Deg2Rad As Double = 2*cPI/360
    Private Const Rad2Deg As Double = 360 / (2 * cPI)
    
    Type eGPSLocation(Latitude As Double, Longitude As Double)
    Type eProjCoordinate(Easting As Double, Northing As Double, Elevation As Double)
    
End Sub

Public Sub Initialize
'    B4XPages.GetManager.LogEvents = True
End Sub

'This event will be called once, before the page becomes visible.
Private Sub B4XPage_Created (Root1 As B4XView)
    Root = Root1
    Root.LoadLayout("MainPage")
    
    
    Log("EPSG:27700 to WGS84")
    Dim Easting As Double = 533717
    Dim Northing As Double = 170818
    Dim eg As eGPSLocation = GetLatLonFromEPSG2770Coordinate(Easting, Northing)
    Log("EPSG:27700: " & Easting & " " & Northing)
    Log("WGS84: " & NumberFormat2(eg.Latitude, 1, 6, 1, False) & " " & NumberFormat2(eg.Longitude, 1, 6, 1, False))
    
    Log(" ")
    
    Log("WGS84 to EPSG:27700")
    Dim Latitude As Double = 51.420611
    Dim Longitude As Double = -0.078301
    Dim ec As eProjCoordinate = GetEPSG2770CoordinateFromLatLon(Latitude, Longitude)
    Log("WGS84: " & Latitude & " " & Longitude)
    Log("EPSG:27700: " & NumberFormat2(ec.Easting, 0, 0, 0, False) & " " & NumberFormat2(ec.Northing, 0, 0, 0, False))
        
    
End Sub

Public Sub GetEPSG2770CoordinateFromLatLon(Latitude As Double, Longitude As Double) As eProjCoordinate
    
    Dim eg As eGPSLocation = HelmertTransformation(Latitude, Longitude, Airy1830_A, Airy1830_B, WGS84_a, WGS84_b, _
                                                   446.448, -125.157, 542.06, _
                                                   0.1502 / 3600 * cPI / 180, 0.247 / 3600 * cPI / 180,  0.8421 / 3600 * cPI / 180, -20.4894 / 1000000 + 1, _
                                                   False)
    Dim e As eProjCoordinate = GPS_GetTransverseMercator_FromLatLon(eg.Latitude, eg.Longitude, 400000, -100000, -2, 0.9996012717, 49, Airy1830_A, Airy1830_B)
    Return e
    
End Sub

Public Sub GetLatLonFromEPSG2770Coordinate(Easting As Double, Northing As Double) As eGPSLocation
    
    Dim e As eGPSLocation = GPS_GetLatLon_FromTransverseMercator(Easting, Northing, 400000, -100000, -2, 0.9996012717, 49, Airy1830_A, Airy1830_B)
    Dim eg As eGPSLocation = HelmertTransformation(e.Latitude, e.Longitude, Airy1830_A, Airy1830_B, WGS84_a, WGS84_b, _
                                                   446.448, -125.157, 542.06, _
                                                   0.1502 / 3600 * cPI / 180, 0.247 / 3600 * cPI / 180,  0.8421 / 3600 * cPI / 180, -20.4894 / 1000000 + 1, _
                                                   True)
    Return eg
End Sub

Private Sub HelmertTransformation(Latitude As Double, _
                                  Longitude As Double, _
                                  a_in As Double, b_in As Double, a_out As Double, b_out As Double, tx As Double, ty As Double, tz As Double, _
                                  rx As Double, ry As Double, rz As Double, s1 As Double, _
                                  bToWGS84 As Boolean) As eGPSLocation
    
    If bToWGS84 = False Then
        Dim temp_a_in As Double = a_in
        Dim temp_b_in As Double = b_in
        Dim temp_a_out As Double = a_out
        Dim temp_b_out As Double = b_out
        a_in = temp_a_out
        b_in = temp_b_out
        a_out = temp_a_in
        b_out = temp_b_in
        rx = (-1)*rx
        ry = (-1)*ry
        rz = (-1)*rz
        tx = (-1)*tx
        ty = (-1)*ty
        tz = (-1)*tz
        s1 = (-(s1 - 1) * 1000000) / 1000000 + 1
    End If

    Dim lat As Double = Latitude * Deg2Rad
    Dim lon As Double = Longitude * Deg2Rad

    Dim a As Double = a_in
    Dim b As Double = b_in
    Dim sinPHI As Double = Sin(lat)
    Dim cosPHI As Double = Cos(lat)
    Dim sinLambda As Double = Sin(lon)
    Dim cosLambda As Double = Cos(lon)
    Dim H As Double = 0 ' Höhe, wird hier nicht verwendet
    Dim eSq As Double = (a * a - b * b) / (a * a)
    Dim nu As Double = a / Sqrt(1 - eSq * sinPHI * sinPHI)
    Dim x1 As Double = (nu + H) * cosPHI * cosLambda
    Dim y1 As Double = (nu + H) * cosPHI * sinLambda
    Dim z1 As Double = ((1 - eSq) * nu + H) * sinPHI

    Dim x2 As Double = tx + x1 * s1 - y1 * rz + z1 * ry
    Dim y2 As Double = ty + x1 * rz + y1 * s1 - z1 * rx
    Dim z2 As Double = tz - x1 * ry + y1 * rx + z1 * s1

    a = a_out
    b = b_out
    Dim precision As Double = 2 / a
    eSq = (a * a - b * b) / (a * a)
    Dim p As Double = Sqrt(x2 * x2 + y2 * y2)
    Dim phi As Double = ATan2(z2, p * (1 - eSq))
    Dim phiP As Double = 2 * cPI
    
    Do While (Abs(phi - phiP) > precision)
        nu = a / Sqrt(1 - eSq * Sin(phi) * Sin(phi))
        phiP = phi
        phi = ATan2(z2 + eSq * nu * Sin(phi), p)
    Loop

    Dim lambda As Double = ATan2(y2, x2)
    H = p / Cos(phi) - nu

    Dim egps As eGPSLocation
    egps.Latitude = phi * Rad2Deg
    egps.Longitude = lambda * Rad2Deg
    
    Return egps
    
End Sub

Private Sub GPS_GetTransverseMercator_FromLatLon(Latitude As Double, _
                                                         Longitude As Double, _
                                                         FalseEasting As Double, _
                                                         FalseNorthing As Double, _
                                                         CentralMeridian As Double, _
                                                         CentralMeridianScaleFactor As Double, _
                                                         LatitudeOfOrigin As Double, _
                                                         Ellipsoid_A As Double, _
                                                         Ellipsoid_B As Double) As eProjCoordinate

    Dim Northing, Easting, eccPrimeSquared, nu, T, C, A, M, M0 As Double
    Dim eccSquared As Double = Ellipsoid_B / Ellipsoid_A
    eccSquared = 1 - eccSquared * eccSquared

    If Longitude < -180 Then Longitude = Longitude + 360
    If Longitude > 180 Then Longitude = Longitude - 360
    Dim LatRad As Double = Latitude * MathPIdiv180
    Dim LongRad As Double = Longitude * MathPIdiv180
    Dim LongOriginRad As Double = CentralMeridian * MathPIdiv180
    Dim LatOriginRad As Double = LatitudeOfOrigin * MathPIdiv180
    eccPrimeSquared = eccSquared / (1 - eccSquared)
    nu = Ellipsoid_A / Sqrt(1 - eccSquared * Sin(LatRad) * Sin(LatRad))
    T = Tan(LatRad) * Tan(LatRad)
    C = eccPrimeSquared * Cos(LatRad) * Cos(LatRad)
    A = Cos(LatRad) * (LongRad - LongOriginRad)
    M = Ellipsoid_A * ((1 - eccSquared / 4 - 3 * eccSquared * eccSquared / 64 - 5 * eccSquared * eccSquared * eccSquared / 256) * LatRad - (3 * eccSquared / 8 + 3 * eccSquared * eccSquared / 32 + 45 * eccSquared * eccSquared * eccSquared / 1024) * Sin(2 * LatRad) + (15 * eccSquared * eccSquared / 256 + 45 * eccSquared * eccSquared * eccSquared / 1024) * Sin(4 * LatRad) - (35 * eccSquared * eccSquared * eccSquared / 3072) * Sin(6 * LatRad))
    M0 = Ellipsoid_A * ((1 - eccSquared / 4 - 3 * eccSquared * eccSquared / 64 - 5 * eccSquared * eccSquared * eccSquared / 256) * LatOriginRad - (3 * eccSquared / 8 + 3 * eccSquared * eccSquared / 32 + 45 * eccSquared * eccSquared * eccSquared / 1024) * Sin(2 * LatOriginRad) + (15 * eccSquared * eccSquared / 256 + 45 * eccSquared * eccSquared * eccSquared / 1024) * Sin(4 * LatOriginRad) - (35 * eccSquared * eccSquared * eccSquared / 3072) * Sin(6 * LatOriginRad))
    Easting = (CentralMeridianScaleFactor * nu * (A + (1 - T + C) * A * A * A / 6 + (5 - 18 * T + T * T + 72 * C - 58 * eccPrimeSquared) * A * A * A * A * A / 120))
    Northing = (CentralMeridianScaleFactor * (M - M0 + nu * Tan(LatRad) * (A * A / 2 + (5 - T + 9 * C + 4 * C * C) * A * A * A * A / 24 + (61 - 58 * T + T * T + 600 * C - 330 * eccPrimeSquared) * A * A * A * A * A * A / 720)))

    Dim Coord As eProjCoordinate
    Coord.Easting = Easting + FalseEasting
    Coord.Northing = Northing + FalseNorthing
    Return Coord

End Sub

Private Sub GPS_GetLatLon_FromTransverseMercator(Easting As Double, _
                                                          Northing As Double, _
                                                          FalseEasting As Double, _
                                                          FalseNorthing As Double, _
                                                          CentralMeridian As Double, _
                                                          CentralMeridianScaleFactor As Double, _
                                                          LatitudeOfOrigin As Double, _
                                                          Ellipsoid_A As Double, _
                                                          Ellipsoid_B As Double) As eGPSLocation

    Dim Lat, Lon, N1, T1, C1, R1, D, M1, mu1, phi1Rad As Double
    Dim LongOriginRad As Double = CentralMeridian * MathPIdiv180
    Dim LatOriginRad As Double = LatitudeOfOrigin * MathPIdiv180
    Dim eccSquared As Double = Ellipsoid_B / Ellipsoid_A
    eccSquared = 1 - eccSquared * eccSquared
    Dim e1 As Double = (1 - Sqrt(1 - eccSquared)) / (1 + Sqrt(1 - eccSquared))
    Dim eccPrimeSquared As Double = (eccSquared) / (1 - eccSquared)
    Dim M0 As Double = Ellipsoid_A * ((1 - eccSquared / 4 - 3 * eccSquared * eccSquared / 64 - 5 * eccSquared * eccSquared * eccSquared / 256) * LatOriginRad - (3 * eccSquared / 8 + 3 * eccSquared * eccSquared / 32 + 45 * eccSquared * eccSquared * eccSquared / 1024) * Sin(2 * LatOriginRad) + (15 * eccSquared * eccSquared / 256 + 45 * eccSquared * eccSquared * eccSquared / 1024) * Sin(4 * LatOriginRad) - (35 * eccSquared * eccSquared * eccSquared / 3072) * Sin(6 * LatOriginRad))
    M1 = M0 + (Northing - FalseNorthing) / CentralMeridianScaleFactor
    mu1 = M1 / (Ellipsoid_A * (1 - eccSquared / 4 - 3 * eccSquared * eccSquared / 64 - 5 * eccSquared * eccSquared * eccSquared / 256))
    phi1Rad = mu1 + (3 * e1 / 2 - 27 * e1 * e1 * e1 / 32) * Sin(2 * mu1) + (21 * e1 * e1 / 16 - 55 * e1 * e1 * e1 * e1 / 32) * Sin(4 * mu1) + (151 * e1 * e1 * e1 / 96) * Sin(6 * mu1)
    N1 = Ellipsoid_A / Sqrt(1 - eccSquared * Sin(phi1Rad) * Sin(phi1Rad))
    T1 = Tan(phi1Rad) * Tan(phi1Rad)
    C1 = eccPrimeSquared * Cos(phi1Rad) * Cos(phi1Rad)
    R1 = Ellipsoid_A * (1 - eccSquared) / Power(1 - eccSquared * Sin(phi1Rad) * Sin(phi1Rad), 1.5)
    D = (Easting - FalseEasting) / (N1 * CentralMeridianScaleFactor)
    Lat = phi1Rad - (N1 * Tan(phi1Rad) / R1) * (D * D / 2 - (5 + 3 * T1 + 10 * C1 - 4 * C1 * C1 - 9 * eccPrimeSquared) * D * D * D * D / 24 + (61 + 90 * T1 + 298 * C1 + 45 * T1 * T1 - 252 * eccPrimeSquared - 3 * C1 * C1) * D * D * D * D * D * D / 720)
    Lon = (D - (1 + 2 * T1 + C1) * D * D * D / 6 + (5 - 2 * C1 + 28 * T1 - 3 * C1 * C1 + 8 * eccPrimeSquared + 24 * T1 * T1) * D * D * D * D * D / 120) / Cos(phi1Rad)

    Dim e As eGPSLocation
    e.Latitude = Lat * (180.0 / cPI)
    e.Longitude = (Lon + LongOriginRad) * (180.0 / cPI)
    Return e

End Sub

Result:

EPSG:27700 to WGS84
EPSG:27700: 533717 170818
WGS84: 51.420611 -0.078301

WGS84 to EPSG:27700
WGS84: 51.420611 -0.078301
EPSG:27700: 533717 170818

RichardN · Nov 4, 2023

@roumei

I just tried plotting the official dataset using your EPSG:27700 to WGS84 code. It works perfectly. Congratulations

B4J Question OSGB36 Offsets To Decimal Lat/Long

RichardN

Well-Known Member

roumei

agraham

Expert

GPStoOSGB library

RichardN

Well-Known Member

MicroDrie

Well-Known Member

MicroDrie

Well-Known Member

agraham

Expert

Brian Dean

Well-Known Member

Attachments

RichardN

Well-Known Member

RichardN

Well-Known Member

Attachments

Brian Dean

Well-Known Member

agraham

Expert

roumei

Active Member

RichardN

Well-Known Member

Similar Threads