Discussion :

[vbain]

Bonjour à tous,

la question qui m'amène aujourd'hui pourra paraître pénible à suivre, en tout cas je crains qu'elle ne soit pas simple à présenter de manière concise (si ça peut faire gagner du temps aux personnes qui courent après... leur temps, et qui ne voudront pas suivre plus avant).

_______________________
1/ Présentation du problème

Je travaille à l'élaboration d'un requêteur de données raster (flux wmts du geoportail) pour le SIG GRASS. Le module dans sa forme actuelle fonctionne a priori normalement, mais dans le détail, certaines imprécisions me font douter quant à sa fiabilité. Pour vérifier la rigueur géodésique des opérations qu'il effectue, ma démarche est la suivante :
- je récupère la fiche géodésique d'un point facilement identifiable sur l'orthophotoplan (dans la limite de la meilleure résolution disponible bien sûr, soit le niveau de zoom 18 de la pyramide). Dans ce cas il s'agit d'un pylône électrique ;
- sur la page d'accueil du Géoportail, je m'assure que les coordonnées du point visé sont cohérentes avec celles fournies par la fiche géodésique ;
- je récupère sur le serveur wmts la tuile où figure le lieu examiné ;
- je procède au géoréférencement de cette tuile, dans un premier temps simplement dans le système natif de projection (pseudo-mercator).

Je constate que les coordonnées du point géodésique :
(i) fournies par l'interface de geoportail.gouv.fr,
(ii) mesurées sur l'image géoréférencée par mes soins
diffèrent de quelques mètres, les premières étant correctes, les secondes divergentes.

Mon problème se rapporte donc à une question de géoréférencement des tuiles du flux wmts.

__________________
2/ Méthode employée

chargement de la tuile :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

wget "https://$user:$pwd@gpp3-wxs.ign.fr/$clef/geoportail/wmts?LAYER=ORTHOIMAGERY.ORTHOPHOTOS&EXCEPTIONS=text/xml&FORMAT=image/jpeg&SERVICE=WMTS&VERSION=1.0.0&REQUEST=GetTile&STYLE=normal&TILEMATRIXSET=PM&TILEMATRIX=18&TILEROW=93501&TILECOL=135775&" -O tilel93501c135775.jpg

détermination des coordonnées de 3 coins de l'image / géoréferencement :
pour être bref, je suis le raisonnement de cette page. Pour la tuile de notre exemple, j'obtiens trois points de contrôle que je fournis à l'utilitaire gdal_translate comme suit :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

gdal_translate -gcp 0 0 718957.177600 5743637.622016 -gcp 256 256 719110.051584 5743484.748032 -gcp 256 0 719110.051584 5743637.622016 tilel93501c135775.jpg tilel93501c135775.tif

Le script gcps2wld permet de créer éventuellement le worldfile associé :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

gcps2wld.py tilel93501c135775.tif>tilel93501c135775.tfw

__________
3/ Résultats

Pour ceux qui m'ont suivi jusque là, je fournis ici en image le résultat. Le point géodésique est situé au centre du pylône marqué par un point jaune. Ses coordonnées en Lambert93 sont les suivantes : 968705.11 6524806.43 (fiche disponible ici). J'obtiens ses coordonnées en pseudo-mercator comme suit :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
cs2cs +init=IGNF:LAMB93 +to +init=epsg:3857
968705.11 6524806.43
718968.86	5743551.33 -0.00

La copie d'écran montre deux fenêtres :
- à droite les coordonnées mercator fournies sur la page geoportail.gouv.fr ;
- à gauche les mêmes coordonnées dans le bas du moniteur de GRASS, correspondant à la position du curseur, soit une petite dizaine de mètres à l'est du point attendu...

Voilà, j'en ai fini ! et me tiens à l'écoute de vos réactions/suggestions pour trouver l'erreur.
Bonne fin de journée,

Vincent.

[mga_geo]

Pour géoréférencer l'image je procède autrement:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
gdal_translate -a_srs EPSG:3857 -co compress=lzw -a_ullr "718967.03" "5743630.84" "719119.90" "5743477.96" GP_PVA_18.jpg GP_PVA_18.tif

je passe ensuite en Lambert-93 avec gdalwarp

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
gdalwarp -co compress=lzw -t_srs EPSG:2154 GP_PVA_18.tif GP_PVA_18_L93.tif

Les gdalinfo donnent:

