[Données] Formule de conversion coordonnées

[arena]

Bien en reprenant le script de Bapt plus haut, j'obtiens ceci en php.

J'aimerai que quelqu'un puisse me valider la conversion avec un jeu d'essai
les paramètres doivent correspondre à du Lambert2 étendu vers WGS84.

nb : tous les algo retranscrits ne sont pas forcément utilisés

Merci de votre aide

<?php
function ALG0001($phi,$e)
{
    $temp = ( 1 - ( $e * sin( $phi ) ) ) / ( 1 + ( $e * sin( $phi ) ) );

    $L = log ( tan ( (pi()/4) + ($phi/2) ) * pow ($temp, ($e/2) ));

    return $L;
}


function ALG0002($L,$e,$epsilon)
{
    $phi[0] = 2 * atan(exp($L)) - (pi()/2);

    $i=0;
    do
    {
        $i++;
        $temp = ( 1 + ( $e * sin( $phi[$i-1] ) ) ) / ( 1 - ( $e * sin( $phi[$i-1] ) ) );
        $phi[$i] = 2 * atan ( pow ($temp, ($e/2)) * exp ($L) ) - pi()/2;
    }
    while (abs($phi[$i] - $phi[$i - 1]) >= $epsilon);

    return $phi[$i];
}

function ALG0004($X,$Y,$n,$c,$Xs,$Ys,$lambdac,$e,$epsilon)
{
    $R = sqrt( pow(($X - $Xs),2) + pow(($Y - $Ys),2) );
    $gamma = atan(($X - $Xs)/($Ys - $Y));

    $lambda = $lambdac + ($gamma / $n);

    $L = (-1 / $n) * log(abs($R/$c));

    $phi = ALG0002($L,$e,$epsilon);

    $coords['lambda']=$lambda;
    $coords['phi']=$phi;
    
    return $coords;
}

function ALG0009($lambda,$phi,$he,$a,$e)
{
    $N = ALG0021($phi,$a,$e);
    
    $X = ($N + $he) * cos($phi) * cos($lambda);

    $Y = ($N + $he) * cos($phi) * sin($lambda);

    $Z = ($N * (1 - $e*$e)  + $he) * sin ($phi);

    $coords['X']=$X;
    $coords['Y']=$Y;
    $coords['Z']=$Z;
    
    return $coords;
}


function ALG0012($X,$Y,$Z,$a,$e,$epsilon)
{
    $lambda = atan ($Y/$X);

    $P = sqrt($X*$X + $Y*$Y);
    $phi[0] = atan ($Z/ ($P * (1 - ( ($a*$e*$e)/sqrt($X*$X + $Y*$Y + $Z*$Z) ) ) ) );

    $i = 0;
    do
    {
        $i++;
        $temp =  pow((1 - ( $a * $e*$e * cos($phi[$i - 1] )/( $P * sqrt(1 - $e*$e * sin($phi[$i - 1])*sin($phi[$i - 1]))))),-1);
        $phi[$i] = atan( $temp * $Z / $P );
    }
    while (abs($phi[$i] - $phi[$i - 1]) >= $epsilon);

    $phix = $phi[$i];

    $he = ($P/cos($phix)) - ($a/sqrt(1 - $e*$e * sin($phix)*sin($phix)));

    $coords['lambda']=$lambda;
    $coords['phi']=$phix;
    $coords['he']=$he;

    return $coords;
}

function ALG0013($Tx,$Ty,$Tz,$D,$Rx,$Ry,$Rz,$U)
{
    $V['X'] = $Tx + $U['X'] * (1 + $D) + $U['Z'] * $Ry - $U['Y'] * $Rz;
    $V['Y'] = $Ty + $U['Y'] * (1 + $D) + $U['X'] * $Rz - $U['Z'] * $Rx;
    $V['Z'] = $Tz + $U['Z'] * (1 + $D) + $U['Y'] * $Rx - $U['X'] * $Ry;

    return $V;
}

function ALG0019($lambda0,$phi0,$k0,$X0,$Y0,$a,$e)
{
    $lambdac = $lambda0;
    $n = sin($phi0);
    $C = $k0 * ALG0021($phi0,$a,$e) * tan (pi()/2 - $phi0) * exp ( $n * ALG0001($phi0,$e) );
    $Xs = $X0;
    $Ys = $Y0 + $k0 * ALG0021($phi0,$a,$e) * tan (pi()/2 - $phi0) ;

    $tab ['e'] = $e;
    $tab ['n'] = $n;
    $tab ['C'] = $C;
    $tab ['lambdac'] = $lambdac;
    $tab ['Xs'] = $Xs;
    $tab ['Ys'] = $Ys;

    return $tab;

}

function ALG0021($phi,$a,$e)
{
    $N = $a/sqrt( 1 - $e * $e * sin($phi) * sin($phi) );
    return $N;

}

function ALG0054($lambda0,$phi0,$X0,$Y0,$phi1,$phi2,$a,$e)
{
    $lambdac = $lambda0;
    $n = ( (log( (ALG0021($phi2,$a,$e)*cos($phi2))/(ALG0021($phi1,$a,$e)*cos($phi1)) ))/(ALG0001($phi1,$e) - ALG0001($phi2,$e) ));
    $C = ((ALG0021($phi1,$a,$e)* cos($phi1))/$n) * exp($n * ALG0001($phi1,$e));


    if ($phi0 == (pi()/2))
    {
        $Xs = $X0;
        $Ys = $Y0;
    }
    else
    {
            echo ('coucou');
        $Xs = $X0;
        $Ys = $Y0 + $C * exp(-1 * $n * ALG0001($phi0,$e));
    }

    $tab ['e'] = $e;
    $tab ['n'] = $n;
    $tab ['C'] = $C;
    $tab ['lambdac'] = $lambdac;
    $tab ['Xs'] = $Xs;
    $tab ['Ys'] = $Ys;

    return $tab;

}

