Algorithms used to convert from Belgian Lambert 1972 to WGS84 and vice versa

Conversion from Belgian Lambert 1972 to spherical coordinates and vice versa are inspired from the following paper: "Syst�mes de r�f�rence et formules de transformation en usage en Belgique".  Institut G�ographique National, 1989, 48 pp. 

Spherical coordinates involved in the algorithms are always expressed in the Belgium datum.  In order to transform spherical coordinates from one datum to another, I use the Molodensky equation (3 parameters; see below). It is not the most accurate way to convert between datums but the precision is quite enough for what I need. Molodensky equation can be found on different web sites (eg here, or here, ...).

The algorithms given below are in vb.net language, the code is not optimized; it is provided "as is", without warranty of any kind, either expressed or implied. The code is distributed as GNU General Public License.

The usual way to convert is the following:

Conversion process

So if you want to convert from Belgian coordinates to WGS84, you have to do it in 2 steps:

  1. use the conversion from Belgian to spherical coordinates
  2. use the Molodensky equation to transform Belgian spherical coordinates to WGS84

On the other hand, if you want to convert from WGS84 to Belgian coordinates, you have to :

  1. use the Molodensky equation to transform WGS84 to Belgian spherical coordinates
  2. use the conversion from Belgian spherical coordinates to Belgian cartesian coordinates

Belgian Lambert 72 to spherical coordinates

    
        '
        '       Belgian Lambert 1972---> Spherical coordinates
        '       Input parameters : X, Y = Belgian coordinates in meters
        '       Output : latitude and longitude in Belgium Datum!
        '

        Const LongRef As Double = 0.076042943        '=4�21'24"983
        Const nLamb As Double = 0.7716421928
        Const aCarre As Double = 6378388 ^ 2
        Const bLamb As Double = 6378388 * (1 - (1 / 297))
        Const eCarre As Double = (aCarre - bLamb ^ 2) / aCarre
        Const KLamb As Double = 11565915.812935

        Dim eLamb As Double = Math.Sqrt(eCarre)
        Dim eSur2 As Double = eLamb / 2

        Dim Tan1 As Double = (X - 150000.01256) / (5400088.4378 - Y)
        Dim Lambda As Double = LongRef + (1 / nLamb) * (0.000142043 + Math.Atan(Tan1))
        Dim RLamb As Double = Math.Sqrt((X - 150000.01256) ^ 2 + (5400088.4378 - Y) ^ 2)

        Dim TanZDemi As Double = (RLamb / KLamb) ^ (1 / nLamb)
        Dim Lati1 As Double = 2 * Math.Atan(TanZDemi)

        Dim eSin As Double
        Dim Mult1, Mult2, Mult As Double
        Dim LatiN, Diff As Double

        Dim lat, lng As Double

        Do
            eSin = eLamb * Math.Sin(Lati1)
            Mult1 = 1 - eSin
            Mult2 = 1 + eSin
            Mult = (Mult1 / Mult2) ^ (eLamb / 2)
            LatiN = (Math.PI / 2) - (2 * (Math.Atan(TanZDemi * Mult)))
            Diff = LatiN - Lati1
            Lati1 = LatiN
        Loop While Math.Abs(Diff) > 0.0000000277777

        lat = (LatiN * 180) / Math.PI
        lng = (Lambda * 180) / Math.PI
        
        

Spherical coordinates to Belgian Lambert 72

    
     '
        '       Conversion from spherical coordinates to Lambert 72
        '       Input parameters : lat, lng (spherical coordinates)
        '       Spherical coordinates are in decimal degrees converted to Belgium datum!
        '

        Const LongRef As Double = 0.076042943        '=4�21'24"983
        Const bLamb As Double = 6378388 * (1 - (1 / 297))
        Const aCarre As Double = 6378388 ^ 2
        Const eCarre As Double = (aCarre - bLamb ^ 2) / aCarre
        Const KLamb As Double = 11565915.812935
        Const nLamb As Double = 0.7716421928

        Dim eLamb As Double = Math.Sqrt(eCarre)
        Dim eSur2 As Double = eLamb / 2

        'conversion to radians
        lat = (Math.PI / 180) * lat
        lng = (Math.PI / 180) * lng

        Dim eSinLatitude As Double = eLamb * Math.Sin(lat)
        Dim TanZDemi As Double = (Math.Tan((Math.PI / 4) - (lat / 2))) * _  
            (((1 + (eSinLatitude)) / (1 - (eSinLatitude))) ^ (eSur2))

        Dim RLamb As Double = KLamb * ((TanZDemi) ^ nLamb)

        Dim Teta As Double = nLamb * (lng - LongRef)

        Dim x, y As Single

        x = 150000 + 0.01256 + RLamb * Math.Sin(Teta - 0.000142043)
        y = 5400000 + 88.4378 - RLamb * Math.Cos(Teta - 0.000142043)

     

