[Données] Formule de conversion coordonnées

[Pierre-Alain]

Envoyé par Bapt

Par le plus grand des hasards, l'un d'entre vous aurait-il implémenté la fonction inverse ?

As tu essayé de réaliser cette conversion en suivant les étapes inverses ?

[Amelnast]

Bonjour Phiz,

tu indiques
//puis on calcule la latitude isométrique
L = ln(tan(PI/4 + phi_n/2) * pow(((1-e_n*sin(phi_n))/(1+e_n*sin(phi_n))),(e_n/2)))
moi je viens de récupérer un programme de conversion de coordonnées lat/long en degrées où on indique
L = 0.5*log((1+sin(Phi))/(1-sin(Phi)))-exc*0.5*log((1+exc*sin(Phi))/(1-exc*sin(Phi)))

D'où vient cette différence ? J'ai du mal à suivre entre toutes les formules

Merci pour ta réponse

[csegonds]

Bon, merci les gars c'est du super boulot.
Mais j'effectue un portage en WinDev de la classe java et je dois avouer que j'ai quelques soucis.
Je n'obtiens pas les résultats escomptés.

Je pense que cela doit venir des fonctions utilisées.

J'avais des formules tirées d'une feuille excel venant de l'IGN, mais malheuresement ces formules ne sont pas super précises contrairement à votre code qui donne quasiment les même résultat que le logiciel CIRCE de l'IGN (je considére que 2 chiffres, voire 3 après la virgule suffit amplement).

Est ce que quelqu'un pourrait me donner les différentes valeurs des calculs intermédiaires pour que je puisse mettre au point mon programme (avec un jeu d'essai).
Merci.
_
[EDIT]
_
Bon, ca y est j'ai trouvé....
WinDev travaille en degree au niveau des fonctions SIN, COS, ATAN, TAN.
Donc ca marche pas....
J'ai implémenté en C et tout fonctionne.

Merci à tous.

Effectivement, je suis interessé par la formule inverse :
Convertir du lambert 2 en WGS84......

A+

[Bapt]

Envoyé par Pierre-Alain

As tu essayé de réaliser cette conversion en suivant les étapes inverses ?

Oui j'ai essayé et je pense que je suis en bonne voie puisque j'ai implémenté les 4 algorithmes inverses. Ceux-ci valident les jeux de test fournis par IGN donc je pense que l'implémentation est correcte. En revanche, je n'arrive pas à obtenir un résultat cohérent en enchaînant les 4. Je pense qu'il s'agit d'un problème d'unités ou de constantes que je choisis mal. C'est pour cela que j'aurais souhaité comparer mon travail avec celui de quelqu'un d'autre.

Ci-dessous, mon implémentation en C#

using System;
using System.Collections.Generic;
using System.Text;