function Lambert2e2WGS84($X,$Y)
{
    $epsilon = 0.00000000001;
        
    $n = 0.7289686274;
    $c = 11745793.39;
    $Xs = 600000;
    $Ys = 8199695.768;
    $lambdac = 0; //0.04079234433 pour greenwich
    
    $e = 0.08248325676; //(première excentricité de l’ellipsoïde Clarke 1880 français)

    
    $coords = ALG0004($X,$Y,$n,$c,$Xs,$Ys,$lambdac,$e,$epsilon);

    $he = 100;
    $a = 6378249.2;
    
    $coords = ALG0009($coords['lambda'],$coords['phi'],$he,$a,$e);
    
    $Tx = -168;
    $Ty = -60;
    $Tz = +320;
    $D  = 0;
    $Rx = $Ry = $Rz = 0;
    $U  = $coords;
    
    $coords = ALG0013($Tx,$Ty,$Tz,$D,$Rx,$Ry,$Rz,$U);
    
    $a = 6378137.0;    // ellipsoïdes WGS84
    $f = 1/298.257223563;
    $b = $a*(1-$f);
    $e = sqrt(($a*$a - $b*$b)/($a*$a));

    $X = $coords['X'];
    $Y = $coords['Y'];
    $Z = $coords['Z'];
    $coords = ALG0012($X,$Y,$Z,$a,$e,$epsilon);
    
    return $coords;
}
//
//http://www.ign.fr/telechargement/MPro/geodesie/CIRCE/transfo.pdf
//http://www.ign.fr/telechargement/MPro/geodesie/CIRCE/NTG_80.pdf
//http://www.ign.fr/telechargement/MPro/geodesie/CIRCE/NTG_71.pdf
//
print_r(Lambert2e2WGS84(591647.56,2426659.65));
?>

[arena]

Bon et bien pour ceux que ca interesse :

code validé avec le méridien de greenwich

<?php
function ALG0001($phi,$e)
{
    $temp = ( 1 - ( $e * sin( $phi ) ) ) / ( 1 + ( $e * sin( $phi ) ) );

    $L = log ( tan ( (pi()/4) + ($phi/2) ) * pow ($temp, ($e/2) ));

    return $L;
}


function ALG0002($L,$e,$epsilon)
{
    $phi[0] = 2 * atan(exp($L)) - (pi()/2);

    $i=0;
    do
    {
        $i++;
        $temp = ( 1 + ( $e * sin( $phi[$i-1] ) ) ) / ( 1 - ( $e * sin( $phi[$i-1] ) ) );
        $phi[$i] = 2 * atan ( pow ($temp, ($e/2)) * exp ($L) ) - pi()/2;
    }
    while (abs($phi[$i] - $phi[$i - 1]) >= $epsilon);

    return $phi[$i];
}

function ALG0004($X,$Y,$n,$c,$Xs,$Ys,$lambdac,$e,$epsilon)
{
    $R = sqrt( pow(($X - $Xs),2) + pow(($Y - $Ys),2) );
    $gamma = atan(($X - $Xs)/($Ys - $Y));

    $lambda = $lambdac + ($gamma / $n);

    $L = (-1 / $n) * log(abs($R/$c));

    $phi = ALG0002($L,$e,$epsilon);

    $coords['lambda']=$lambda;
    $coords['phi']=$phi;
    
    return $coords;
}

function ALG0009($lambda,$phi,$he,$a,$e)
{
    $N = ALG0021($phi,$a,$e);
    
    $X = ($N + $he) * cos($phi) * cos($lambda);

    $Y = ($N + $he) * cos($phi) * sin($lambda);

    $Z = ($N * (1 - $e*$e)  + $he) * sin ($phi);

    $coords['X']=$X;
    $coords['Y']=$Y;
    $coords['Z']=$Z;
    
    return $coords;
}


function ALG0012($X,$Y,$Z,$a,$e,$epsilon)
{
    $lambda = atan ($Y/$X);

    $P = sqrt($X*$X + $Y*$Y);
    $phi[0] = atan ($Z/ ($P * (1 - ( ($a*$e*$e)/sqrt($X*$X + $Y*$Y + $Z*$Z) ) ) ) );

    $i = 0;
    do
    {
        $i++;
        $temp =  pow((1 - ( $a * $e*$e * cos($phi[$i - 1] )/( $P * sqrt(1 - $e*$e * sin($phi[$i - 1])*sin($phi[$i - 1]))))),-1);
        $phi[$i] = atan( $temp * $Z / $P );
    }
    while (abs($phi[$i] - $phi[$i - 1]) >= $epsilon);

    $phix = $phi[$i];

    $he = ($P/cos($phix)) - ($a/sqrt(1 - $e*$e * sin($phix)*sin($phix)));

    $coords['lambda']=$lambda;
    $coords['phi']=$phix;
    $coords['he']=$he;

    return $coords;
}

