Comment géoréférencer correctement une tuile web mercator en utilisant gdal?

16

À titre d'exemple, je prendrai la vignette suivante http://a.tile.openstreetmap.org/3/4/2.png et l'enregistrer sous "4_2.png".

Les coordonnées WGS84 de cette tuile peuvent être calculées ou lire là en cliquant sur le carreau correspondant:

0 66.51326044311185 45 40.97989806962013 (West North East South)

Comment géoréférencer correctement la tuile (en utilisant gdal pour générer un géotiff ou un autre format géoréférencé) afin que:

  • Le bitmap n'a pas besoin d'être étiré (= les pixels dans le géotiff sont exactement les mêmes que dans le bitmap d'origine)
  • L'image résultante sera ouverte au bon endroit dans un visualiseur / éditeur SIG (comme par exemple dans TatukGIS Free Viewer )?

(Modifié le 19 septembre 2011 pour clarifier ma question et inclure mes conclusions)

Ma conclusion:

J'ai d'abord pensé que la troisième idée (voir ci-dessous) était la bonne. J'ai ouvert le géotiff dans une visionneuse SIG et comparé les coordonnées affichées avec ce que j'attendais. Le géotiff de la deuxième idée semble décalé de 2 pixels vers le nord. C'est pourquoi j'ai considéré l'idée 3 (ou 4) comme la bonne.
Mais si vous essayez avec une tuile à un niveau de zoom beaucoup plus élevé, le géotiff de l'idée 3 est définitivement décalé vers le sud. Il était stupide de comparer les coordonnées sur une tuile de niveau de zoom 3. Les frontières des pays à un tel niveau de zoom sont mus simplifiées afin que la comparaison ne donne pas de bons résultats.

Dan S. avait raison, l'image de la tuile est déjà en EPSG: 3857. La deuxième idée est alors la bonne (et donne également de bons résultats à des niveaux de zoom élevés)

Première idée: EPSG: 4326
Le code EPSG pour les coordonnées WGS84 est EPSG: 4326 . J'utilise donc simplement les coordonnées WGS84 pour géoréférencer la tuile en tant que géotif en utilisant gdal_translate :

gdal_translate -of Gtiff -co tfw=yes -a_ullr 0 66.51326044311185 45 40.97989806962013 -a_srs EPSG:4326 4_2.png t4326.tif

La carte résultante est affichée au bon endroit mais je crains que la projection ne soit pas correcte et qu'il puisse y avoir un décalage au milieu de la tuile. Après avoir longtemps essayé de vérifier qu'en reprojetant la carte avec gdalwarp, j'ai téléchargé une version de démonstration de Global Mapper et cela semble être le cas (même frontières que par l'idée 3 mais un décalage à l'intérieur de la tuile). L'image doit être étirée pour pouvoir utiliser les coordonnées EPSG: 4326.

Deuxième idée: EPSG: 3857
Cette tuile utilise une projection "web mercator" (alias projection google map), qui a maintenant un code EPSG: EPSG: 3857 (alias EPSG: 900913). Je convertis simplement les coordonnées en utilisant gdaltransform :

gdaltransform -s_srs EPSG:4326 -t_srs EPSG:3857
0 66.51326044311185
0 10018754.1713946 0
45 40.97989806962013
5009377.08569731 5009377.08569731 0

Mes coordonnées en mètres sont: 

0 10018754.1713946 5009377.08569731 5009377.08569731 (West North East South)

Maintenant, je peux utiliser gdal_translate pour générer un géotiff:

gdal_translate -of Gtiff -co tfw=yes -a_ullr 0 10018754.1713946 5009377.08569731 5009377.08569731 -a_srs EPSG:3857 4_2.png t3857.tif

J'ai l'impression que ce n'est pas correct car les frontières des cartes sont décalées vers le nord. Cela semble être la bonne idée.

Troisième idée: EPSG: 3857 à EPSG: 4055
J'ai lu que "web mercator" utilise les coordonnées WGS84 mais les considère comme si elles étaient des coordonnées sphériques. En raison de la différence entre une latitude géodésique et une latitude géocentrique (voir Wikipedia sur la latitude ), les valeurs de latitude ne seront pas les mêmes sur un ellipsoïde ou sur une sphère. J'ai trouvé que EPSG: 4055 est le code pour les coordonnées sphériques sur une sphère basée sur WGS84.

Conversion des coordonnées en EPSG: 4055:

gdaltransform -s_srs EPSG:4326 -t_srs EPSG:4055
0 66.51326044311185
0 66.3722684317026 -17964.0621483233
45 40.97989806962013
45 40.7894557844857 -9152.84527519904

Les coordonnées sphériques correspondantes sont alors:

0 66.3722684317026 45 40.7894557844857 (West North East South)

Ensuite, je fais comme si ces coordonnées étaient toujours sur l'ellipsoïde (EPSG: 4326) et les convertis en web mercator:

gdaltransform -s_srs EPSG:4326 -t_srs EPSG:3857
0 66.3722684317026
0 9979483.26733298 0
45 40.7894557844857
5009377.08569731 4981335.86590183 0

Les coordonnées résultantes diffèrent de celles de idea2:

0 9979483.26733298 5009377.08569731 4981335.86590183 (West North East South)

Il ne me reste plus qu'à écrire les coordonnées sur la carte:

gdal_translate -of Gtiff -co tfw=yes -a_ullr 0 9979483.26733298 5009377.08569731 4981335.86590183 -a_srs EPSG:3857 4_2.png t3857_through_4055.tif

Cette troisième idée semble donner les meilleurs résultats. Mais je ne sais pas si c'est correct. Si l'idée 3 est correcte, existe-t-il un code EPSG pour effectuer cette opération en une seule étape?

Quatrième idée: EPSG: 3857 à travers towgs84 = 0,0,0,0,0,0,0

gdal (et apparemment epsg aussi) définit EPSG: 3857 comme ceci:

+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

alors que spatialreference.org comme ça:

+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

Si j'utilise la définition de spatialreference.org, j'ai obtenu les coordonnées correctes en une seule étape (Eh bien, je ne le fais toujours pas si ce sont les coordonnées "correctes" mais au moins ce sont les mêmes que dans l'idée 3):

gdaltransform -s_srs EPSG:4326 -t_srs "+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
0 66.51326044311185
0 9979483.26733298 -17964.0621483233
45 40.97989806962013
5009377.08569731 4981335.86590183 -9152.84527519904

Pourquoi y a-t-il une telle différence dans les définitions de l'EPSG: 3857?

web-mapping gdal epsg tiles Nom
la source

Réponses:

11

L'image de tuile est déjà en EPSG: 3857. Pourquoi ne pas simplement créer un fichier mondial pour le référencer?

http://en.wikipedia.org/wiki/World_file

Pour la tuile qui couvre l'Amérique du Nord au zoom 1, vous regarderiez le contenu du fichier mondial suivant:

78271.517
0
0
-78271.517
-19998372.6
19998372.6

D'où provenaient ces chiffres:

  • Ligne 1: largeur d'un pixel d'image en coordonnées universelles = 20037508.342789244 mètres / 256 pixels.
  • Ligne 2 et 3: rotation, donc n / a.
  • Ligne 4: hauteur d'un pixel d'image en coordonnées universelles. Identique à la ligne 1 mais négatif, car dans les fichiers image, l'augmentation de y correspond à «bas» tandis que dans le système de coordonnées, l'augmentation de y correspond à «haut».
  • Ligne 5: coordonnée X en coordonnées universelles du centre du pixel supérieur gauche. Il s'agit de -20037508.342789244, comme indiqué par le lien à la carte des tuiles, plus 1/2 de pixel pour l'amener au centre.
  • Ligne 6: Idem, seule la coordonnée Y en haut à gauche.

GDAL devrait récupérer votre fichier mondial (.pgw pour le png); vous devrez toujours lui indiquer EPSG: 3857 sur la ligne de commande.

(Remarque: je n'ai pas eu le temps de tester cela, donc tout est hors du jeu ... mais j'espère que tout sera correct au premier essai!;)

Dan S.
la source
Merci et désolé pour la réponse tardive. Mais en fait, je pense que cela conduirait à la même chose que ma deuxième idée, où gdaltransform n'est vraiment utilisé que pour géoréférencer l'image.
Nom du
Vous aviez raison sur l'image de tuile déjà en EPSG: 3857. La solution était en fait d'utiliser simplement EPSG: 3857. C'est la façon dont je vérifiais les résultats qui était faux.
Nom du
0

Fonctions nécessaires à la génération de tuiles globales utilisées sur le web. Il contient des classes implémentant des conversions de coordonnées pour:

  • GlobalMercator (basé sur EPSG: 900913 = EPSG: 3785 ) pour les cartes compatibles Google Maps, Yahoo Maps, Microsoft Bing Maps

  • GlobalGeodetic (basé sur EPSG: 4326) pour OpenLayers Base Map et les tuiles compatibles avec Google Earth

http://svn.osgeo.org/gdal/sandbox/klokan/gdal2tiles.py

Mapperz
la source
Merci mais cela ne répond pas à ma question. Je ne veux pas générer de tuiles. Je veux savoir quel est / serait le géoréférencement correct des tuiles.
Nom du
J'ai fait "Si c'est correct, existe-t-il un code EPSG pour effectuer cette opération en une seule étape?" Réponse EPSG: 900913
Mapperz
Vous ne l'avez pas fait: EPSG: 900913 fonctionne comme EPSG: 3857 (ou comme EPSG: 3785) et ne permet pas de faire l'opération en une seule étape, vous devez toujours parcourir EPSG: 4055. De plus, j'avais déjà mentionné EPSG: 900913 comme alias pour EPSG: 3857 dans ma question initiale.
Nom du
0

À propos de ma question secondaire dans les quatrièmes idées: Pourquoi y a-t-il une telle différence dans les définitions d'EPSG: 3857 entre la définition dans gdal et sur spatialreference.org :

La principale différence est que gdal utilise "+ nadgrids = @ null" et la référence spatiale "+ towgs84 = 0,0,0,0,0,0,0". Selon la documentation ou le format PROJ.4, les deux paramètres sont utilisés pour les transformations de datum. J'ai trouvé un commentaire intéressant de Mikael Rittri sur le serveur de liste Proj4 :

"The +to definition you suggest is not correct for WGS84 Web Mercator, though.
This is because the datum shift you use, 

   +towgs84=0,0,0,0,0,0,0

is not the same as unmodified lat/long, or "without any datum shift" as you say.
It means "no datum shift" only if the +from and +to definitions use the same ellipsoid,
and in your case the +from definition uses the WGS84 ellipsoid, while the +to definition
uses a sphere instead.  

So, you have to use a trick: use 

   +nadgrids=@null 

instead."

L'utilisation de "+ towgs84 = 0,0,0,0,0,0,0,0" ne semble pas être correcte ou du moins uniquement dans des conditions spécifiques.

Nom
la source