An Englishman In Google Maps - Plotting Easting/Northing on Google Maps

Plotting Easting/Northing on Google Maps or Virtual Earth

I work for a public sector organisation in the UK and wanted to be able to display various bits of information using Google Maps. (Nothing personal Mr Gates but the first map-related RSS article I read was about how to add Google Maps to your site in 10 minutes - not Virtual Earth.)
Great I thought: we've got the OS Easting/Northing for every local property so putting that on the maps will be easy.

As if.

Anyone who's tried this before will know that it's actually quite tedious and is marginally preferable to drilling your own teeth with a Black & Decker. Personally I wished that I learned more javascript earlier in my career and have advised my kids that they should pay attention at school when sines and cosines are mentioned.

To get a point on Google Maps it seemed I would have to convert the OS Grid Ref to a Latitude/Longitude and plot that on the map. Almost correct except that this would plot an OS Lat/Lon against a WGS84 map, so the OS Lat/Long needs to be converted to a WGS84 Lat/Lon first.

Luckily the conversions were already available (thanks to Chris Veness) on the web but I couldn't locate a VB based version. So I decided to write my own VB.Net stuff to do the conversions so I could do these calculations server-side rather than in Javascript. To do this I converted Chris Veness' scripts which I barely understood. Once I'd written the converter I used it to bulk convert all property easting/northings and save the WGS84 Lat/Long in the databases, meaning that I'd never have to convert them again.

Caveat: I've tested this for various properties within my organisation's local area and it works for me. If it works for you too, great, if not then let me know where the errors are.

I created a .Net class project named GISFunctions.

How To Test The Class Project

    Dim strEasting as String = "462533"
    Dim strNorthing as String = "104496"
    Dim objEastingNorthing As New GISFunctions.EastingNorthing(strEasting & strNorthing)
    ' create OS Lat/Long
    Dim latlngOS As GISFunctions.LatitudeLongitude
    latlngOS = GISFunctions.Conversions.OSGridToLatLong(dr("Xref").ToString & dr("Yref").ToString)
    ' create WGS84 Lat/Long
    Dim latlngWGS84 As GISFunctions.LatitudeLongitude
    latlngWGS84 = GISFunctions.Conversions.OSGB36toWGS84(latlngOS)

The Code For The Class Project

Public Class EastingNorthing
    Private _easting As String
    Private _northing As String
    Private _gridSquare As String
    Private _gridReference As String
    Private _gridReferenceDigitsOnly As String
    Public Sub New(ByVal reference As String)
        _gridReference = reference
        ConvertGridReference(reference)
    End Sub

    Public ReadOnly Property Easting() As String
        Get
            Return _easting
        End Get
    End Property

    Public ReadOnly Property Northing() As String
        Get
            Return _northing
        End Get
    End Property

    Public ReadOnly Property GridSquare() As String
        Get
            Return _gridSquare
        End Get
    End Property

    Public ReadOnly Property EastingNorthing() As String
        Get
            Return _easting & _northing
        End Get
    End Property

    Public ReadOnly Property GridReference() As String
        Get
            Return _gridReference
        End Get
    End Property

    Public ReadOnly Property GridReferenceDigitsOnly() As String
        Get
            Return _gridReferenceDigitsOnly
        End Get
    End Property

    Private Sub ConvertGridReference(ByVal reference As String)
        If IsNumeric(reference.Substring(0, 2)) = False Then
            _gridSquare = reference.Substring(0, 2)
            _easting = reference.Substring(2, (reference.Length - 2) / 2)
            _northing = reference.Substring(2 + ((reference.Length - 2) / 2))
            Dim tmpString As String = Conversions.ConvertOSGridLetterToNumber(_gridSquare)
            _easting = tmpString.Substring(0, 1) + _easting
            _northing = tmpString.Substring(1, 1) + _northing
        Else
            _easting = reference.Substring(0, reference.Length / 2)
            _northing = reference.Substring(reference.Length / 2)
            _gridSquare = Conversions.ConvertOSGridNumberToLetter(_easting, _northing)
        End If
        _gridReferenceDigitsOnly = _easting.ToString & _northing.ToString
        _gridReference = _gridSquare & _easting.Substring(1) & _northing.Substring(1)
    End Sub
