Como georreferenciar um bloco de mercador da web corretamente usando gdal?

Como exemplo, vou pegar o seguinte bloco http://a.tile.openstreetmap.org/3/4/2.png e salvá-lo como "4_2.png".

As coordenadas WGS84 desta peça pode ser calculado ou ler lá clicando na telha correspondente:

0 66.51326044311185 45 40.97989806962013 (West North East South)

Como georreferenciar o bloco corretamente (usando gdal para gerar um geotiff ou outro formato georreferenciado) para que:

• O bitmap não precisa ser esticado (= os pixels no geotiff são exatamente iguais aos do bitmap original)
• A imagem resultante será aberta no local certo em um visualizador / editor de GIS (como por exemplo no TatukGIS Free Viewer )?

(Editado em 19 de setembro de 2011 para tornar minha pergunta mais clara e incluir minhas conclusões)

Minha conclusão:

Primeiro pensei que a terceira ideia (veja abaixo) era a correta. Abri o geotiff em um GIS Viewer e comparei as coordenadas exibidas com o que eu esperava. O geotiff fora da segunda idéia parece ser deslocado 2 pixels para o norte. Por isso, considerei a idéia 3 (ou 4) a correta.
Mas se você tentar com um bloco com um nível de zoom muito maior, o geotiff da idéia 3 será definitivamente deslocado para o sul. Foi tolo comparar coordenadas em um bloco de nível de zoom 3. Os limites do país nesse nível de zoom são simplificados para que a comparação não dê bons resultados.

Dan S. estava certo, a imagem do bloco já está no EPSG: 3857. A segunda idéia é a correta (e também fornece bons resultados em altos níveis de zoom)

Primeira idéia: EPSG: 4326
O código EPSG para as coordenadas WGS84 é EPSG: 4326 . Então, eu simplesmente uso as coordenadas WGS84 para georreferenciar o bloco como geotif usando 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

O mapa resultante é exibido no lugar certo, mas temo que a projeção não esteja correta e que possa haver uma mudança no meio do bloco. Depois de muito tempo para verificar que, reprojetando o mapa com gdalwarp, baixei uma versão demo do Global Mapper e parece ser esse o caso (sames border como na idéia 3, mas uma mudança dentro do tile). A imagem deve ser esticada para poder usar as coordenadas EPSG: 4326.

Segunda ideia: EPSG: 3857
Este bloco usa uma projeção "web mercator" (também conhecida como projeção do google map), que agora possui um código EPSG: EPSG: 3857 (também conhecido por EPSG: 900913). Simplesmente converto as coordenadas usando 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

Minhas coordenadas em metros são:

0 10018754.1713946 5009377.08569731 5009377.08569731 (West North East South)

Agora eu posso usar gdal_translate para gerar um geotiff:

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

Minha impressão é que isso não está correto, porque as bordas dos mapas estão voltadas para o norte. Parece ser a ideia certa.

Terceira idéia: EPSG: 3857 a EPSG: 4055
Li que "web mercator" usa coordenadas WGS84, mas as considero como se fossem coordenadas esféricas. Devido à diferença entre uma latitude geodésica e geocêntrica (consulte a Wikipedia sobre latitude ), os valores de latitude não serão os mesmos em um elipsóide ou em uma esfera. Descobri que o EPSG: 4055 é o código para coordenadas esféricas em uma esfera baseada no WGS84.

Convertendo as coordenadas em 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

As coordenadas esféricas correspondentes são então:

0 66.3722684317026 45 40.7894557844857 (West North East South)

Então, faço como se essas coordenadas ainda estivessem no elipsóide (EPSG: 4326) e as converti em 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

As coordenadas resultantes diferem da coordenada idea2:

0 9979483.26733298 5009377.08569731 4981335.86590183 (West North East South)

Agora só tenho que escrever as coordenadas no mapa:

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

Essa terceira idéia parece dar os melhores resultados. Mas não tenho certeza se está correto. Se a idéia 3 estiver correta, existe um código EPSG para executar esta operação em uma única etapa?

Quarta ideia: EPSG: 3857 através de towgs84 = 0,0,0,0,0,0,0

gdal (e aparentemente epsg também) define o EPSG: 3857 assim:

+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

considerando spatialreference.org assim:

+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

Se eu usar a definição de spatialreference.org, obtive as coordenadas corretas em uma etapa (ainda não sei se são as coordenadas "corretas", mas pelo menos são as mesmas da idéia 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

Por que existe tanta diferença nas definições de EPSG: 3857?

A imagem do bloco já está em EPSG: 3857. Por que não criar um arquivo mundial para fazer referência a ele?

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

Para o bloco que cobre a América do Norte no zoom 1, você deve observar o seguinte conteúdo do arquivo mundial:

78271.517
0
0
-78271.517
-19998372.6
19998372.6

De onde vieram esses números:

• Linha 1: largura de um pixel de imagem em coordenadas mundiais = 20037508.342789244 metros / 256 pixels.
• Linha 2 e 3: rotação, então n / a.
• Linha 4: altura de um pixel de imagem nas coordenadas do mundo. Igual à linha 1, mas negativa, porque nos arquivos de imagem o aumento de y corresponde a 'inativo' enquanto no sistema de coordenadas, o aumento de y corresponde a 'up'.
• Linha 5: coordenada X em coordenadas mundiais do centro do pixel superior esquerdo. Isso é -20037508.342789244, conforme relatado pelos links à la carte dos blocos, mais 1/2 de pixel para trazê-lo para o centro.
• Linha 6: Idem, apenas a coordenada Y do canto superior esquerdo.

O GDAL deve pegar seu arquivo mundial (.pgw para o png); você ainda precisará informar o EPSG: 3857 na linha de comando.

(Nota: não teve tempo de testar isso, então está tudo fora do manguito ... mas espero que esteja correto na primeira tentativa de qualquer maneira!;)

Funções necessárias para a geração de blocos globais usados ​​na web. Ele contém classes que implementam conversões de coordenadas para:

• GlobalMercator (baseado em EPSG: 900913 = EPSG: 3785 ) para blocos compatíveis com Google Maps, Yahoo Maps, Microsoft Bing Maps

• GlobalGeodetic (baseado em EPSG: 4326) para o OpenLayers Base Map e blocos compatíveis com Google Earth

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

Sobre minha pergunta secundária nas quarta idéias: Por que existe tanta diferença nas definições de EPSG: 3857 entre a definição em gdal e em spatialreference.org :

A principal diferença é que o gdal usa "+ nadgrids = @ null" e referência espacial "+ towgs84 = 0,0,0,0,0,0,0". De acordo com a documentação ou o formato PROJ.4, ambos os parâmetros são usados ​​para transformações de dados. Encontrei um comentário interessante de Mikael Rittri no servidor de lista do 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."

O uso de "+ towgs84 = 0,0,0,0,0,0,0" não parece correto ou pelo menos apenas sob condições específicas.