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)!

Did this help you out? Please consider making a donation:

Original Content , , , , , ,

17 Responses to Convert LAMBERT 1972 to WGS84

  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

    • 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.

      • 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 ?

    • 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

      Now about the latitudes
      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.

  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!

  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,
    Thank you for your tutorial,
    Please, can you add the code to convert WGS84 to Lambert 72 in java

    Thanks

Leave a Reply

Your email address will not be published. Required fields are marked *