function ALG0013($Tx,$Ty,$Tz,$D,$Rx,$Ry,$Rz,$U)
{
    $V['X'] = $Tx + $U['X'] * (1 + $D) + $U['Z'] * $Ry - $U['Y'] * $Rz;
    $V['Y'] = $Ty + $U['Y'] * (1 + $D) + $U['X'] * $Rz - $U['Z'] * $Rx;
    $V['Z'] = $Tz + $U['Z'] * (1 + $D) + $U['Y'] * $Rx - $U['X'] * $Ry;

    return $V;
}

function ALG0019($lambda0,$phi0,$k0,$X0,$Y0,$a,$e)
{
    $lambdac = $lambda0;
    $n = sin($phi0);
    $C = $k0 * ALG0021($phi0,$a,$e) * tan (pi()/2 - $phi0) * exp ( $n * ALG0001($phi0,$e) );
    $Xs = $X0;
    $Ys = $Y0 + $k0 * ALG0021($phi0,$a,$e) * tan (pi()/2 - $phi0) ;

    $tab ['e'] = $e;
    $tab ['n'] = $n;
    $tab ['C'] = $C;
    $tab ['lambdac'] = $lambdac;
    $tab ['Xs'] = $Xs;
    $tab ['Ys'] = $Ys;

    return $tab;

}

function ALG0021($phi,$a,$e)
{
    $N = $a/sqrt( 1 - $e * $e * sin($phi) * sin($phi) );
    return $N;

}

function ALG0054($lambda0,$phi0,$X0,$Y0,$phi1,$phi2,$a,$e)
{
    $lambdac = $lambda0;
    $n = ( (log( (ALG0021($phi2,$a,$e)*cos($phi2))/(ALG0021($phi1,$a,$e)*cos($phi1)) ))/(ALG0001($phi1,$e) - ALG0001($phi2,$e) ));
    $C = ((ALG0021($phi1,$a,$e)* cos($phi1))/$n) * exp($n * ALG0001($phi1,$e));


    if ($phi0 == (pi()/2))
    {
        $Xs = $X0;
        $Ys = $Y0;
    }
    else
    {
            echo ('coucou');
        $Xs = $X0;
        $Ys = $Y0 + $C * exp(-1 * $n * ALG0001($phi0,$e));
    }

    $tab ['e'] = $e;
    $tab ['n'] = $n;
    $tab ['C'] = $C;
    $tab ['lambdac'] = $lambdac;
    $tab ['Xs'] = $Xs;
    $tab ['Ys'] = $Ys;

    return $tab;

}

function Lambert2WGS84($orig,$X,$Y)
{
    $epsilon = 0.00000000001;

    switch ($orig)
    {
        case 'LII' :
            $n = 0.7289686274;
            $c = 11745793.39;
            $Xs = 600000;
            $Ys = 6199695.768;
            $lambdac = 0.04079234433; // pour greenwich

            $e = 0.08248325676; //(première excentricité de l’ellipsoïde Clarke 1880 français)

            $he = 100;
            $a = 6378249.2;

            $Tx = -168;
            $Ty = -60;
            $Tz = +320;
            $D  = 0;
            $Rx = $Ry = $Rz = 0;
        break;
        case 'LIIe' :
            $n = 0.7289686274;
            $c = 11745793.39;
            $Xs = 600000;
            $Ys = 8199695.768;
            $lambdac = 0.04079234433; // pour greenwich

            $e = 0.08248325676; //(première excentricité de l’ellipsoïde Clarke 1880 français)

            $he = 100;
            $a = 6378249.2;

            $Tx = -168;
            $Ty = -60;
            $Tz = +320;
            $D  = 0;
            $Rx = $Ry = $Rz = 0;
        break;
        case 'L93' :
            $n = 0.7256077650;
            $c = 11745255.426;
            $Xs = 700000;
            $Ys = 12655612.050;
            $lambdac = 0.04079234433; // pour greenwich

            $e = 0.08248325676; //(première excentricité de l’ellipsoïde Clarke 1880 français)

            $he = 100;
            $a = 6378249.2;

            $Tx = -168;
            $Ty = -60;
            $Tz = +320;
            $D  = 0;
            $Rx = $Ry = $Rz = 0;
        break;
    }

    $coords = ALG0004($X,$Y,$n,$c,$Xs,$Ys,$lambdac,$e,$epsilon);

    $coords = ALG0009($coords['lambda'],$coords['phi'],$he,$a,$e);
    
    $coords = ALG0013($Tx,$Ty,$Tz,$D,$Rx,$Ry,$Rz,$coords);
    
    $a = 6378137.0;    // ellipsoïdes WGS84
    $f = 1/298.257223563;
    $b = $a*(1-$f);
    $e = sqrt(($a*$a - $b*$b)/($a*$a));

    $X = $coords['X'];
    $Y = $coords['Y'];
    $Z = $coords['Z'];
    $coords = ALG0012($X,$Y,$Z,$a,$e,$epsilon);
    
    $coords['lambda']=rad2deg($coords['lambda']);
    $coords['phi']   =rad2deg($coords['phi']);
    return $coords;
}
//
//http://www.ign.fr/telechargement/MPro/geodesie/CIRCE/transfo.pdf
//http://www.ign.fr/telechargement/MPro/geodesie/CIRCE/NTG_80.pdf
//http://www.ign.fr/telechargement/MPro/geodesie/CIRCE/NTG_71.pdf
//
print_r(Lambert2WGS84('LIIe',591647.56,2426659.65));
?>

[Sylvain33]

Comment on utilise le fichier php?
Dans un navigateur web, ca ne fonctionne pas !!!