namespace LT2versWGS84
{
    class Program
    {
        static void Main(string[] args)
        {

            /*
             * Quelques constantes ...
             */ 

            double n = 0.7289686274;            
            double c = 11745793.39;             // En mètres

            double Xs = 600000;                 // En mètres
            double Ys = 6199695.768;            // En mètres

            double l0 = 0;// 0.04079234433;     // 0 pour Paris, l'autre pour Greenwich, en radians

            double e = 0.08248325676;
            double eps = Math.Pow(10,-11);      // précision


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


            Console.WriteLine("X = " + X);
            Console.WriteLine("Y = " + Y);


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

            Console.WriteLine("\n-----------------------------------------------\n");
            Console.WriteLine("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
            /*X = 1029705.083;
            Y = 272723.849;
            n = 0.760405966;
            c = 11603796.9767;
            Ys = 5657616.674;*/
            // fin test

            R = Math.Sqrt(Math.Pow((X - Xs),2) + Math.Pow((Y - Ys),2));
            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);
            phiprec = phi0;
            phii= 2 * Math.Atan((Math.Pow(((1 + e * Math.Sin(phiprec)) / (1 - e * Math.Sin(phiprec))), e / 2)
                     * Math.Exp(L))) - (Math.PI / 2);

            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)
                        * Math.Exp(L))) - (Math.PI / 2);
            }

            phi = phii;

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



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

            Console.WriteLine("\n-----------------------------------------------\n");
            Console.WriteLine("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
            /*l = 0.01745329248;
            phi = 0.02036217457;
            h = 100;*/
            // fin test

            N = a / (Math.Pow((1 - Math.Pow(e, 2) * Math.Pow(Math.Sin(phi), 2)), 1 / 2));
            X_cart = (N + h) * Math.Cos(phi) * Math.Cos(l);
            Y_cart = (N + h) * Math.Cos(phi) * Math.Sin(l);
            Z_cart = ((N * (1 - Math.Pow(e,2))) + h) * Math.Sin(phi);


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



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


            Console.WriteLine("\n-----------------------------------------------\n");
            Console.WriteLine("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;


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


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

            Console.WriteLine("\n-----------------------------------------------\n");
            Console.WriteLine("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.8181919106;               // On change e pour le mettre dans le système WGS84 au lieu de NTF

            // Pour test : jeu de tests d'IGN. Cf. Algo 0012
            /*XWGS84 = 6376064.695;
            YWGS84 = 111294.623;
            ZWGS84 = 128984.725;
            l840 = 0;
            e = 0.08248325679;*/
            // Fin tests


            P = Math.Sqrt(Math.Pow(XWGS84, 2) + Math.Pow(YWGS84, 2));

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

            phi840 = Math.Atan(ZWGS84 / (P * (1 - ((a * Math.Pow(e, 2))
                                        / Math.Sqrt(Math.Pow(XWGS84, 2) + Math.Pow(YWGS84, 2) + Math.Pow(ZWGS84, 2))))));
            
            phi84prec = phi840;

            phi84i = Math.Atan((ZWGS84 / P) * Math.Pow(1 - ((a * Math.Pow(e, 2) * Math.Cos(phi84prec))
                            / (P * Math.Sqrt(1 - Math.Pow(e, 2) * Math.Pow(Math.Sin(phi84prec), 2)))), -1));

            while (!(Math.Abs(phi84i - phi84prec) < eps))
            {
                phi84prec = phi84i;
                phi84i = Math.Atan((ZWGS84 / P) * Math.Pow(1 - ((a * Math.Pow(e, 2) * Math.Cos(phi84prec))
                            / (P * Math.Sqrt(1 - (Math.Pow(e, 2) * Math.Pow(Math.Sin(phi84prec), 2))))), -1));

            }

            phi84 = phi84i;

            Console.WriteLine("lat WGS84  = " + l84 + "rad = " + l84 * 200 / Math.PI + "gr");
            Console.WriteLine("long WGS84 = " + phi84 + "rad = " + phi84 * 200 / Math.PI + "gr");



            Console.ReadLine();
        }
    }
}

[Pierre-Alain]

Envoyé par Bapt

Oui j'ai essayé et je pense que je suis en bonne voie puisque j'ai implémenté les 4 algorithmes inverses. Ceux-ci valident les jeux de test fournis par IGN donc je pense que l'implémentation est correcte. En revanche, je n'arrive pas à obtenir un résultat cohérent en enchaînant les 4. Je pense qu'il s'agit d'un problème d'unités ou de constantes que je choisis mal.

Bonsoir,

Au niveau logique, je ne te suis pas : Tu dis que chaque algorithme marche séparement(validé par des tests avec des données IGN), mais quand tu les regroupes, cela ne fonctionne plus ?En traçant chaque variable, tu n'y arrives pas ?Tu essayé de séparer tes quatres étapes en fonction ?

ca parait bizarre

Envoyé par Bapt

C'est pour cela que j'aurais souhaité comparer mon travail avec celui de quelqu'un d'autre.

Ci-dessous, mon implémentation en C#

J'peux pas trop t'aider, j'connais le C, C++ mais le C# ...

