# Convert LAMBERT 1972 to WGS84 Needed to convert LAMBERT 1972 (EPSG:31370) coordinates to WGS84 (EPSG:4326) for a project I’m working on in order to use the coordinates with Google Maps (which uses a SphericalMercator (EPSG:900913) projection).

Took me quite a while to find the needed algorithm, so hereby I’m listing them (and am providing the JavaScript implementation):

JavaScript:

``````var lambert72toWGS84 = function(x, y){

var newLongitude, newLatitude;

var n = 0.77164219,
F = 1.81329763,
thetaFudge = 0.00014204,
e = 0.08199189,
a = 6378388,
xDiff = 149910,
yDiff = 5400150,
theta0 = 0.07604294;

var xReal = xDiff - x,
yReal = yDiff - y;

var rho = Math.sqrt(xReal * xReal + yReal * yReal),
theta = Math.atan(xReal / -yReal);

newLongitude = (theta0 + (theta + thetaFudge) / n) * 180 / Math.PI;
newLatitude = 0;

for (var i = 0; i < 5 ; ++i) {
newLatitude = (2 * Math.atan(Math.pow(F * a / rho, 1 / n) * Math.pow((1 + e * Math.sin(newLatitude)) / (1 - e * Math.sin(newLatitude)), e / 2))) - Math.PI / 2;
}
newLatitude *= 180 / Math.PI;
return [newLatitude, newLongitude];

}``````

Java:

``````public static Point lambert72toWGS84(double x, double y) {

double newLongitude;
double newLatitude;

double n = 0.77164219;
double F = 1.81329763;
double thetaFudge = 0.00014204;
double e = 0.08199189;
double a = 6378388;
double xDiff = 149910;
double yDiff = 5400150;

double theta0 = 0.07604294;

double xReal = xDiff - x;
double yReal = yDiff - y;

double rho = Math.sqrt(xReal * xReal + yReal * yReal);
double theta = Math.atan(xReal / -yReal);

newLongitude = (theta0 + (theta + thetaFudge) / n) * 180 / Math.PI;
newLatitude = 0;

for (int i = 0; i < 5 ; ++i) {
newLatitude = (2 * Math.atan(Math.pow(F * a / rho, 1 / n) * Math.pow((1 + e * Math.sin(newLatitude)) / (1 - e * Math.sin(newLatitude)), e / 2))) - Math.PI / 2;
}
newLatitude *= 180 / Math.PI;
return new Point(newLatitude, newLongitude);
}``````

C/C++:

``````void LambertToLongLat(double x,double y, double& longt, double& lat){
double n = 0.77164219;
double F = 1.81329763;
double thetaFudge = 0.00014204;
double e= 0.08199189;
double a= 6378388;
double xDiff= 150000;
double yDiff = 5400088.44;
double theta0 = 0.07604294;

double xReal = xDiff-x;
double yReal = yDiff-y;

double rho = sqrt(xReal*xReal + yReal * yReal);
double theta = atan(xReal/-yReal);

longt = (theta0 + (theta+thetaFudge) /n)*180/PI;
lat =0;
for(size_t i=0; i<5; ++i){
lat =(2*atan(pow(F*a/rho,1/n)* pow((1+ e*sin(lat))/(1-e*sin(lat)),e/2)))-PI/2;
}
lat *=180/PI;
}``````

Hope this helps!

Related: if you want to do other transforms (or flip the transorm), check out Transformation between datum’s in JavaScript. Thanks for the tip, Ibrahim (via e-mail)!

Consider donating.

I don't run ads on my blog nor do I do this for profit. A donation however would always put a smile on my face though. Thanks! Bramus is a frontend web developer from Belgium, working as a Chrome Developer Relations Engineer at Google. From the moment he discovered view-source at the age of 14 (way back in 1997), he fell in love with the web and has been tinkering with it ever since (more …)

## Join the Conversation

1. 2. 3. 4. 5. 6. 1. venu says:

Hi,
I’ve used your java version for the conversion but unfortunately i am not getting correct latitude and longitude . what is the package that you are using for the Poin class. Any help on this would be appreciated .

Venu

1. Bramus! says:

Point is just a basic class with a latitude and a longitude datamember (type: double). That class can easily be created, or you could let the function return an Array instead of a Point instance.

1. venu says:

Hi,
Thank you for your email. I’ve applied your algorithm for the below input and i am able to get the output but it isn’t coming as expected as we want to convert Lambert II Étendu Coordinate system (EPSG:27572) to WGS84. Is there any way that we can manipulate the above algorithm that can fit our requirement. Thanks.

Regards,
Venu

2. Wim Michiels says:

Dear Bram,

thanks for this wonderful tool!

I rewrote the conversion in VBA, useful if one has an access database with co-ordinates… The longitudinal works fine, but the lateral is a bit off… Mechelen ends up in Holland…

Am i doing something wrong here?

Cheers,

Wim

Option Compare Database

Function lambert72toWGS84_X(x, y) As Double

Dim newLongitude, newLatitude As Double
Dim n As Double
Dim F As Double
Dim thetaFudge As Double
Dim e As Double
Dim a As Double
Dim xDiff As Double
Dim yDiff As Double
Dim theta0 As Double
Dim xReal As Double
Dim yReal As Double
Dim PI As Double

n = 0.77164219
F = 1.81329763
thetaFudge = 0.00014204
e = 0.08199189
a = 6378388
xDiff = 149910
yDiff = 5400150
theta0 = 0.07604294

xReal = xDiff – x
yReal = yDiff – y
PI = 3.14159265358979

rho = Sqr((xReal * xReal) + (yReal * yReal))
theta = Atn(xReal / -yReal)
newLongitude = (theta0 + ((theta + thetaFudge) / n)) * 180 / PI
lambert72toWGS84_X = newLongitude

End Function

Function lambert72toWGS84_Y(x, y) As Double

Dim newLongitude, newLatitude As Double

Dim n As Double
Dim F As Double
Dim thetaFudge As Double
Dim e As Double
Dim a As Double
Dim xDiff As Double
Dim yDiff As Double
Dim theta0 As Double
Dim xReal As Double
Dim yReal As Double
Dim PI As Double
Dim rho As Double
Dim theta As Double

n = 0.77164219
F = 1.81329763
thetaFudge = 0.00014204
e = 0.08199189
a = 6378388
xDiff = 149910
yDiff = 5400150
theta0 = 0.07604294
PI = 3.141592653

xReal = xDiff – x
yReal = yDiff – y

rho = Sqr((xReal * xReal) + (yReal * yReal))
theta = Atn(xReal / -yReal)
newLatitude = 0
For i = 0 To 4
newLatitude = (2 * (Atn(((F * a) / rho) ^ (1 / n))) * (((1 + (e * Sin(newLatitude))) / (1 – (e * Sin(newLatitude)))) ^ (e / 2))) – (PI / 2)

Next i
newLatitude = newLatitude * (180 / PI)
lambert72toWGS84_Y = newLatitude

End Function

3. Joost says:

Thanks for sharing, this was really useful for me.

4. Phil De Troy says:

Hello