[Pascal Boulerie]

Envoyé par Sylvain33

Comment utiliser un fichier php ?
Dans un navigateur web, ca ne fonctionne pas !!!

Si, en installant, sur son poste de travail un serveur AMP : Apache MySQL PHP. Sur un poste Linux, un LAMP (Linux AMP), et sur un poste Windows, un WAMP (Windows AMP).

Référence : http://fr.wikipedia.org/wiki/PHP
http://fr.wikipedia.org/wiki/PHP:_Hy...t_Preprocessor
"langage de programmation côté serveur"

[jm.2]

Un petit up pour ce superbe post qui m'a sauvé la mise.
Cela fait au moins un an que je reviens régulierement dessus, pour voire si l'implémentaiton de Bapt de Lambert2 etendu vers WGS84 était fonctionelle, et rien. Je ne sais pas si l'implémentation PHP est dèjà ok, car je n'ai rien pour la tester.

Aujourd'hui, j'ai trouvé les 3 erreurs dans l'implémentation de Bapt (Merci à toi et aux autres)
Les 2 premieres étaient dèjà connues:

• Ys pour Lambert2 étendu est 8199695.768 et non 6199695.768 qui est en lambert2 zone.
• e= 0.08181919106 (facteur 10)

Merci Pfpillon.
La dernière erreur était lors du passage WGS cartésien à WGS Géo (Algo 0012). La constante 'a' doit etre modifiée car ce n'est plus celle utilisée en Lambert, mais celle du WGS84:

• a = 6378137.0

Et maintenant ca marche!
J'ai aussi ajouté une conversion en degrés décimaux au lieu de donner le résultat en radian
Je donne mon implémentation ci-dessous (en Java)

import org.apache.log4j.Logger;

/**
 * 
 * @author Phiz
 *
 * http://www.forumsig.org/showthread.php?t=7418
 */
public class GeoTranslator 
{

    public static Logger logger = Logger.getLogger(GeoTranslator.class.getName());