[pfpillon]

Bonjour,

Envoyé par Bapt

            double Ys = 6199695.768;            // En mètres

déjà, il me semble que tu n'utilises pas le 'bon' Ys (celui de Lambert II +) qui est 8199695.768.

en utilisant cette valeur, on trouve la bonne latitude.

Toutefois, la longitude reste erronée.

Si quelqu'un a une idée...
_
[EDIT]
_

Envoyé par Bapt

            e = 0.8181919106;               // On change e pour le mettre dans le système WGS84 au lieu de NTF

            // Pour test : jeu de tests d'IGN. Cf. Algo 0012
            /*XWGS84 = 6376064.695;
            YWGS84 = 111294.623;
            ZWGS84 = 128984.725;
            l840 = 0;
            e = 0.08248325679;*/
            // Fin tests

En utilisant comme e la valeur des jeux de tests de l'IGN (0.08248325679), je trouve la bonne longitude et la bonne latitude.

[EDIT]

FAUX --> on a quelques centaines de mètres d'erreur en latitude

[pfpillon]

En fait, le e (e = 0.8181919106) était presque bon mais erroné d'un facteur 10 -->
e = 0.08181919106;

Avec cette valeur, on a une réponse quasi parfaite (quelques mètres d'erreur)

[Phiz]

Envoyé par Amelnast

Bonjour Phiz,

tu indiques
//puis on calcule la latitude isométrique
L = ln(tan(PI/4 + phi_n/2) * pow(((1-e_n*sin(phi_n))/(1+e_n*sin(phi_n))),(e_n/2)))
moi je viens de récupérer un programme de conversion de coordonnées lat/long en degrées où on indique
L = 0.5*log((1+sin(Phi))/(1-sin(Phi)))-exc*0.5*log((1+exc*sin(Phi))/(1-exc*sin(Phi)))

D'où vient cette différence ? J'ai du mal à suivre entre toutes les formules

Merci pour ta réponse

Tout ce que je peux te dire c'est que la formule que j'utilise est extraite directement d'un document trouvé sur le site de l'IGN
( http://www.ign.fr/telechargement/MPr...RCE/NTG_71.pdf )

Désolé, je ne vois pas le lien entre nos deux formules et pourtant j'ai retourné la mienne dans tout les sens pour faire apparaitre la tienne et je n'ai retrouvé que la moitié de la formule (-exc*0.5*......)

[Aurelien59]

Bonjour tout le monde,

J'ai moi même aussi repris les algos de l'IGN et ai les mêmes problèmes que Bapt. Les algos sont bons avec les jeux d'essai de l'IGN, mais dès que l'on test avec ses propres donnée, les résultats sont... "amusants"...

Voici mon code:

using System;
using System.Drawing;
using System.Collections;
using System.ComponentModel;
using System.Windows.Forms;
using System.Data;

namespace ConversionGEO
{
/// <summary>
/// Description résumée de Form1.
/// </summary>
public class Form1 : System.Windows.Forms.Form
{
private System.Windows.Forms.TextBox textBox1;
private System.Windows.Forms.TextBox textBox2;
private System.Windows.Forms.Button button1;
private System.Windows.Forms.Label label1;
private System.Windows.Forms.Label label2;
private System.Windows.Forms.Label label3;
private System.Windows.Forms.Label label4;
private System.Windows.Forms.ComboBox comboBox1;
/// <summary>
/// Variable nécessaire au concepteur.
/// </summary>
private System.ComponentModel.Container components = null;
private double latitude, longitude, n, c, Xs, Ys, e, lambda0, a, epsilon, nw, cw, Xw, Yw, ew, lambda0w, aw;
private System.Windows.Forms.Label label5;
private System.Windows.Forms.Label label6;
private int typeConversion=3;

public Form1()
{
//
// Requis pour la prise en charge du Concepteur Windows Forms
//
InitializeComponent();

//
// TODO : ajoutez le code du constructeur après l'appel à InitializeComponent
//
}

/// <summary>
/// Nettoyage des ressources utilisées.
/// </summary>
protected override void Dispose( bool disposing )
{
if( disposing )
{
if (components != null)
{
components.Dispose();
}
}
base.Dispose( disposing );
}

#region Code généré par le Concepteur Windows Form
/// <summary>
/// Méthode requise pour la prise en charge du concepteur - ne modifiez pas
/// le contenu de cette méthode avec l'éditeur de code.
/// </summary>
private void InitializeComponent()
{
this.textBox1 = new System.Windows.Forms.TextBox();
this.textBox2 = new System.Windows.Forms.TextBox();
this.button1 = new System.Windows.Forms.Button();
this.label1 = new System.Windows.Forms.Label();
this.label2 = new System.Windows.Forms.Label();
this.label3 = new System.Windows.Forms.Label();
this.label4 = new System.Windows.Forms.Label();
this.comboBox1 = new System.Windows.Forms.ComboBox();
this.label5 = new System.Windows.Forms.Label();
this.label6 = new System.Windows.Forms.Label();
this.SuspendLayout();
//
// textBox1
//
this.textBox1.Location = new System.Drawing.Point(128, 22);
this.textBox1.Name = "textBox1";
this.textBox1.Size = new System.Drawing.Size(160, 20);
this.textBox1.TabIndex = 0;
this.textBox1.Text = "";
//
// textBox2
//
this.textBox2.Location = new System.Drawing.Point(128, 46);
this.textBox2.Name = "textBox2";
this.textBox2.Size = new System.Drawing.Size(160, 20);
this.textBox2.TabIndex = 1;
this.textBox2.Text = "";
//
// button1
//
this.button1.Location = new System.Drawing.Point(280, 80);
this.button1.Name = "button1";
this.button1.Size = new System.Drawing.Size(96, 24);
this.button1.TabIndex = 2;
this.button1.Text = "Convertir";
this.button1.Click += new System.EventHandler(this.button1_Click);
//
// label1
//
this.label1.Location = new System.Drawing.Point(8, 24);
this.label1.Name = "label1";
this.label1.Size = new System.Drawing.Size(112, 16);
this.label1.TabIndex = 3;
this.label1.Text = "latitude à convertir";
//
// label2
//
this.label2.Location = new System.Drawing.Point(8, 48);
this.label2.Name = "label2";
this.label2.Size = new System.Drawing.Size(112, 16);
this.label2.TabIndex = 4;
this.label2.Text = "longitude à convertir";
//
// label3
//
this.label3.Location = new System.Drawing.Point(304, 24);
this.label3.Name = "label3";
this.label3.Size = new System.Drawing.Size(184, 16);
this.label3.TabIndex = 5;
//
// label4
//
this.label4.Location = new System.Drawing.Point(304, 48);
this.label4.Name = "label4";
this.label4.Size = new System.Drawing.Size(184, 16);
this.label4.TabIndex = 6;
//
// comboBox1
//
this.comboBox1.Items.AddRange(new object[] {
"WGS84 -> Lambert I",
"WGS84 -> Lambert II",
"WGS84 -> Lambert II étendu",
"WGS84 -> Lambert III",
"WGS84 -> Lambert IV",
"WGS84 -> Lambert -93",
"Lambert I -> WGS84 ",
"Lambert II -> WGS84",
"Lambert II étendu -> WGS84",
"Lambert III -> WGS84",
"Lambert IV -> WGS84",
"Lambert -93 -> WGS84"});
this.comboBox1.Location = new System.Drawing.Point(80, 80);
this.comboBox1.Name = "comboBox1";
this.comboBox1.Size = new System.Drawing.Size(176, 21);
this.comboBox1.TabIndex = 9;
this.comboBox1.Text = "choisissez le sens de conversion";
//
// label5
//
this.label5.Location = new System.Drawing.Point(152, 168);
this.label5.Name = "label5";
this.label5.Size = new System.Drawing.Size(184, 16);
this.label5.TabIndex = 10;
//
// label6
//
this.label6.Location = new System.Drawing.Point(144, 200);
this.label6.Name = "label6";
this.label6.Size = new System.Drawing.Size(184, 16);
this.label6.TabIndex = 11;
//
// Form1
//
this.AutoScaleBaseSize = new System.Drawing.Size(5, 13);
this.ClientSize = new System.Drawing.Size(496, 245);
this.Controls.Add(this.label6);
this.Controls.Add(this.label5);
this.Controls.Add(this.comboBox1);
this.Controls.Add(this.label4);
this.Controls.Add(this.label3);
this.Controls.Add(this.label2);
this.Controls.Add(this.label1);
this.Controls.Add(this.button1);
this.Controls.Add(this.textBox2);
this.Controls.Add(this.textBox1);
this.Name = "Form1";
this.Text = "Form1";
this.ResumeLayout(false);

}
#endregion

/// <summary>
/// Point d'entrée principal de l'application.
/// </summary>
[STAThread]
static void Main()
{
Application.Run(new Form1());
}

/**
* phi : latitude.
* a : demi-grand axe de l’ellipsoïde.
* e : première excentricité de l’ellipsoïde.
*
* renvoi :
* N : la grande normale
* */
private double calculDeLaGrandeNormale(double phi, double a, double e)
{
return a / Math.Sqrt(1 - e * e * Math.Sin(phi) * Math.Sin(phi));
}

/**
* lambda : longitude par rapport au méridien origine.
* phi : latitude.
* he : hauteur au-dessus de l’ellipsoïde.
* a : demi-grand axe de l’ellipsoïde.
* e : première excentricité de l’ellipsoïde.
*
* Renvoi:
* resultat[0] -> X
* resultat[1] -> Y
* resultat[2] -> Z
* */
private double[] conversionGeographiqueVersCartesienne(double lambda, double phi, double he, double a, double e)
{
double[] resultat = new double[3];
double N = calculDeLaGrandeNormale(phi,a,e);
resultat[0] = (N + he) * Math.Cos(phi) * Math.Cos(lambda);
resultat[1] = (N + he) * Math.Cos(phi) * Math.Sin(lambda);
resultat[2] = (N * (1 - e*e) + he) * Math.Sin(phi);
return resultat;
}

/**
* X, Y, Z : coordonnées cartésiennes.
* a : demi-grand axe de l’ellipsoïde.
* e : première excentricité de l’ellipsoïde.
* epsilon : tolérance de convergence.
*
* renvoi:
* resultat[0] -> lambda : longitude par rapport au méridien origine.
* resultat[1] -> phi : latitude.
* resultat[2] -> he : hauteur au-dessus de l’ellipsoïde.
* */
private double[] conversionCartesienneVersGeorgraphique(double X, double Y, double Z, double a, double e, double epsilon)
{
double[] resultat = new double[3];
resultat[0] = Math.Atan(Y / X);
double temp = Math.Atan(Z / (Math.Sqrt(X * X + Y * Y) * (1 - (a * e * e) / Math.Sqrt(X * X + Y * Y + Z * Z))));
double tempo = temp + 1;
while(Math.Abs(temp - tempo) >= epsilon)
{
tempo = temp;
temp = Math.Atan(Z / (Math.Sqrt(X * X + Y * Y)) * Math.Pow(1 - ((a * e * e * Math.Cos(tempo)) / (Math.Sqrt(X * X + Y * Y) * Math.Sqrt(1 - e * e * Math.Sin(tempo)* Math.Sin(tempo)))),-1));
}
resultat[1] = temp;
resultat[2] = (Math.Sqrt(X * X + Y * Y) / Math.Cos(temp)) - (a / Math.Sqrt(1 - e * e * Math.Sin(temp)* Math.Sin(temp)));
return resultat;
}

/**
* phi : latitude.
* e : première excentricité de l’ellipsoïde.
*
* renvoi la latitude isométrique correspondante
* */
private double calculLatitudeIsometrique(double phi, double e)
{
return Math.Log(Math.Tan((Math.PI / 4) + (phi / 2)) * Math.Pow((1 - e * Math.Sin(phi)) / (1 + e * Math.Sin(phi)),e / 2));
}

/**
* L : latitude isométrique.
* e : première excentricité de l’ellipsoïde.
* epsilon : tolérance de convergence.
*
* renvoi:
* phi : latitude en radian.
**/
private double calculLatitudeAPartirLatitudeIsometrique(double L, double e, double epsilon)
{
double lambdaA = 2 * Math.Atan(Math.Exp(L)) - Math.PI / 2;
double lambda1 = lambda0 + 1;
while(Math.Abs(lambdaA - lambda1) >= epsilon)
{
lambda1 = lambdaA;
lambdaA = 2 * Math.Atan(Math.Pow((1 + e * Math.Sin(lambda1)) / (1 - e * Math.Sin(lambda1)), e / 2) * Math.Exp(L)) - Math.PI / 2;
}
return lambdaA;
}

/**
* lambda : longitude par rapport au méridien origine.
* phi : latitude.
* n : exposant de la projection.
* c : constante de la projection.
* e : première excentricité de l’ellipsoïde.
* lambdaC : longitude de l’origine par rapport au méridien origine.
* Xs, Ys : coordonnées en projection du pôle.
*
* renvoi:
* resultat[0] -> X
* resultat[1] -> Y
* */
private double[] conversionCoordGeographiquesVersCoordLambert(doubl e lambda, double phi, double n, double c, double e, double lambdaC, double Xs, double Ys)
{
double[] resultat = new double[2];
double latIso = calculLatitudeIsometrique(phi, e);
resultat[0] = Xs + c * Math.Exp(-n * latIso) * Math.Sin(n * (lambda - lambdaC));
resultat[1] = Ys - c * Math.Exp(-n * latIso) * Math.Cos(n * (lambda - lambdaC));
return resultat;
}

/**
* X, Y : coordonnées en projection conique conforme de Lambert du point.
* n : exposant de la projection.
* c : constante de la projection.
* e : première excentricité de l’ellipsoïde.
* lambdaC : longitude de l’origine par rapport au méridien origine.
* Xs, Ys : coordonnées en projection du pôle.
* epsilon : tolérance de convergence
*
* renvoi:
* resultat[0] -> lambda : longitude par rapport au méridien origine.
* resultat[1] -> phi : latitude.
**/
private double[] conversionCoordLambertVersCoordLambertGeographique (double X, double Y, double n, double c, double e, double lambdaC, double Xs, double Ys, double epsilon)
{
double[] resultat = new double[2];
double R = Math.Sqrt(Math.Pow((X - Xs), 2) + Math.Pow((Y - Ys), 2));
label6.Text = Ys.ToString();
double gamma = Math.Atan((X - Xs) / (Ys - Y));
resultat[0] = lambdaC + gamma / n;
double L = (1 / -n) * Math.Log(Math.Abs(R / c));
label5.Text = Xs.ToString();
resultat[1] = calculLatitudeAPartirLatitudeIsometrique(L, e, epsilon);
return resultat;
}

/**
* X, Y, Z: coordonnées cartésiennes en WGS84
*
* renvoi:
* resultat[0] -> X coordonnée cartésienne NTF
* resultat[1] -> Y coordonnée cartésienne NTF
* resultat[2] -> Z coordonnée cartésienne NTF
* */
private double[] translationCoordCartWGS84VersNTF(double X, double Y, double Z)
{
double[] resultat = new double[3];
resultat[0] = X + 168.0;
resultat[1] = Y + 60.0;
resultat[2] = Z - 320.0;
return resultat;
}

/**
* X, Y, Z: coordonnées cartésiennes en Lambert
*
* renvoi:
* resultat[0] -> X coordonnée cartésienne WGS84
* resultat[1] -> Y coordonnée cartésienne WGS84
* resultat[2] -> Z coordonnée cartésienne WGS84
* */
private double[] translationCoordCartNTFVersWGS84(double X, double Y, double Z)
{
double[] resultat = new double[3];
resultat[0] = X - 168.0;
resultat[1] = Y - 60.0;
resultat[2] = Z + 320.0;
return resultat;
}

private void initilisationDonneesLambert()
{
if(comboBox1.SelectedIndex == 1 || comboBox1.SelectedIndex == 7) //si Lambert I
{
n = 0.7604059656;
c = 11603796.98;
Xs = 600000.0;
Ys = 5657616.674;
}
else
if(comboBox1.SelectedIndex == 2 || comboBox1.SelectedIndex == 8) //si Lambert II
{
n = 0.7289686274;
c = 11745793.39;
Xs = 600000.0;
Ys = 6199695.768;
}
else
if(comboBox1.SelectedIndex == 3 || comboBox1.SelectedIndex == 9) //si Lambert II étendu
{
n = 0.7289686274;
c = 11745793.39;
Xs = 600000.0;
Ys = 8199695.768;
}
else
if(comboBox1.SelectedIndex == 4 || comboBox1.SelectedIndex == 10) //si Lambert III
{
n = 0.6959127966;
c = 11947992.52;
Xs = 600000.0;
Ys = 6791905.085;
}
else
if(comboBox1.SelectedIndex == 5 || comboBox1.SelectedIndex == 11) //si Lambert IV
{
n = 0.6712679322;
c = 12136281.99;
Xs = 234.358;
Ys = 7239161.542;
}
else
if(comboBox1.SelectedIndex == 6 || comboBox1.SelectedIndex == 12) //si Lambert -93
{
n = 0.7256077650;
c = 11754255.426;
Xs = 700000.0;
Ys = 12655612.050;
}
else //par défaut, Lambert II étendu
{
n = 0.7289686274;
c = 11745793.39;
Xs = 600000.0;
Ys = 8199695.768;
}
lambda0 = 0.04079234433198; //correspond à la longitude en radian de Paris (2°20'14.025" E) par rapport à Greenwich
e = 0.08248325676;
a = 6378249.2;
}

private void initilisationDonneesWGS84()
{
double a_w = 6378137.0;
double b_w = 6356752.314;
double e2_w = (a_w*a_w-b_w*b_w)/(a_w*a_w);
nw = -1;
cw = -1;
Xw = 500000;
Yw = 0;
lambda0w = -1;
ew = Math.Sqrt(e2_w);
aw = 6378137.0;
}

private void conversionWGS84VersLambert()
{
initilisationDonneesLambert();
initilisationDonneesWGS84();
double[] coordcartwgs84 = conversionGeographiqueVersCartesienne(longitude, latitude, 0, aw, ew); //he = 0: on considère etre à 0m d'altitude par rapport à l'ellipsoïde
double[] coordcartntf = translationCoordCartWGS84VersNTF(coordcartwgs84[0],coordcartwgs84[1],coordcartwgs84[2]);
double[] coordgeontf = conversionCartesienneVersGeorgraphique(coordcartnt f[0], coordcartntf[1], coordcartntf[2], a, e, epsilon);
double[] solutioncoordlambert = conversionCoordGeographiquesVersCoordLambert(coord geontf[0], coordgeontf[1], n, c, e, lambda0, Xs, Ys);
label3.Text = solutioncoordlambert[0].ToString();
label4.Text = solutioncoordlambert[1].ToString();
}

private void conversionLambertVersWGS84()
{
initilisationDonneesLambert();
initilisationDonneesWGS84();
double[] coordgeontf = conversionCoordLambertVersCoordLambertGeographique (latitude, longitude, n, c, e, lambda0, Xs, Ys, epsilon);
double[] coordcartntf = conversionGeographiqueVersCartesienne(coordgeontf[0], coordgeontf[1], 0, a, e);
double[] coordcartwgs84 = translationCoordCartNTFVersWGS84(coordcartntf[0],coordcartntf[1],coordcartntf[2]);
double[] solutionwgs84 = conversionCartesienneVersGeorgraphique(coordcartwg s84[0], coordcartwgs84[1], coordcartwgs84[2], aw, ew, epsilon);
label3.Text = solutionwgs84[0].ToString();
label4.Text = solutionwgs84[1].ToString();
}

private void button1_Click(object sender, System.EventArgs e)
{
try
{
epsilon = Math.Pow(10,-10); //tolérance de convergence
latitude = Convert.ToDouble(textBox1.Text)*0.0174532925; // conversion en radians
longitude = Convert.ToDouble(textBox2.Text)*0.0174532925;// conversion en radians
typeConversion = comboBox1.SelectedIndex+1;
label3.Text = "Conversion";
label4.Text = "en cours...";
if(this.typeConversion<=6) {
conversionWGS84VersLambert();
}
else {
conversionLambertVersWGS84();
}
}
catch (Exception exc)
{
MessageBox.Show(exc.Message, "Mauvais format de coordonnées", MessageBoxButtons.OK, MessageBoxIcon.Error);
}
}


}
}

Cependant, j'ai un e encore différent des votres, quelqu'un pourrait-il m'éclairer? SVP!!! merciiii

[n314]

Bonjour,
je ressors ce topic en posant la question: si tous vos calculs sont bons:

peut on faire confiance à circe ?
si la réponse est non, d'où vient l'erreur? Y a t il des précautions à prendre ?

Bien cordialement,
n314

[Rimi S]

Salut Monsieur N.

Je suis très intéressé par la réponse, également...
Merci encore, pour ce partage de vos expériences.

[ReadWrite]

Envoyé par n314

peut on faire confiance à circe ?

j'ai utilisé le logiciel Circé pour évaluer sur plusieurs zones de travail la précision de la transformation par grille de l'IGN

la comparaison portait sur l'écart en planimétrie de points observés en wgs 84 et transformés en lambert zone

• la première transformation était une simililtude 3D issue de l'observation de points de la NTF
• la seconde était la grille planimétrique de l'IGN utilisée dans Circé

sur les différentes zones de travail, d'une quinzaine de kms de côté, disséminées sur un seul département,
les écarts étaient tous inférieurs à 10 cm, et la plupart du temps d'environ 5 cm



tout porte à croire selon moi que les résultats donnés par Circé sont corrects
ils l'ont été en tout cas pour les tests décrits ci-dessus


NB : certains gps permettent de stocker les transformations par grille et de les appliquer aux coordonnées wgs 84 observées directement sur le terrain
si l'on souhaite travailler dans un lambert zone de très bonne qualité, sans avoir à déterminer les similitudes, c'est vraiment pas mal

[yjacolin]

Bonsoir,

Merci pour vos contribution. Voici ci-dessous de quoi transformer des coordonnées en WGS84 (de google par exemple) en des coordonnées lambert II étendue. La fonction est écrite en php.


Y.

[Le Docteur]

Merci Yves

[arena]

Envoyé par yjacolin

Bonsoir,

Merci pour vos contribution. Voici ci-dessous de quoi transformer des coordonnées en WGS84 (de google par exemple) en des coordonnées lambert II étendue. La fonction est écrite en php.


Y.

bonjour a tous, je suis complétement débutant en geomatique etc ...

je cherche le script php qui fasse la conversion inverse (qqchose du type lambert2wgs84.php)

a l'avance merci.

ps : j'ai essayé de m'en sortir avec les doc de l'ign mais je nage ...
http://www.ign.fr/telechargement/MPr...CE/transfo.pdf
http://www.ign.fr/telechargement/FAQ/FAQ1.pdf
http://www.ign.fr/telechargement/MPr...RCE/NTG_71.pdf