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

16

Como ejemplo, tomaré el siguiente mosaico enlace y lo guardaré como "4_2.png" .

Las coordenadas WGS84 de este mosaico se pueden calcular o leer allí haciendo clic en la ficha 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) de modo 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 (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 GIS 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 lo intentas con un mosaico a 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 de los países en dicho nivel de zoom se simplifican 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 un buen resultado en niveles de zoom altos)

Primera idea: EPSG: 4326
El código EPSG para las coordenadas WGS84 es EPSG: 4326 . Así que simplemente utilizo las coordenadas WGS84 para georreferenciar 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 de la baldosa. Después de intentar comprobarlo durante mucho tiempo reproyectando el mapa con gdalwarp, descargué una versión demo de Global Mapper y este parece ser el caso (sames bordes como por idea 3 pero un cambio dentro de la baldosa). La imagen debe estar estirada para poder utilizar las coordenadas EPSG: 4326.

Segunda idea: EPSG: 3857
Este mosaico usa 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 través de 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 geocéntrica (consulte Wikipedia sobre la latitud ), los valores de latitud no serán los mismos en un elipsoide o en una esfera. Descubrí que EPSG: 4055 es el código para las coordenadas esféricas en una esfera basada en WGS84.

Convertir 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 estuvieran aún en el elipsoide (EPSG: 4326) y las convirtiera a mercator web:

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 de 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 de si es correcto. Si la idea 3 es correcta, ¿hay un código EPSG para realizar 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 también epsg) 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 [email protected] +wktext +no_defs

mientras que spatialreference.org es 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, obtuve las coordenadas correctas en un solo paso (bueno, todavía no las tengo si son las coordenadas "correctas", pero al menos son las mismas como por 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 tal diferencia en las definiciones de EPSG: 3857?

    
pregunta Name 25.10.2010 - 23:56

3 respuestas

11

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

enlace

Para el mosaico que cubre N. América en el zoom 1, verás los siguientes contenidos de archivo mundial:

78271.517
0
0
-78271.517
-19998372.6
19998372.6

De dónde vinieron esos números:

  • Línea 1: ancho de un píxel de imagen en coordenadas mundiales = 20037508.342789244 metros / 256 píxeles.
  • Líneas 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 negativo, porque en los archivos de imagen y y corresponde a "abajo" mientras que en el sistema de coordenadas, y 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 a la carta de las baldosas, más la mitad de un píxel para llevarlo al centro.
  • Línea 6: Ídem, solo coordenada Y de la parte superior izquierda.

GDAL debería recoger tu archivo de mundo (.pgw para el png); todavía tendrás que indicarlo EPSG: 3857 en la línea de comandos.

(Nota: no tuve tiempo de probar esto, por lo que todo está descartado ... ¡pero espero que sea correcto en el primer intento de todos modos!;)

    
respondido por el Dan S. 30.11.2010 - 21:35
0

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 los mapas de Google Maps, Yahoo Maps, Microsoft Bing Maps compatibles

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

enlace

    
respondido por el Mapperz 26.10.2010 - 03:24
0

Acerca de mi pregunta secundaria en la cuarta idea: ¿Por qué existe tal diferencia en las definiciones de EPSG: 3857 entre la definición en gdal y en spatialreference.org :

La principal diferencia es que gdal utiliza "+ nadgrids = @ null" y spatialreference "+ towgs84 = 0,0,0,0,0,0,0". Según la documentación o el formato PROJ.4 , ambos parámetros se utilizan para las transformaciones de referencia. 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 

   [email protected] 

instead."

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

    
respondido por el Name 19.09.2011 - 14:30

Lea otras preguntas en las etiquetas