    public static double[] Lambert2eToWGS84(double x, double y)
    {    
        /*
         * Quelques constantes ...
         */ 
        
        double n = 0.7289686274;            
        double c = 11745793.39;             // En mètres
        double Xs = 600000;                 // En mètres
        double Ys = 8199695.768;         // En mètres
        double l0 = 0.0;   //correspond à la longitude en radian de Paris (2°20'14.025" E) par rapport à Greenwich
        double e = 0.08248325676; //e du NTF (on le change après pour passer en WGS)
        double eps = Math.pow(10,-10);      // précision


        /***********************************************************
        *  coordonnées dans la projection de Lambert 2 à convertir *
        ************************************************************/
        double X = x;
        double Y = y;

        logger.debug("X = " + X);
        logger.debug("Y = " + Y);


        /*
         * Conversion Lambert 2 -> NTF géographique : ALG0004
         */

        logger.trace("\n-----------------------------------------------\n");
        logger.trace("Conversion Lambert 2 -> NTF géographique");

        double l;
        double L;
        double phi;
        double phi0;
        double phii;
        double phiprec;
        double R;
        double g;

        // Pour test : jeu de test d'ign. Cf. Algo 0004
        if(TEST)
        {    X = 1029705.083;
            Y = 272723.849;
            n = 0.760405966;
            l0 = 0.04079234433198;
            c = 11603796.9767;
            Ys = 5657616.674;
        } // fin test

        R = Math.sqrt(((X - Xs)*(X - Xs)) + ((Y - Ys)*(Y - Ys)));
        g = Math.atan((X - Xs) / (Ys - Y));

        l = l0 + (g / n);
        L = -(1 / n) * Math.log(Math.abs(R / c));

        
        phi0 = 2 * Math.atan(Math.exp(L)) - (Math.PI / 2.0);
        phiprec = phi0;
        phii= 2 * Math.atan((Math.pow(((1 + e * Math.sin(phiprec)) / (1 - e * Math.sin(phiprec))), e / 2.0) * Math.exp(L))) - (Math.PI / 2.0);

        while (!(Math.abs(phii - phiprec) < eps))
        {
            phiprec = phii;
            phii = 2 * Math.atan((Math.pow(((1 + e * Math.sin(phiprec)) / (1 - e * Math.sin(phiprec))), e / 2.0) * Math.exp(L))) - (Math.PI / 2.0);
        }

        phi = phii;

        logger.trace("Lambda = " + l + "rad = " + l * 200 /Math.PI + "gr");
        logger.trace("Phi    = " + phi + "rad = " + phi * 200 /Math.PI + "gr");



        /*
         * Conversion NTF géographique -> NTF cartésien : ALG0009
         */

        logger.trace("\n-----------------------------------------------\n");
        logger.trace("Conversion NTF géographique -> NTF cartésien");

        double N;
        double X_cart;
        double Y_cart;
        double Z_cart;

        double a = 6378249.2;
        double h = 100;         // En mètres

        double XWGS84;
        double YWGS84;
        double ZWGS84;


        // Pour test : jeu de test d'IGN. Cf. Algo 0009
        if(TEST)
        {    l = 0.01745329248;
            phi = 0.02036217457;
            h = 100;
        } // fin test

        N = a / (Math.pow((1 - (e * e) * (Math.sin(phi) * Math.sin(phi))), 0.5));
        X_cart = (N + h) * Math.cos(phi) * Math.cos(l);
        Y_cart = (N + h) * Math.cos(phi) * Math.sin(l);
        Z_cart = ((N * (1 - (e * e))) + h) * Math.sin(phi);


        logger.trace("X cartésien NTF = " + X_cart);
        logger.trace("Y cartésien NTF = " + Y_cart);
        logger.trace("Z cartésien NTF = " + Z_cart);



        /*
         * Conversion NTF cartésien -> WGS84 cartésien : ALG0013
         */


        logger.trace("\n-----------------------------------------------\n");
        logger.trace("Conversion NTF cartésien -> WGS84 cartésien");

        // Il s'agit d'une simple translation
        XWGS84 = X_cart - 168;
        YWGS84 = Y_cart - 60;
        ZWGS84 = Z_cart + 320;


        logger.trace("X cartésien WGS84 = " + XWGS84);
        logger.trace("Y cartésien WGS84 = " + YWGS84);
        logger.trace("Z cartésien WGS84 = " + ZWGS84);


        /*
         * Conversion WGS84 cartésien -> WGS84 géographique : ALG0012
         */

        logger.trace("\n-----------------------------------------------\n");
        logger.trace("Conversion WGS84 cartésien -> WGS84 géographique");

        double P;
        double phi840;
        double phi84prec;
        double phi84i;
        double phi84;
        double l840 = 0.04079234433;    // 0.04079234433 pour passer dans un référentiel par rapport au méridien
                                        // de Greenwich, sinon mettre 0
        double l84;
        e = 0.08181919106;              // On change e pour le mettre dans le système WGS84 au lieu de NTF
        a = 6378137.0;

        // Pour test : jeu de tests d'IGN. Cf. Algo 0012
        if(TEST)
        {    XWGS84 = 6376064.695;
            YWGS84 = 111294.623;
            ZWGS84 = 128984.725;
            l840 = 0;
            e = 0.08248325679;
        } // Fin tests


        P = Math.sqrt((XWGS84*XWGS84) + (YWGS84*YWGS84));

        l84 = l840 + Math.atan(YWGS84 / XWGS84);

        phi840 = Math.atan(ZWGS84 / (P * (1 - ((a * e * e))
                                    / Math.sqrt((XWGS84*XWGS84) + (YWGS84*YWGS84) + (ZWGS84*ZWGS84)))));
        
        phi84prec = phi840;

        //phi84i = Math.atan((ZWGS84 / P) * Math.pow(1 - ((a * e * e * Math.cos(phi84prec))
        //        / (P * Math.sqrt(1 - e * e * (Math.sin(phi84prec)*Math.sin(phi84prec))))), -1));
        phi84i = Math.atan((ZWGS84 / P) / (1 - ((a * e * e * Math.cos(phi84prec))
                / (P * Math.sqrt(1 - e * e * (Math.sin(phi84prec)*Math.sin(phi84prec)))))));

        while (!(Math.abs(phi84i - phi84prec) < eps))
        {
            phi84prec = phi84i;
            //phi84i = Math.atan((ZWGS84 / P) * Math.pow(1 - ((a * e * e * Math.cos(phi84prec))
            //            / (P * Math.sqrt(1 - ((e*e) * (Math.sin(phi84prec)*Math.sin(phi84prec)))))), -1));
            phi84i = Math.atan((ZWGS84 / P) / (1 - ((a * e * e * Math.cos(phi84prec))
                    / (P * Math.sqrt(1 - ((e*e) * (Math.sin(phi84prec)*Math.sin(phi84prec))))))));

        }

        phi84 = phi84i;

        logger.debug("lat WGS84  = " + l84 + "rad = " + l84 * 180.0 / Math.PI + " deg");
        logger.debug("long WGS84 = " + phi84 + "rad = " + phi84 * 180.0 / Math.PI + " deg");

        //return new double[]{phi84, l84};
        return new double[]{phi84 * 180.0 / Math.PI, l84 * 180.0 / Math.PI};
    }
}

PS: Les traces log4j peuvent etre effacées ou remplacées par du System.out.println()
PPS: Jai également changé les mises au carré:

pow(x, 2)

est optimisé par

(x*x)

Les resultats sont précis à 4 chiffres après la virgule lors d'une transformation de WGS84ToLambert2e puis Lambert2eToWGS84.
Encore merci pour ce forum.

[ataupin]

Bonjour,

j'ai essayé l'algorithme de Bapt pour la conversion du format Lambert 93 en WSG84.
et, visiblement j'ai des soucis de calcul car il y a des décalages.

Est-ce que quelqu'un qui a déjà l'algorithme de conversion Lambert 93 en Wsg84 ?

D'avance merci.

[cuq]

Merci .... beaucoup

Ca m'a permis de résoudre mon problème de conversion de Lambert vers WGS84 .

Si ca intéresse du monde j'ai fais le code en VB6 (marche en VBA) et l'algo marche apparement avec tous les lambert I II

Suffit de changer les valeurs de

Const n = 0.7289686274
Const c = 11745793.39
Const Xs = 600000
Const Ys = 8199695.768

J'ai pris les valeurs dans le document de l'ign NTG_71.Pdf ( Projection cartographique conique conforme de Lambert ) en dernière page il y a les valeurs constantes lambert France

Testé avec du Lambert I et Lambert II ca marche le reste je sais pas.

Merci a tous

[cuq]

bon J'ai utilisé le code de Jm.2 pour passer des coordonnées Lambert vers WGS84

J'ai repris les constantes suivantes données par L'iGN