+ gdalinfo GP_PVA_18.tif
Driver: GTiff/GeoTIFF
Files: GP_PVA_18.tif
Size is 256, 256
Coordinate System is:
PROJCS["WGS 84 / Pseudo-Mercator",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]],
PROJECTION["Mercator_1SP"],
PARAMETER["central_meridian",0],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],
AUTHORITY["EPSG","3857"]]
Origin = (718967.030000000030000,5743630.839999999900000)
Pixel Size = (0.597148437499982,-0.597187499999563)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left ( 718967.030, 5743630.840) ( 6d27'30.93"E, 45d46'13.42"N)
Lower Left ( 718967.030, 5743477.960) ( 6d27'30.93"E, 45d46' 9.97"N)
Upper Right ( 719119.900, 5743630.840) ( 6d27'35.87"E, 45d46'13.42"N)
Lower Right ( 719119.900, 5743477.960) ( 6d27'35.87"E, 45d46' 9.97"N)
Center ( 719043.465, 5743554.400) ( 6d27'33.40"E, 45d46'11.70"N)
Band 1 Block=256x10 Type=Byte, ColorInterp=Red
Band 2 Block=256x10 Type=Byte, ColorInterp=Green
Band 3 Block=256x10 Type=Byte, ColorInterp=Blue

gdalinfo GP_PVA_18_L93.tif
Driver: GTiff/GeoTIFF
Files: GP_PVA_18_L93.tif
Size is 267, 267
Coordinate System is:
PROJCS["RGF93 / Lambert-93",
GEOGCS["RGF93",
DATUM["Reseau_Geodesique_Francais_1993",
SPHEROID["GRS 1980",6378137,298.2572221010002,
AUTHORITY["EPSG","7019"]],
TOWGS84[0,0,0,0,0,0,0],
AUTHORITY["EPSG","6171"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4171"]],
PROJECTION["Lambert_Conformal_Conic_2SP"],
PARAMETER["standard_parallel_1",49],
PARAMETER["standard_parallel_2",44],
PARAMETER["latitude_of_origin",46.5],
PARAMETER["central_meridian",3],
PARAMETER["false_easting",700000],
PARAMETER["false_northing",6600000],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AUTHORITY["EPSG","2154"]]
Origin = (968701.412163856670000,6524866.320167277900000)
Pixel Size = (0.416225728463330,-0.416225728463330)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left ( 968701.412, 6524866.320) ( 6d27'30.94"E, 45d46'13.57"N)
Lower Left ( 968701.412, 6524755.188) ( 6d27'30.71"E, 45d46' 9.97"N)
Upper Right ( 968812.544, 6524866.320) ( 6d27'36.08"E, 45d46'13.41"N)
Lower Right ( 968812.544, 6524755.188) ( 6d27'35.85"E, 45d46' 9.82"N)
Center ( 968756.978, 6524810.754) ( 6d27'33.39"E, 45d46'11.69"N)
Band 1 Block=267x10 Type=Byte, ColorInterp=Red
Band 2 Block=267x10 Type=Byte, ColorInterp=Green
Band 3 Block=267x10 Type=Byte, ColorInterp=Blue

Pour vérifier les coordonnées du pylône, j'ai utilisé IGNMap, et cela me semble tout bon !

[tcoupin]

La résolution que vous avez utilisé (0.597164) est donnée avec assez de chiffres après la virgules pour que ça fonctionne dans le sens Coord->TileCOL/ROW.

En revanche, dans l'autre sens, on faut faire attention et recalculer la résolution. Il faut utiliser le scaleDenominator donné dans le GetCapabilities. Pour le niveau 18 en EPSG:3857, il vaut 1:2132.7295838497840572. La norme WMTS affirme que un pixel vaut 0.28 mm.
On obtient donc une résolution de 0.28/1000*2132.72958[...]572 = 0.597164283477940 m/px

L'écart entre la résolution utilisée et la résolution recalculée, multiplié par 256 pixels, multiplié par les 135775 colonnes donne bien un écart de 10 m environ.

On obtient bien les coordonnées données par mga_geo :
topleft : (718967.031, 5743630.837)
bottomrigth : (719119.905, 5743477.963)

[vbain]

Bonjour mga_geo,

merci pour votre réponse. Je ne connaissais pas l'option -a_ullr de gdal_translate. C'est en regardant vos coordonnées, différentes de quelques mètres des miennes, que j'ai compris mon erreur : mon script calcule la résolution de chaque pyramide en suivant la formule :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

resol=resolzero / 2^tilematrix

avec
resolzero=156543.033928,
Je m'étais contenté d'un résultat à 6 chiffres après la virgule*. D'où cette erreur de quelques mètres...

Merci, et bonne journée.
Vincent

*
écrit en bash ça donnait :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

resol=`echo "scale=6;${resolzero} / 2^${tilematrix}" | bc`

[vbain]

tcoupin, désolé : nos réponses se sont croisées !
on est bien d'accord.
Merci à vous.

Vincent.