End Class


''' <summary>
''' A class which contains latitude, longitude and height
''' </summary>
''' <remarks></remarks>
Public Class LatitudeLongitude
    Private _latitude As Double
    Private _longitude As Double
    Private _height As Double
    Public Sub New(ByVal latitude As Double, ByVal longitude As Double, ByVal height As Double)
        _latitude = latitude
        _longitude = longitude
        _height = height
    End Sub

    Public Property latitude() As Double
        Get
            Return _latitude
        End Get
        Set(ByVal value As Double)
            _latitude = value
        End Set
    End Property

    Public Property longitude() As Double
        Get
            Return _longitude
        End Get
        Set(ByVal value As Double)
            _longitude = value
        End Set
    End Property

    Public Property height() As Double
        Get
            Return _height
        End Get
        Set(ByVal value As Double)
            _height = value
        End Set
    End Property
End Class


''' <summary>
''' EllipseParameter contains ellipse definition parameters used for conversion
''' </summary>
''' <remarks></remarks>
Public Class EllipseParameter
    Private _semiMajorAxisA As Double
    Private _semiMinorAxisB As Double
    Private _f As Double
    Public Sub New(ByVal EllipseType As String)
        Select Case LCase(EllipseType)
            Case "wgs84"
                _semiMajorAxisA = 6378137
                _semiMinorAxisB = 6356752.3142
                _f = 1 / 298.257223563
            Case "airy1830"
                _semiMajorAxisA = 6377563.396
                _semiMinorAxisB = 6356256.91
                _f = 1 / 299.3249646
            Case Else
                Throw New Exception("EllipseType should be either WGS84 or Airy1830")
        End Select
    End Sub

    Public ReadOnly Property SemiMajorAxisA() As Double
        Get
            Return _semiMajorAxisA
        End Get
    End Property

    Public ReadOnly Property SemiMinorAxisB() As Double
        Get
            Return _semiMinorAxisB
        End Get
    End Property

    Public ReadOnly Property f() As Double
        Get
            Return _f
        End Get
    End Property
End Class


''' <summary>
''' ConversionParameter is a class of parameters used when converting from OSGB36 to WGS84 and vice versa
''' </summary>
''' <remarks></remarks>
Public Class ConversionParameter
    Private _tx As Double
    Private _ty As Double
    Private _tz As Double
    Private _rx As Double
    Private _ry As Double
    Private _rz As Double
    Private _s As Double
    Private _conversionType As String

    Public Sub New(ByVal ConversionType As String)
        _conversionType = ConversionType
        Select Case LCase(ConversionType)
            Case "osgb36towgs84"
                _tx = 446.448
                _ty = -125.157
                _tz = 542.06
                _rx = 0.1502
                _ry = 0.247
                _rz = 0.8421
                _s = -20.4894
            Case "wgs84toosgb36"
                _tx = -446.448
                _ty = 125.157
                _tz = -542.06
                _rx = -0.1502
                _ry = -0.247
                _rz = -0.8421
                _s = 20.4894
            Case Else
                Throw New Exception("ConversionType should be either WGS84toOSGB36 or OSGB36toWGS84")
        End Select
    End Sub

    Public ReadOnly Property ConversionType() As String
        Get
            Return _conversionType
        End Get
    End Property

    Public ReadOnly Property tx() As Double
        Get
            Return _tx
        End Get
    End Property

    Public ReadOnly Property ty() As Double
        Get
            Return _ty
        End Get
    End Property

    Public ReadOnly Property tz() As Double
        Get
            Return _tz
        End Get
    End Property

    Public ReadOnly Property rx() As Double
        Get
            Return _rx
        End Get
    End Property

    Public ReadOnly Property ry() As Double
        Get
            Return _ry
        End Get
    End Property

    Public ReadOnly Property rz() As Double
        Get
            Return _rz
        End Get
    End Property

    Public ReadOnly Property s() As Double
        Get
            Return _s
        End Get
    End Property
End Class