' |---------------------------------------------------------------------------------------------------------------|
' | Const | 1 'Lambert I | 2 'Lambert II | 3 'Lambert III | 4 'Lambert IV | 5 'Lambert II Etendue | 6 'Lambert 93 |
' |-------|--------------|---------------|----------------|---------------|-----------------------|---------------|
' |    n  | 0.7604059656 |  0.7289686274 |   0.6959127966 | 0.6712679322  |    0.7289686274       |  0.7256077650 |
' |-------|--------------|---------------|----------------|---------------|-----------------------|---------------|
' |    c  | 11603796.98  |  11745793.39  |   11947992.52  | 12136281.99   |    11745793.39        |  11754255.426 |
' |-------|--------------|---------------|----------------|---------------|-----------------------|---------------|
' |    Xs |   600000.0   |    600000.0   |   600000.0     |      234.358  |    600000.0           |     700000.0  |
' |-------|--------------|---------------|----------------|---------------|-----------------------|---------------|
' |    Ys | 5657616.674  |  6199695.768  |   6791905.085  |  7239161.542  |    8199695.768        | 12655612.050  |
' |---------------------------------------------------------------------------------------------------------------|

Et le calcul marche très bien pour Lambert I , II , III , IV et II etendues mais je n'ai pas les bonnes valeurs pour Lambert 93.

Il y a t'il une différence en plus pour le calcul en lambert 93 en plus des constantes de départ.

Par avance merci

Laurent
_
[EDIT]
_
Bon je m'approche du résultat mais j'ai encore une imprécision.

Donc Lambert -93 la longitude d'origine c'est 3°E Greenwich est plus Paris comme les autres donc :

'Pour lambert 93 Longitude origine 3° Est Greenwich
l0 = (3 / 180 * Pi) - 0.04079234433198

Et je pense qu'il n'est pas nécessaire de faire la Conversion NTF cartésien - WGS84 cartésien

Vu que lambert 93 est associé au système géodésique RGF93 alors que les autres sont associé au systèmes NTF

Dans ce cas je trouve un résultat qui est "bon" à 8 ou 10 mètres près ..

Si j'ai commis une erreur merci de me le dire.

A+

LAurent

[Bude]

Ta formule de conversion m'interresse beaucoup peux-tu la publier et l'expliquer ?

merci

daniel

Envoyé par cuq

Merci .... beaucoup

Ca m'a permis de résoudre mon problème de conversion de Lambert vers WGS84 .

Si ca intéresse du monde j'ai fais le code en VB6 (marche en VBA) et l'algo marche apparement avec tous les lambert I II

Suffit de changer les valeurs de

Const n = 0.7289686274
Const c = 11745793.39
Const Xs = 600000
Const Ys = 8199695.768

J'ai pris les valeurs dans le document de l'ign NTG_71.Pdf ( Projection cartographique conique conforme de Lambert ) en dernière page il y a les valeurs constantes lambert France

Testé avec du Lambert I et Lambert II ca marche le reste je sais pas.

Merci a tous

[cuq]

Voilà le code en VB6

Option Explicit

Global System_Coord As Integer

'---------------------------------------------------------------------------
' Point en 2D
'---------------------------------------------------------------------------
Public Type Point2
    X As Double
    Y As Double
End Type


' Convertion de données en lambert II vers position GPS
' X et Y en mètres
' Option Système de coordonnées
' |---------------------------------------------------------------------------------------------------------------|
' | Const | 1 'Lambert I | 2 'Lambert II | 3 'Lambert III | 4 'Lambert IV | 5 'Lambert II Etendue | 6 'Lambert 93 |
' |-------|--------------|---------------|----------------|---------------|-----------------------|---------------|
' |    n  | 0.7604059656 |  0.7289686274 |   0.6959127966 | 0.6712679322  |    0.7289686274       |  0.7256077650 |
' |-------|--------------|---------------|----------------|---------------|-----------------------|---------------|
' |    c  | 11603796.98  |  11745793.39  |   11947992.52  | 12136281.99   |    11745793.39        |  11754255.426 |
' |-------|--------------|---------------|----------------|---------------|-----------------------|---------------|
' |    Xs |   600000.0   |    600000.0   |   600000.0     |      234.358  |    600000.0           |     700000.0  |
' |-------|--------------|---------------|----------------|---------------|-----------------------|---------------|
' |    Ys | 5657616.674  |  6199695.768  |   6791905.085  |  7239161.542  |    8199695.768        | 12655612.050  |
' |---------------------------------------------------------------------------------------------------------------|
'
Function LTversWGS84(ByVal X As Double, ByVal Y As Double, Optional SystemCoord As Integer = 5) As Point2

Dim longitude As Double
Dim latitude As Double
Dim Longi As Double
Dim L As Double
Dim phi As Double
Dim phi0 As Double
Dim phii As Double
Dim phiprec As Double
Dim R As Double
Dim g As Double
Dim VarN As Double
Dim X_cart As Double
Dim Y_cart As Double
Dim Z_cart As Double
Dim XWGS84 As Double
Dim YWGS84 As Double
Dim ZWGS84 As Double
Dim p As Double
Dim phi840 As Double
Dim phi84prec As Double
Dim phi84i As Double
Dim phi84 As Double
Dim l84 As Double
Dim l840 As Double
Dim Const_a As Double
Dim Const_e As Double

Dim Decalage_X_cart As Double
Dim Decalage_Y_cart As Double
Dim Decalage_Z_cart As Double

Dim l0 As Double