Belgian Datum to WGS84 conversion (Molodensky 3 parameters)

        '
        'Input parameters : Lat, Lng : latitude / longitude in decimal degrees and in Belgian 1972 datum
        'Output parameters : LatWGS84, LngWGS84 : latitude / longitude in decimal degrees and in WGS84 datum
        '

        Const Haut = 0      'Altitude
        Dim LatWGS84, LngWGS84 As Double
        Dim DLat, DLng As Double
        Dim Dh As Double
        Dim dy, dx, dz As Double
        Dim da, df As Double
        Dim LWa, Rm, Rn, LWb As Double
        Dim LWf, LWe2 As Double
        Dim SinLat, SinLng As Double
        Dim CoSinLat As Double
        Dim CoSinLng As Double

        Dim Adb As Double

        'conversion to radians
        Lat = (Math.PI / 180) * Lat
        Lng = (Math.PI / 180) * Lng

        SinLat = Math.Sin(Lat)
        SinLng = Math.Sin(Lng)
        CoSinLat = Math.Cos(Lat)
        CoSinLng = Math.Cos(Lng)

        dx = -125.8
        dy = 79.9
        dz = -100.5
        da = -251.0
        df = -0.000014192702

        LWf = 1 / 297
        LWa = 6378388
        LWb = (1 - LWf) * LWa
        LWe2 = (2 * LWf) - (LWf * LWf)
        Adb = 1 / (1 - LWf)

        Rn = LWa / System.Math.Sqrt(1 - LWe2 * SinLat * SinLat)
        Rm = LWa * (1 - LWe2) / (1 - LWe2 * Lat * Lat) ^ 1.5

        DLat = -dx * SinLat * CoSinLng - dy * SinLat * SinLng + dz * CoSinLat
        DLat = DLat + da * (Rn * LWe2 * SinLat * CoSinLat) / LWa
        DLat = DLat + df * (Rm * Adb + Rn / Adb) * SinLat * CoSinLat
        DLat = DLat / (Rm + Haut)

        DLng = (-dx * SinLng + dy * CoSinLng) / ((Rn + Haut) * CoSinLat)
        Dh = dx * CoSinLat * CoSinLng + dy * CoSinLat * SinLng + dz * SinLat
        Dh = Dh - da * LWa / Rn + df * Rn * Lat * Lat / Adb

        LatWGS84 = ((Lat + DLat) * 180) / Math.PI
        LngWGS84 = ((Lng + DLng) * 180) / Math.PI
        
        

WGS84 to Belgian Datum conversion (Molodensky 3 parameters)

        '
        'Input parameters : Lat, Lng : latitude / longitude in decimal degrees and in WGS84 datum
        'Output parameters : LatBel, LngBel : latitude / longitude in decimal degrees and in Belgian datum
        '

        Const Haut = 0      'Altitude
        Dim LatBel, LngBel As Double
        Dim DLat, DLng As Double
        Dim Dh As Double
        Dim dy, dx, dz As Double
        Dim da, df As Double
        Dim LWa, Rm, Rn, LWb As Double
        Dim LWf, LWe2 As Double
        Dim SinLat, SinLng As Double
        Dim CoSinLat As Double
        Dim CoSinLng As Double

        Dim Adb As Double

        'conversion to radians
        Lat = (Math.PI / 180) * Lat
        Lng = (Math.PI / 180) * Lng

        SinLat = Math.Sin(Lat)
        SinLng = Math.Sin(Lng)
        CoSinLat = Math.Cos(Lat)
        CoSinLng = Math.Cos(Lng)

        dx = 125.8
        dy = -79.9
        dz = 100.5
        da = 251.0
        df = 0.000014192702

        LWf = 1 / 297
        LWa = 6378388
        LWb = (1 - LWf) * LWa
        LWe2 = (2 * LWf) - (LWf * LWf)
        Adb = 1 / (1 - LWf)

        Rn = LWa / System.Math.Sqrt(1 - LWe2 * SinLat * SinLat)
        Rm = LWa * (1 - LWe2) / (1 - LWe2 * Lat * Lat) ^ 1.5

        DLat = -dx * SinLat * CoSinLng - dy * SinLat * SinLng + dz * CoSinLat
        DLat = DLat + da * (Rn * LWe2 * SinLat * CoSinLat) / LWa
        DLat = DLat + df * (Rm * Adb + Rn / Adb) * SinLat * CoSinLat
        DLat = DLat / (Rm + Haut)

        DLng = (-dx * SinLng + dy * CoSinLng) / ((Rn + Haut) * CoSinLat)
        Dh = dx * CoSinLat * CoSinLng + dy * CoSinLat * SinLng + dz * SinLat
        Dh = Dh - da * LWa / Rn + df * Rn * Lat * Lat / Adb

        LatBel = ((Lat + DLat) * 180) / Math.PI
        LngBel = ((Lng + DLng) * 180) / Math.PI
        
Ancien site web ©Yvan Barbier, UMons, DEMNA