Public Class Conversions
    Public Shared Function OSGB36toWGS84(ByVal p1) As LatitudeLongitude
        Dim latlngNew = ConvertCoordinates(p1, New EllipseParameter("Airy1830"), New ConversionParameter("OSGB36toWGS84"), _
            New EllipseParameter("WGS84"))
        Return latlngNew
    End Function

    Public Shared Function WGS84toOSGB36(ByVal p1) As LatitudeLongitude
        Dim latlngNew = ConvertCoordinates(p1, New EllipseParameter("WGS84"), New ConversionParameter("WGS84toOSGB36"), _
            New EllipseParameter("Airy1830"))
        Return latlngNew
    End Function

    Private Shared Function ConvertCoordinates(ByVal OriginalLatLong As LatitudeLongitude, ByVal EllipseParamSource As EllipseParameter, _
        ByVal ConversionParam As ConversionParameter, ByVal ElliseParamTarget As EllipseParameter)
        ' convert polar to cartesian coordinates (using source ellipse)
        Dim dblOriginalLatitudeInRadians As Double = ConvertDegreesToRadians(OriginalLatLong.latitude)
        Dim dblOriginalLongitudeInRadians As Double = ConvertDegreesToRadians(OriginalLatLong.longitude)
        Dim dblTempSemiMajorAxis As Double = EllipseParamSource.SemiMajorAxisA
        Dim dblTempSemiMinorAxis As Double = EllipseParamSource.SemiMinorAxisB
        Dim f As Double = EllipseParamSource.f
        Dim dblsinPhi = Math.Sin(dblOriginalLatitudeInRadians)
        Dim dblcosPhi = Math.Cos(dblOriginalLatitudeInRadians)
        Dim dblsinLambda = Math.Sin(dblOriginalLongitudeInRadians)
        Dim dblcosLambda = Math.Cos(dblOriginalLongitudeInRadians)
        Dim H As Integer = 1 'p1.height
        Dim eSq = (dblTempSemiMajorAxis * dblTempSemiMajorAxis - dblTempSemiMinorAxis * dblTempSemiMinorAxis) / (dblTempSemiMajorAxis * dblTempSemiMajorAxis)
        Dim nu = dblTempSemiMajorAxis / Math.Sqrt(1 - eSq * dblsinPhi * dblsinPhi)
        Dim x1 = (nu + H) * dblcosPhi * dblcosLambda
        Dim y1 = (nu + H) * dblcosPhi * dblsinLambda
        Dim z1 = ((1 - eSq) * nu + H) * dblsinPhi
        ' apply helmert transform using appropriate params
        Dim rx As Double = ConversionParam.rx / 3600 * Math.PI / 180  ' normalise seconds to radians
        Dim ry As Double = ConversionParam.ry / 3600 * Math.PI / 180
        Dim rz As Double = ConversionParam.rz / 3600 * Math.PI / 180
        Dim s1 As Double = ConversionParam.s / 1000000.0 + 1              ' normalise ppm to (s+1)
        ' apply transform
        Dim x2 As Double = ConversionParam.tx + x1 * s1 - y1 * rz + z1 * ry
        Dim y2 As Double = ConversionParam.ty + x1 * rz + y1 * s1 - z1 * rx
        Dim z2 As Double = ConversionParam.tz - x1 * ry + y1 * rx + z1 * s1
        ' convert cartesian to polar coordinates (using target ellipse)
        dblTempSemiMajorAxis = ElliseParamTarget.SemiMajorAxisA
        dblTempSemiMinorAxis = ElliseParamTarget.SemiMinorAxisB
        Dim precision = 4 / dblTempSemiMajorAxis  ' results accurate to around 4 metres
        eSq = (dblTempSemiMajorAxis * dblTempSemiMajorAxis - dblTempSemiMinorAxis * dblTempSemiMinorAxis) / _
            (dblTempSemiMajorAxis * dblTempSemiMajorAxis)
        Dim p = Math.Sqrt(x2 * x2 + y2 * y2)
        Dim phi = Math.Atan2(z2, p * (1 - eSq)), phiP = 2 * Math.PI
        Do While (Math.Abs(phi - phiP) > precision)
            nu = dblTempSemiMajorAxis / Math.Sqrt(1 - eSq * Math.Sin(phi) * Math.Sin(phi))
            phiP = phi
            phi = Math.Atan2(z2 + eSq * nu * Math.Sin(phi), p)
        Loop
        Dim lambda = Math.Atan2(y2, x2)
        H = p / Math.Cos(phi) - nu
        Dim LatLng As New LatitudeLongitude(ConvertRadiansToDegrees(phi), ConvertRadiansToDegrees(lambda), H)
        Return LatLng
    End Function

    ''' <summary>
    ''' convert geodesic co-ordinates to OS grid reference
    ''' </summary>
    ''' <param name="OriginalLatLong"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Shared Function LatLongToOSGrid(ByVal OriginalLatLong As LatitudeLongitude)
        Dim lat = ConvertDegreesToRadians(OriginalLatLong.latitude)
        Dim lon = ConvertDegreesToRadians(OriginalLatLong.longitude)
        Dim dblSemiMajorAxisA As Double = 6377563.396
        Dim dblSemiMinorAxisB As Double = 6356256.91          ' Airy 1830 major & minor semi-axes
        Dim dboScaleFactorF0 As Double = 0.9996012717                         ' NatGrid scale factor on central meridian
        Dim dblTrueOriginLat0 As Double = ConvertDegreesToRadians(49)
        Dim dblTrueOriginLon0 As Double = ConvertDegreesToRadians(-2)  ' NatGrid true origin
        Dim dblMapCoordTrueOriginN0 As Double = -100000
        Dim dblMapCoordTrueOriginE0 As Double = 400000        ' northing & easting of true origin, metres
        Dim dblEccentricitySquared As Double = 1 - (dblSemiMinorAxisB * dblSemiMinorAxisB) / (dblSemiMajorAxisA * dblSemiMajorAxisA)  ' eccentricity squared
        Dim n As Double = (dblSemiMajorAxisA - dblSemiMinorAxisB) / (dblSemiMajorAxisA + dblSemiMinorAxisB), n2 = n * n, n3 = n * n * n
        Dim cosLat As Double = Math.Cos(lat), sinLat = Math.Sin(lat)
        Dim nu As Double = dblSemiMajorAxisA * dboScaleFactorF0 / Math.Sqrt(1 - dblEccentricitySquared * sinLat * sinLat)              ' transverse radius of curvature
        Dim rho As Double = dblSemiMajorAxisA * dboScaleFactorF0 * (1 - dblEccentricitySquared) / Math.Pow(1 - dblEccentricitySquared * sinLat * sinLat, 1.5)  ' meridional radius of curvature
        Dim eta2 As Double = nu / rho - 1
        Dim Ma As Double = (1 + n + (5 / 4) * n2 + (5 / 4) * n3) * (lat - dblTrueOriginLat0)
        Dim Mb As Double = (3 * n + 3 * n * n + (21 / 8) * n3) * Math.Sin(lat - dblTrueOriginLat0) * Math.Cos(lat + dblTrueOriginLat0)
        Dim Mc As Double = ((15 / 8) * n2 + (15 / 8) * n3) * Math.Sin(2 * (lat - dblTrueOriginLat0)) * Math.Cos(2 * (lat + dblTrueOriginLat0))
        Dim Md As Double = (35 / 24) * n3 * Math.Sin(3 * (lat - dblTrueOriginLat0)) * Math.Cos(3 * (lat + dblTrueOriginLat0))
        Dim M As Double = dblSemiMinorAxisB * dboScaleFactorF0 * (Ma - Mb + Mc - Md)              ' meridional arc
        Dim cos3lat As Double = cosLat * cosLat * cosLat
        Dim cos5lat As Double = cos3lat * cosLat * cosLat
        Dim tan2lat As Double = Math.Tan(lat) * Math.Tan(lat)
        Dim tan4lat As Double = tan2lat * tan2lat
        Dim I As Double = M + dblMapCoordTrueOriginN0
        Dim II As Double = (nu / 2) * sinLat * cosLat
        Dim III As Double = (nu / 24) * sinLat * cos3lat * (5 - tan2lat + 9 * eta2)
        Dim IIIA As Double = (nu / 720) * sinLat * cos5lat * (61 - 58 * tan2lat + tan4lat)
        Dim IV As Double = nu * cosLat
        Dim V As Double = (nu / 6) * cos3lat * (nu / rho - tan2lat)
        Dim VI As Double = (nu / 120) * cos5lat * (5 - 18 * tan2lat + tan4lat + 14 * eta2 - 58 * tan2lat * eta2)
        Dim dLon As Double = lon - dblTrueOriginLon0
        Dim dLon2 As Double = dLon * dLon
        Dim dLon3 As Double = dLon2 * dLon
        Dim dLon4 As Double = dLon3 * dLon
        Dim dLon5 As Double = dLon4 * dLon
        Dim dLon6 As Double = dLon5 * dLon
        Dim dblNorthing = I + II * dLon2 + III * dLon4 + IIIA * dLon6
        Dim dblEasting = dblMapCoordTrueOriginE0 + IV * dLon + V * dLon3 + VI * dLon5
        Return ConvertOSGridNumberToLetter(dblEasting, dblNorthing) ', 8)
    End Function

    ''' <summary>
    ''' Convert OS Grid Reference to a Latitude/Longitude
    ''' </summary>
    ''' <param name="GridReference">An OS Grid Reference in the form TG1234512345 or 612345312345</param>
    ''' <returns>A LatitudeLongitude class containing the OS Lat/Long</returns>
    ''' <remarks>For a lat/long that can be used in GoogleMaps or VirtualEarth the returned class needs to be converted using ConvertOSGB36toWGS84</remarks>
    Public Shared Function OSGridToLatLong(ByVal GridReference As String) As LatitudeLongitude
        ' convert grid reference to Easting/Northing
        Dim gr As New EastingNorthing(GridReference)
        ' convert Easting/Northing to Latitude/Longitude
        Return OSGridToLatLong(gr.Easting, gr.Northing)
    End Function

    ''' <summary>
    ''' Convert OS Grid Easting and Northing to Latitude/Longitude
    ''' </summary>
    ''' <param name="Easting">An OS Grid Reference Easting in the form nnn or nnnnnn etc.</param>
    ''' <param name="Northing">An OS Grid Reference Northing in the form nnn or nnnnnn etc.</param>
    ''' <returns>A LatitudeLongitude class containing the OS Lat/Long</returns>
    ''' <remarks></remarks>
    Public Shared Function OSGridToLatLong(ByVal Easting As Double, ByVal Northing As Double) As LatitudeLongitude
        Dim dblEasting As Double = CType((Easting.ToString & "000000").Substring(0, 6), Double)
        Dim dblNorthing As Double = CType((Northing.ToString & "000000").Substring(0, 6), Double)
        Dim dblSemiMajorAxisA As Double = 6377563.396               ' Airy 1830 major semi-axis
        Dim dblSemiMinorAxisB As Double = 6356256.91                ' Airy 1830 minor semi-axis
        Dim dboScaleFactorF0 As Double = 0.9996012717               ' NatGrid scale factor on central meridian
        Dim dblTrueOriginLat0 As Double = 49 * Math.PI / 180        ' NatGrid true origin latitude
        Dim dblTrueOriginLon0 As Double = -2 * Math.PI / 180        ' NatGrid true origin longitude
        Dim dblMapCoordTrueOriginN0 As Double = -100000             ' northing of true origin, metres
        Dim dblMapCoordTrueOriginE0 As Double = 400000              ' easting of true origin, metres
        Dim dblEccentricitySquared As Double = 1 - (dblSemiMinorAxisB * dblSemiMinorAxisB) / (dblSemiMajorAxisA * dblSemiMajorAxisA) ' eccentricity squared
        Dim n As Double = (dblSemiMajorAxisA - dblSemiMinorAxisB) / (dblSemiMajorAxisA + dblSemiMinorAxisB)
        Dim n2 As Double = n * n
        Dim n3 As Double = n * n * n
        Dim lat As Double = dblTrueOriginLat0
        Dim M As Double = 0
        Do
            lat = (dblNorthing - dblMapCoordTrueOriginN0 - M) / (dblSemiMajorAxisA * dboScaleFactorF0) + lat
            Dim Ma As Double = (1 + n + (5 / 4) * n2 + (5 / 4) * n3) * (lat - dblTrueOriginLat0)
            Dim Mb As Double = (3 * n + 3 * n * n + (21 / 8) * n3) * Math.Sin(lat - dblTrueOriginLat0) * Math.Cos(lat + dblTrueOriginLat0)
            Dim Mc As Double = ((15 / 8) * n2 + (15 / 8) * n3) * Math.Sin(2 * (lat - dblTrueOriginLat0)) * Math.Cos(2 * (lat + dblTrueOriginLat0))
            Dim Md As Double = (35 / 24) * n3 * Math.Sin(3 * (lat - dblTrueOriginLat0)) * Math.Cos(3 * (lat + dblTrueOriginLat0))
            M = dblSemiMinorAxisB * dboScaleFactorF0 * (Ma - Mb + Mc - Md)                ' meridional arc
        Loop While (dblNorthing - dblMapCoordTrueOriginN0 - M >= 0.00001)  ' ie until < 0.01mm
        Dim sinLat As Double = Math.Sin(lat)
        Dim nu As Double = dblSemiMajorAxisA * dboScaleFactorF0 / Math.Sqrt(1 - dblEccentricitySquared * sinLat * sinLat)                     ' transverse radius of curvature
        Dim rho As Double = dblSemiMajorAxisA * dboScaleFactorF0 * (1 - dblEccentricitySquared) / Math.Pow(1 - dblEccentricitySquared * sinLat * sinLat, 1.5)     ' meridional radius of curvature
        Dim eta2 As Double = nu / rho - 1
        Dim tanLat As Double = Math.Tan(lat)
        Dim tan2lat As Double = tanLat * tanLat
        Dim tan4lat As Double = tan2lat * tan2lat
        Dim tan6lat As Double = tan4lat * tan2lat
        Dim secLat As Double = 1 / Math.Cos(lat)
        Dim VII As Double = tanLat / (2 * rho * nu)
        Dim VIII As Double = tanLat / (24 * rho * (nu ^ 3)) * (5 + 3 * (tanLat ^ 2) + eta2 - 9 * (tanLat ^ 2) * eta2)
        Dim IX As Double = tanLat / (720 * rho * (nu ^ 5)) * (61 + 90 * (tanLat ^ 2) + 45 * (tanLat ^ 4))
        Dim X As Double = secLat / nu
        Dim XI As Double = secLat / (6 * (nu ^ 3)) * (nu / rho + 2 * (tanLat ^ 2))
        Dim XII As Double = secLat / (120 * (nu ^ 5)) * (5 + 28 * (tanLat ^ 2) + 24 * (tanLat ^ 4))
        Dim XIIA As Double = secLat / (5040 * (nu ^ 7)) * (61 + 662 * (tanLat ^ 2) + 1320 * (tanLat ^ 4) + 720 * (tanLat ^ 6))
        Dim dE As Double = (dblEasting - dblMapCoordTrueOriginE0)
        lat = lat - VII * (dE ^ 2) + VIII * (dE ^ 4) - IX * (dE ^ 6)
        Dim lon As Double = dblTrueOriginLon0 + X * dE - XI * (dE ^ 3) + XII * (dE ^ 5) - XIIA * (dE ^ 7)
        Return New LatitudeLongitude(ConvertRadiansToDegrees(lat), ConvertRadiansToDegrees(lon), 0)
    End Function

    ''' <summary>Convert standard grid reference ('SU6253304496') to fully numeric ref ([462533,104496])</summary>
    ''' <param name="GridReference">Ordnance Survey grid reference in the form SUen, e.g. SU123123 or TG1234512345</param>
    ''' <returns>The converted grid offset digits, e.g. SU123123 returns 41 (4 squares along, 1 grid up)</returns>
    Public Shared Function ConvertOSGridLetterToNumber(ByVal GridReference As String) As String
        ' get numeric values of letter references, mapping A->0, B->1, C->2, etc:
        Dim intFirstChar As Integer = Asc(GridReference.ToUpper.Chars(0)) - 65
        Dim intSecondChar As Integer = Asc(GridReference.ToUpper.Chars(1)) - 65
        ' shuffle down letters after 'I' since 'I' is not used in grid:
        If (intFirstChar > 7) Then intFirstChar -= 1
        If (intSecondChar > 7) Then intSecondChar -= 1
        ' convert grid letters into 100km-square indexes from false origin (grid square SV):
        Dim intEasting As Integer = ((intFirstChar - 2) Mod 5) * 5 + (intSecondChar Mod 5)
        Dim intNorthing As Integer = (19 - Math.Floor(intFirstChar / 5) * 5) - Math.Floor(intSecondChar / 5)
        Return intEasting.ToString & intNorthing.ToString
    End Function

    ''' <summary>Convert fully numeric ref ([462533104496]) to standard grid reference ('SU6253304496')</summary>
    ''' <param name="Easting">Ordnance Survey Easting in the form n,  e.g. 123 or 12345</param>
    ''' <param name="Northing">Ordnance Survey Northing in the form n, e.g. 123 or 12345</param>
    ''' <returns>The grid square equivalent of the input parameters, e.g. 412345 and 1e.g. SU123123 returns 41 (4 squares along, 1 grid up)</returns>
    Public Shared Function ConvertOSGridNumberToLetter(ByVal Easting As Integer, ByVal Northing As Integer) As String
        ' get the 100km-grid indices
        Dim e100k As Integer = Math.Floor(Easting / 100000)
        Dim n100k As Integer = Math.Floor(Northing / 100000)
        If (e100k < 0 OrElse e100k > 6 OrElse n100k < 0 OrElse n100k > 12) Then Return ""
        ' translate those into numeric equivalents of the grid letters
        Dim intFirstChar = (19 - n100k) - (19 - n100k) Mod 5 + Math.Floor((e100k + 10) / 5)
        Dim intSecondChar = (19 - n100k) * 5 Mod 25 + e100k Mod 5
        ' compensate for skipped 'I' and build grid letter-pairs
        If (intFirstChar > 7) Then intFirstChar += 1
        If (intSecondChar > 7) Then intSecondChar += 1
        Return Chr(intFirstChar + 65) + Chr(intSecondChar + 65) ' String.fromCharCode(l1+'A'.charCodeAt(0), l2+'A'.charCodeAt(0))
    End Function

    ''' <summary>
    ''' convert degrees to radians
    ''' </summary>
    ''' <param name="Degrees">Value to be converted</param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Shared Function ConvertDegreesToRadians(ByVal Degrees As Double) As Double
        Return Degrees * Math.PI / 180
    End Function

    ''' <summary>
    ''' convert radians to degrees (signed)
    ''' </summary>
    ''' <param name="Radians">Value to be converted</param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Shared Function ConvertRadiansToDegrees(ByVal Radians As Double) As Double
        Return Radians * 180 / Math.PI
    End Function

    ''' <summary>
    ''' Return co-ordinate in degrees, minutes and seconds from a decimal value
    ''' </summary>
    ''' <param name="NumericDegrees">Numeric degrees</param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Shared Function ConvertNumericDegreesToDegreesMinutesSeconds(ByVal NumericDegrees As Double) As String
        Dim dblTemp As Double = Math.Abs(NumericDegrees)
        dblTemp += 1 / 720000000  ' add .005 millisecond for rounding
        Dim deg As Double = Math.Floor(dblTemp)
        Dim min As Double = Math.Floor((dblTemp - deg) * 60)
        Dim sec As Double = Math.Floor((dblTemp - deg - min / 60) * 3600)
        Return deg.ToString("000") + Chr("176") + min.ToString("00") + "'" + sec.ToString("0") + """"
    End Function

    ''' <summary>
    ''' convert numeric degrees to deg/min/sec longitude
    ''' </summary>
    ''' <param name="NumericDegrees">Numeric degrees</param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Shared Function ConvertNumericDegreesToLatitude(ByVal NumericDegrees As Double) As String
        Dim strTemp As String = ""
        strTemp = ConvertNumericDegreesToDegreesMinutesSeconds(NumericDegrees).Substring(1)
        If NumericDegrees < 0 Then strTemp &= "S" Else strTemp &= "N"
        Return strTemp
    End Function

    ''' <summary>
    ''' convert numeric degrees to deg/min/sec longitude
    ''' </summary>
    ''' <param name="NumericDegrees">Numeric degrees</param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Shared Function ConvertNumericDegreesToLongitude(ByVal NumericDegrees As Double) As Double
        Dim strTemp As String = ""
        strTemp = ConvertNumericDegreesToDegreesMinutesSeconds(NumericDegrees)
        If NumericDegrees < 0 Then strTemp &= "W" Else strTemp &= "E"
        Return strTemp
    End Function
End Class