'Debug.Print "Conversion Lambert 2 - NTF géographique"

'--  Quelques constantes ...  -->
Dim n As Double
Dim C As Double           '--   En mètres   --
Dim Xs As Double          '--   En mètres   --
Dim Ys As Double          '--   En mètres   --
    
    

Const eps = 0.000000000001        '--   précision   --
Const h = 100                   '--   En mètres   --

l0 = 0                          '--  correspond à la longitude en radian de Paris (2°20'14.025" E) par rapport à Greenwich   --
l840 = 2.337229167 / 180 * Pi '0.04079234433            '--   0.04079234433 pour passer dans un référentiel par rapport au méridien   --
                                '--   de Greenwich, sinon mettre 0   --

Const_a = 6378249.2
Const_e = 0.08248325676         '--  Const_e du NTF (on le change après pour passer en WGS)   --
                                '-- (première excentricité de l’ellipsoïde Clarke 1880 français). ---
                                
 Decalage_X_cart = -168         '--   En mètres   --
 Decalage_Y_cart = -60          '--   En mètres   --
 Decalage_Z_cart = 320          '--   En mètres   --


'Definition selon systèmes de coordonnées Lambert
Select Case SystemCoord

Case 1 'Lambert I

     n = 0.7604059656
     C = 11603796.98
     Xs = 600000
     Ys = 5657616.674
     
Case 2 'Lambert II

     n = 0.7289686274
     C = 11745793.39
     Xs = 600000
     Ys = 6199695.768

Case 3 'Lambert III
     n = 0.6959127966
     C = 11947992.52
     Xs = 600000
     Ys = 6791905.085
     
Case 4 'Lambert IV
     n = 0.6712679322
     C = 12136281.99
     Xs = 234.358
     Ys = 7239161.542
     
Case 5 'Lambert II Etendue
     n = 0.7289686274
     C = 11745793.39
     Xs = 600000
     Ys = 8199695.768
     
Case 6 'Lambert 93
     
     
        n = 0.725607765
        C = 11754255.426
        Xs = 700000
        Ys = 12655612.05


        'Pour lambert 93 Longitude origine 3° Est Greenwich
        l0 = (3 / 180 * Pi) ' - 0.04079234433198
        l840 = 0
        
        
        Decalage_X_cart = 0
        Decalage_Y_cart = 0
        Decalage_Z_cart = 0

End Select

R = Sqr(((X - Xs) * (X - Xs)) + ((Y - Ys) * (Y - Ys)))
g = Atn((X - Xs) / (Ys - Y))

Longi = l0 + (g / n)
L = -(1 / n) * Log(Abs(R / C))


