¿Cómo georreferenciar un mosaico web mercator correctamente usando gdal?

Como ejemplo, tomaré el siguiente mosaico http://a.tile.openstreetmap.org/3/4/2.png y lo guardaré como "4_2.png".

Las coordenadas WGS84 de este mosaico se pueden calcular o leer allí haciendo clic en el mosaico correspondiente:

0 66.51326044311185 45 40.97989806962013 (West North East South)

Cómo georreferenciar el mosaico correctamente (usando gdal para generar un geotiff u otro formato georreferenciado) para que:

• No es necesario estirar el mapa de bits (= los píxeles en el geotiff son exactamente los mismos que en el mapa de bits original)
• ¿La imagen resultante se abrirá en el lugar correcto en un visor / editor SIG (como, por ejemplo, en TatukGIS Free Viewer )?

(Editado el 19 de septiembre de 2011 para aclarar mi pregunta e incluir mis conclusiones)

Mi conclusión:

Primero pensé que la tercera idea (ver más abajo) era la correcta. Abrí el geotiff en un visor SIG y comparé las coordenadas mostradas con lo que esperaba. El geotiff de la segunda idea parece estar desplazado 2 píxeles hacia el norte. Es por eso que consideré la idea 3 (o 4) como la correcta.
Pero si intenta con un mosaico en un nivel de zoom mucho más alto, el geotiff de la idea 3 se desplaza definitivamente hacia el sur. Fue una tontería comparar las coordenadas en un mosaico de nivel de zoom 3. Los límites del país en ese nivel de zoom deben simplificarse para que la comparación no dé buenos resultados.

Dan S. tenía razón, la imagen del mosaico ya está en EPSG: 3857. La segunda idea es la correcta (y también da buenos resultados a niveles de zoom altos)

Primera idea: EPSG: 4326
El código EPSG para las coordenadas WGS84 es EPSG: 4326 . Entonces, simplemente uso las coordenadas WGS84 para georreflexionar el mosaico 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

El mapa resultante se muestra en el lugar correcto, pero me temo que la proyección no es correcta y que podría haber un cambio en el centro del mosaico. Después de intentar mucho tiempo para comprobar eso al volver a proyectar el mapa con gdalwarp, descargué una versión de demostración de Global Mapper y este parece ser el caso (las fronteras de sames como en la idea 3 pero un cambio dentro del mosaico). La imagen debe extenderse para poder utilizar las coordenadas EPSG: 4326.

Segunda idea: EPSG: 3857
Este mosaico utiliza una proyección "web mercator" (alias google map projection), que ahora tiene un código EPSG: EPSG: 3857 (alias EPSG: 900913). Simplemente convierto las 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

Mis coordenadas en metros son:

0 10018754.1713946 5009377.08569731 5009377.08569731 (West North East South)

Ahora puedo usar gdal_translate para generar un 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

Mi impresión es que esto no es correcto porque los bordes de los mapas se desplazan hacia el norte. Parece ser la idea correcta.

Tercera idea: EPSG: 3857 a EPSG: 4055
Leí que "web mercator" usa coordenadas WGS84 pero las considero como si fueran coordenadas esféricas. Debido a la diferencia entre una latitud geodésica y una latitud geocéntrica (ver Wikipedia sobre la latitud ), los valores de latitud no serán los mismos en un elipsoide o en una esfera. Encontré que EPSG: 4055 es el código para coordenadas esféricas en una esfera basada en WGS84.

Convirtiendo las coordenadas a 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

Las coordenadas esféricas correspondientes son entonces:

0 66.3722684317026 45 40.7894557844857 (West North East South)

Luego hago como si esas coordenadas aún estuvieran en el elipsoide (EPSG: 4326) y las convierto 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

Las coordenadas resultantes difieren de la idea2:

0 9979483.26733298 5009377.08569731 4981335.86590183 (West North East South)

Ahora solo tengo que escribir las coordenadas en el 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

Esta tercera idea parece dar los mejores resultados. Pero no estoy seguro si es correcto. Si la idea 3 es correcta, ¿hay un código EPSG para hacer esta operación en un solo paso?

Cuarta idea: EPSG: 3857 a través de towgs84 = 0,0,0,0,0,0,0

gdal (y aparentemente epsg también) define EPSG: 3857 así:

+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

mientras que spatialreference.org así:

+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 uso la definición de spatialreference.org, obtengo las coordenadas correctas en un solo paso (bueno, todavía no lo hago si son las coordenadas "correctas", pero al menos son las mismas que en la idea 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 qué hay tanta diferencia en las definiciones de EPSG: 3857?

La imagen del mosaico ya está en EPSG: 3857. ¿Por qué no simplemente crear un archivo mundial para referenciarlo?

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

Para el mosaico que cubre América del Norte en el zoom 1, estaría viendo los siguientes contenidos de archivos mundiales:

78271.517
0
0
-78271.517
-19998372.6
19998372.6

De dónde provienen esos números:

• Línea 1: ancho de un píxel de imagen en coordenadas mundiales = 20037508.342789244 metros / 256 píxeles.
• Línea 2 y 3: rotación, entonces n / a.
• Línea 4: altura de un píxel de imagen en coordenadas mundiales. Igual que la línea 1 pero negativa, porque en los archivos de imagen, el aumento de y corresponde a 'abajo' mientras que en el sistema de coordenadas, el aumento de y corresponde a 'arriba'.
• Línea 5: coordenada X en coordenadas mundiales del centro del píxel superior izquierdo. Esto es -20037508.342789244, según lo informado por el enlace de mosaicos a la carta, más 1/2 de un píxel para llevarlo al centro.
• Línea 6: Ídem, solo coordenada Y de la esquina superior izquierda.

GDAL debería recoger su archivo mundial (.pgw para el png); aún tendrá que decirle EPSG: 3857 en la línea de comando.

(Nota: no tuve tiempo de probar esto, así que todo está fuera del alcance ... ¡pero espero que sea correcto en el primer intento de todos modos!)

Funciones necesarias para la generación de mosaicos globales utilizados en la web. Contiene clases que implementan conversiones de coordenadas para:

• GlobalMercator (basado en EPSG: 900913 = EPSG: 3785 ) para mosaicos compatibles con Google Maps, Yahoo Maps, Microsoft Bing Maps

• GlobalGeodetic (basado en EPSG: 4326) para OpenLayers Base Map y mosaicos compatibles con Google Earth

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

Acerca de mi pregunta secundaria en las cuarta ideas: ¿Por qué hay tanta diferencia en las definiciones de EPSG: 3857 entre la definición en gdal y en spaceialreference.org :

La principal diferencia es que gdal usa "+ nadgrids = @ null" y referencia espacial "+ towgs84 = 0,0,0,0,0,0,0". De acuerdo con la documentación o el formato PROJ.4, ambos parámetros se utilizan para las transformaciones de datos. Encontré un comentario interesante de Mikael Rittri en el servidor de listas 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."

El uso de "+ towgs84 = 0,0,0,0,0,0,0" no parece ser correcto o al menos solo bajo condiciones específicas.