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:
So if you want to convert from Belgian coordinates to WGS84, you have to do it in 2 steps:
-
use the conversion from Belgian to spherical coordinates
-
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 :
-
use the Molodensky equation to transform WGS84 to Belgian spherical coordinates
-
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