A first comment about the method for Latitude, you may always sart from an approximate latitude (0,87 for Southern Belgium, 0,89 for Flanders, it reduced the required iteration number and increases accurracy.
Second comment is that your C/C++ script uses slightly different values for xDiff and yDiff
Being involved in studies about management of the territory, source info is aways given in Lambert 72, and points on Google maps are geographical.
I found a page http://www.econet.ulg.ac.be/biogeonet/index.php?pg=951 using a similar script and results are tested and found compliant. But copy-paste method takes a lot of time for maps with more than 500 points to convert.
The formula as proposed by Wim, inserted in an Excel file, give the same result for longitude, difference 1.1 metre being probably due to the reference Google Spherical vs WGS84. But as he points the problem for latitude, I encountered the same, with a difference of 38.65 km, anyway when I tried to write it from Bram’s formulas, and the difference was higher.
Source Lambert 72 141332,132433 (somewhere at the end of Rue des Verreries, Manage). Biogeonet returns 50,502627°N – 004,246559°E.
Excel table returns 51,049445°N – 004,246574° E

Anybody can help ?

1. Phil De Troy says:

New information found on NGI Web site see http://www.ngi.be/Common/Lambert2008/Transformation_Geographic_Lambert_NL.pdf

Values for xDiff and yDiff in C/C++ script is the original, supposed to give co-ordinates as Hayford 1924 ellipsoïd.
Values used in the other scripts are adjusted to match WGS84 co-ordinates. I found better results, matching closely the values of Biogeonet with 149911,1 and 5400149,9

NGI suggest an initial value of NewLatitude.
Calculation for longitude are giving exactly the same results aus your scripts but with another formula not using thetaFudge.

For latitude, an intermediate variable is calculated, based on the initial values.
Using these formulas, I found the expected results.

5. aart says:

Thanks for the usefull code Bram.us en Wim.

@Phil De Troy, could you please post the code for the latitude calculation for those of us who are mathematically impaired and not able to translate the formulas into code.

1. Gelini says:

i encounter the same problem…also latitude is way off the expected value.
Does anyone has the code of Phill?

6. Hilde says:

I tried the code of Wim in Access as well and also remarked the problem with the latitude. I wondered if there exists un update of the code. Can I please get it ? I would be very greatful !

7. Jelle says:

Hello. Did anyone found a way to get the latitude working in Excel or Access? The longitude seems to be ok more or less, but the latitude is way off. Thank you very much!

1. Axel Van den Bosch says:

I’ve made an excel with VBA for converting multiple WGS84 coordinates to Lambert 72

(Thanks to: http://zoologie.umh.ac.be/tc/algorithms.aspx#c2 (Yvan Barbier, UMons, DEMNA))

1. John De Meyer says:

Axel,
I know it has been a while since you to post the link here. I am at the moment with a problem that I need the same conversion. Would you like to make the excel file back available or send it to mail.

Mvg John

8. Kjel Dupon says:

Hi Axel

Your excel file is very useful but I was wondering if there exists a file to convert Lambert 72 into WGS84, the other way around.

Thanks very much.

Kjel Dupon

9. Robin Mulkers says:

If you are looking after a tool to convert bunch of coordinates in a CSV file (or in a shapefile), you can always use QGIS. http://qgis.org
You have to import a csv file as a vector layer and specify it is in the Lambert 72 coordinate system (EPSG:31370) and then export it in the WGS84 coordinate system (EPSG:4326)

10. thomas says:

You are my goddamn hero.

11. Jan Albert says:

For the delphi/pascal users
Writes the coordinates to a Tpointf (single-precision).

function lambert72toWGS84(const x, y : double): Tpointf;
var XReal, Yreal, rho, theta : double;
i : integer;
const n = 0.77164219;
F = 1.81329763;
thetaFudge = 0.00014204;
e = 0.08199189;
a = 6378388;
xDiff = 149910;
yDiff = 5400150;
theta0 = 0.07604294;
begin
xReal := xDiff – x;
yReal := yDiff – y;

rho := sqrt(xReal * xReal + yReal * yReal);
theta := arctan(xReal/(yReal * -1));
result.X := (theta0 + (theta + thetaFudge)/n) * 180 / system.PI; //Newlongitude

result.Y := 0;
for i := 0 to 4 do
result.Y := (2 * arctan(power(F * a/rho, 1/n) * power((1 + e*sin(result.Y)) / (1-e * sin(result.Y)), e/2))) – System.PI/2;
result.Y := (result.Y * 180)/system.PI; //Newlatitude
end;

12. Yaz says:

Hi,