phi0 = 2 * Atn(Exp(L)) - (Pi / 2#)
phiprec = phi0
phii = 2 * Atn(((((1 + Const_e * Sin(phiprec)) / (1 - Const_e * Sin(phiprec))) ^ (Const_e / 2#)) * Exp(L))) - (Pi / 2#)

While Not (Abs(phii - phiprec) < eps)
        phiprec = phii
        phii = 2 * Atn((((((1 + Const_e * Sin(phiprec)) / (1 - Const_e * Sin(phiprec))) ^ (Const_e / 2#)) * Exp(L)))) - (Pi / 2#)

Wend


phi = phii

'Debug.Print "Lambda = " & Longi & " rad = " & Longi * 200 / Pi & "gr"
'Debug.Print "Phi    = " & phi & " rad = " & phi * 200 / Pi & "gr"



'-- Conversion NTF géographique - NTF cartésien : ALG0009   --
'Debug.Print "Conversion NTF géographique - NTF cartésien"


VarN = Const_a / ((1 - (Const_e * Const_e) * (Sin(phi) * Sin(phi))) ^ 0.5)
X_cart = (VarN + h) * Cos(phi) * Cos(Longi)
Y_cart = (VarN + h) * Cos(phi) * Sin(Longi)
Z_cart = ((VarN * (1 - (Const_e * Const_e))) + h) * Sin(phi)

'Debug.Print "X cartésien NTF = " & X_cart
'Debug.Print "Y cartésien NTF = " & Y_cart
'Debug.Print "Z cartésien NTF = " & Z_cart

'-- Conversion NTF cartésien - WGS84 cartésien : ALG0013   --
'Debug.Print "Conversion NTF cartésien - WGS84 cartésien"

'--   Il s'agit d'une simple translation   --
XWGS84 = X_cart + Decalage_X_cart
YWGS84 = Y_cart + Decalage_Y_cart
ZWGS84 = Z_cart + Decalage_Z_cart


'Debug.Print "X cartésien WGS84 = " & XWGS84
'Debug.Print "Y cartésien WGS84 = " & YWGS84
'Debug.Print "Z cartésien WGS84 = " & ZWGS84


'-- Conversion WGS84 cartésien - WGS84 géographique : ALG0012   --
'Debug.Print "Conversion WGS84 cartésien - WGS84 géographique"

Const_e = 0.08181919106              '--   On change Const_e pour le mettre dans le système WGS84 au lieu de NTF   --
Const_a = 6378137


p = Sqr((XWGS84 * XWGS84) + (YWGS84 * YWGS84))

l84 = l840 + Atn(YWGS84 / XWGS84)

phi840 = Atn(ZWGS84 / (p * (1 - ((Const_a * Const_e * Const_e)) / Sqr((XWGS84 * XWGS84) + (YWGS84 * YWGS84) + (ZWGS84 * ZWGS84)))))

phi84prec = phi840


phi84i = Atn((ZWGS84 / p) / (1 - ((Const_a * Const_e * Const_e * Cos(phi84prec)) / (p * Sqr(1 - Const_e * Const_e * (Sin(phi84prec) * Sin(phi84prec)))))))

While Not (Abs(phi84i - phi84prec) < eps)
        phi84prec = phi84i
        phi84i = Atn((ZWGS84 / p) / (1 - ((Const_a * Const_e * Const_e * Cos(phi84prec)) / (p * Sqr(1 - ((Const_e * Const_e) * (Sin(phi84prec) * Sin(phi84prec))))))))
Wend

phi84 = phi84i

LTversWGS84.X = l84 * 180 / Pi
LTversWGS84.Y = phi84 * 180 / Pi

'Debug.Print "latitude WGS84  = " & l84 & " rad = " & l84 * 180 / Pi & " deg"
'Debug.Print "longitude WGS84 = " & phi84 & " rad = " & phi84 * 180 / Pi & " deg"
End Function

'conversion longitude latitude en texte
Function Texte_Position(longitude As Double, latitude As Double) As String
Dim TexteTempo As String
Dim Temp As Double
Dim degres As Integer
Dim minutes As Integer
Dim secondes As Integer
Dim Signe As String

             Temp = Abs(longitude)
             degres = Int(Temp)
             minutes = Int((Temp - degres) * 60)
             secondes = Round((Temp - degres - minutes / 60) * 3600)
                
            If (Temp > 0) Then
                Signe = " E "
            Else
                Signe = " O "
            End If
                
             TexteTempo = degres & "° " & minutes & "' " & secondes & "'' " & Signe
                       
            Temp = Abs(latitude)
            degres = Int(Temp)
            minutes = Int((Temp - degres) * 60)
            secondes = Round((Temp - degres - minutes / 60) * 3600)
            
            If (Temp > 0) Then
                Signe = " N "
            Else
                Signe = " S "
            End If
                
            TexteTempo = TexteTempo & vbCrLf & degres & "° " & minutes & "' " & secondes & "'' " & Signe
            
            Texte_Position = TexteTempo

End Function

La méthode est la même que celle des postes précédent en Java

[Bude]

Merci , le code fonctionne pratiquement à l' identique en vba.
a+
daniel

[Dam81]

Bonjour à tous,

Pour un projet informatique, j'ai utilisé le code de Jm.2 (conversion Lambert II Etendu --> WGS84).

Cependant, au lieu d'utiliser Xs = 600000 et Ys = 8199695.768, j'ai utilisé Xs = 599993.15 et Ys = 8199694.3 déterminés de manière empirique afin de m'approcher des résultats de Circé (l'outil de l'IGN).

J'ai testé la transformation sur environ 280000 points situés un peu partout en France métropolitaine et en Corse et j'ai systématiquement obtenu des résultats plus proches de ceux de Circé.

Je souhaiterais savoir le danger (au niveau des résultats de conversion) qui peut résulter de l'utilisation de Xs et Ys biaisés ainsi que les problèmes que cela peut engendrer.

Merci par avance.

[Jérôme C]

Dam, quand tu compares avec Circée, tu utilises la transformation avec paramètres ou la transformation grille ?

Le problème dans cette transformation n'est pas dans le passage de X, Y --> lambda, phi, mais plutôt le passage de (lambda, phi)NTF --> (lambda, phi)WGS84. Jm.2 utilise des paramètres globaux alors que l'IGN conseille l'utilisation de la grille GRD97A, où pour chaque région d'une "seconde²" il y a des paramètres de corrections différents, cf http://lambert93.ign.fr/index.php?id=30#c71 La carte des déformations http://georezo.net/forum/viewtopic.php?pid=90590#p90590 aide à voir la différence et le danger d'utiliser les paramètres globaux.

[Dam81]

Bonjour Jérôme C et merci pour ta réponse.

J'utilise Circé v3.2 en mode "Transformation standard" pour réaliser une conversion depuis le système "NTF (Paris) / Planes / Lambert 2 étendu" vers le système "WGS84 / Géographiques".
D'après la description faite dans Circé :
"La transformation standard est effectuée entre les systèmes WGS84, NTF et ED50 avec des paramètres standard (translation, rotation et facteur d'échelle) uniques pour toute la France."
Cette transformation ne semble donc pas utiliser de grille.
Pourtant, il y a un écart entre les résultats issus du script de jm.2 et ceux issus de Circé en "Transformation standard" (J'ai besoin de 6 décimales).

Cependant, si j'ai bien compris ta remarque, tu m'indiques que je modifie les paramètres de la transformation "Lambert II Etendu --> NTF Géographique" alors que les variations sont issues des étapes "NTF Géographique --> NTF Cartésien --> WGS84 Cartésien --> WGS84 Géographique".
Mon raisonnement totalement empirique est donc d'un point de vue théorique parfaitement faux.
Je vais voir s'il n'est pas mieux de m'accomoder de ces petites variations entre le script de Jm.2 et Circé plutôt que de chercher à les corriger à tout prix.

[Jérôme C]

(J'ai besoin de 6 décimales).

Est-ce que tu es en train de déclarer que tu as des différences de l'ordre du micromètre que tu essaies de corriger ? Même s'il y a une erreur, qui s'en soucit ? Tu sais ce que c'est qu'un micromètre ?

Pour essayer de savoir ce que je voulais t'indiquer, tu as tout à fait raison.