Le portail francophone de la géomatique

ClementB · Wed 03 July 2013 16:22

Bonjour à tous,

je tente d'exploiter une dalle de l'ortho littorale 2000 (dalle 23, http://geolittoral.application.equipeme … alle23.zip) avec gdal (en python), mais je suis confronté à un décalage au niveau de la projection.

L'ortho littorale est projetée en Lambert II étendu.

La projection déterminée par gdal à partir de l'ecw est la suivante :

Code:

PROJCS["LM2FRANC",
    GEOGCS["N.T.F.",
        DATUM["NTF",
            SPHEROID["CLA80IGN",6378249.2,293.4660213]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",45.89891888888889],
    PARAMETER["standard_parallel_2",47.69601444444444],
    PARAMETER["latitude_of_origin",46.80000000000000],
    PARAMETER["central_meridian",2.337229169999754],
    PARAMETER["false_easting",600000],
    PARAMETER["false_northing",2200000],
    UNIT["Meter",1]]

Code:

+proj=lcc +lat_1=45.89891888888889 +lat_2=47.69601444444444 +lat_0=46.8 +lon_0=2.337229169999754 +x_0=600000 +y_0=2200000 +a=6378249.2 +b=6356515.000000472 +units=m +no_defs

Pour tester je me base sur ce point géodésique (château d'eau): http://geodesie.ign.fr/fiches/pdf/56162C.pdf
lon 3° 23' 51.8943'' O
lat 47° 44' 58.4950'' N

Avec circé:
NFT(Paris) / Lambert II Étendu
x = 170383.255 m
y = 2321251.802 m
précision: 2-5 m
ce qui rapporté sur le raster désigne bien le château d'eau

Avec gdal (ou proj4) on obtient :
x = 170312.760151 m
y = 2321245.57081 m
deltaX: 70.494848759 m
deltaY: 6.23118673917 m

D'où vient ce décalage ?

Je remarque que la projection obtenue à partir de l'EPSG (27572) est différente :

Code:

PROJCS["NTF (Paris) / Lambert zone II",
    GEOGCS["NTF (Paris)",
        DATUM["Nouvelle_Triangulation_Francaise_Paris",
            SPHEROID["Clarke 1880 (IGN)",6378249.2,293.4660212936265,
                AUTHORITY["EPSG","7011"]],
            TOWGS84[-168,-60,320,0,0,0,0],
            AUTHORITY["EPSG","6807"]],
        PRIMEM["Paris",2.33722917,
            AUTHORITY["EPSG","8903"]],
        UNIT["grad",0.01570796326794897,
            AUTHORITY["EPSG","9105"]],
        AUTHORITY["EPSG","4807"]],
    PROJECTION["Lambert_Conformal_Conic_1SP"],
    PARAMETER["latitude_of_origin",52],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",0.99987742],
    PARAMETER["false_easting",600000],
    PARAMETER["false_northing",2200000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    AUTHORITY["EPSG","27572"]]

Code:

+proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742 +x_0=600000 +y_0=2200000 +a=6378249.2 +b=6356515 +towgs84=-168,-60,320,0,0,0,0 +pm=paris +units=m +no_defs

De plus cette projection me fourni les bons résultats.

Y-a-t'il une erreur dans le fichier ecw ? Ce serait étonnant, je dois donc rater quelque chose.
Savez vous d'où vient ce problème ?

Merci

ps: je joins le fichier de test en python
[img]http://img20.imageshack.us/img20/6342/8nri.png[/img]

Dernière modification par ClementB (Wed 03 July 2013 23:39)

Fichier(s) joint(s) :
Pour accéder aux fichiers vous devez vous inscrire.

dominique.lys · Wed 03 July 2013 18:27

Salut,

Ce sont exactement les mêmes projections sauf que l'une est exprimée en tant que conique tangente et l'autre en tant conique sécante mais ça revient au même.

Par contre dans la définition proj4 issue du ECW il faut que tu ajoutes +towgs84=-168,-60,320,0,0,0,0, ce paramètre permet de calculer correctement les translations entre datum differents (de WGS84 pour les coordonnées geo vers NTF pour le lambert 2).

Vérification:

Code:

### -*- coding:Latin-1 -*-
from osgeo import osr

geo = osr.SpatialReference()
geo.ImportFromEPSG(4326)

lb2 = osr.SpatialReference()
lb2.ImportFromEPSG(27572)

transfo = osr.CoordinateTransformation(geo,lb2)
result=transfo.TransformPoint(-3.39774842, 47.74958194)
print(result)#170383.25519322767, 2321251.801505601 (idem circé)

#Avec la définition de projection du fichier ECW
lb2 = osr.SpatialReference()
lb2.ImportFromProj4("+proj=lcc +lat_1=45.89891888888889 +lat_2=47.69601444444444 +lat_0=46.8 +lon_0=2.337229169999754 +x_0=600000 +y_0=2200000 +a=6378249.2 +b=6356515.000000472 +towgs84=-168,-60,320,0,0,0,0 +units=m +no_defs")
transfo = osr.CoordinateTransformation(geo,lb2)
result=transfo.TransformPoint(-3.39774842, 47.74958194)
print(result)#170383.2549358232, 2321251.8015073068 --> précision centimètrique

Dernière modification par dominique.lys (Wed 03 July 2013 18:28)

jcr83 · Wed 03 July 2013 22:10

Avec PROJ4, j'obtiens exactement le même résultat que Circé :
$ cs2cs -f "%.3f" +init=epsg:4326 +to +init=epsg:27572
3d23'51.8943"W 47d44'58.4950"N
170383.255 2321251.802 -47.488

Pour plus de précision, utilise une transformation à grille :
$ cs2cs -f "%.3f" +init=epsg:4326 +to +init=IGNF:LAMB2C
3d23'51.8943"W 47d44'58.4950"N
170382.771 2321247.719 0.000

Ces coordonnées devraient être encore plus proches du château d'eau. Ce serait intéressant si tu pouvais poster la photo, pour comparer.

ClementB · Wed 03 July 2013 23:38

Merci pour ces précisions, ça marche effectivement mieux avec la transformation de datum. Dommage que cette information ne soit pas contenue dans l'ecw.

jcr83, voici la photo avec les coordonnées obtenue avec la transformation à grille:
[img]http://img16.imageshack.us/img16/4671/r2c.png[/img]
c'est en effet plus au centre du château d'eau.

jcr83 · Thu 04 July 2013 00:08

Tu devrais ouvrir un ticket sur http://trac.osgeo.org/gdal/newticket
Pour info, si tu veux l'ajouter toi-même, c'est dans le fichier ecw_cs.wkt qu'il faut le faire.

Communauté	Ressources	L'emploi	Le marché	L'association	Fils RSS
Forum Blogs Wiki Planet Sites amis	Agenda Biblio Fiches AMO-CNIG Paris PCGIS Téléchargements	Offres d'emploi Banque de CV Formations Démarche Métiers	Geo-Entreprises Geo-Communiqué Appels d'offres	L'association Qui sommes-nous ? Nous aider Bienfaiteurs Partenaires Historique Statuts	Tout le site Job Géomatique Services GeoBlogs ...	Faire un don Mentions légales

Le portail francophone de la géomatique

Annonce

#1 Wed 03 July 2013 16:22

Ortho littorale, Lambert II étendu et gdal/proj4: décallage

Code:

Code:

Code:

Code:

#2 Wed 03 July 2013 18:27

Re: Ortho littorale, Lambert II étendu et gdal/proj4: décallage

Code:

#3 Wed 03 July 2013 22:10

Re: Ortho littorale, Lambert II étendu et gdal/proj4: décallage

#4 Wed 03 July 2013 23:38

Re: Ortho littorale, Lambert II étendu et gdal/proj4: décallage

#5 Thu 04 July 2013 00:08

Re: Ortho littorale, Lambert II étendu et gdal/proj4: décallage

Pied de page des forums