gis web-mapping – ¿Cómo georreferenciar un mosaico web mercator correctamente usando gdal?

Pregunta:

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 GIS (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 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. Por eso consideré la idea 3 (o 4) como la correcta.
Pero si intenta con un mosaico con 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 se han simplificado 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 entonces la correcta (y también da un buen resultado con niveles de zoom altos)


Primera idea: EPSG: 4326
El código EPSG para las coordenadas WGS84 es EPSG: 4326 . Así que simplemente uso las coordenadas WGS84 para georrefenciar 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 medio del mosaico. Después de intentar comprobarlo durante mucho tiempo volviendo a proyectar el mapa con gdalwarp, descargué una versión de demostración de Global Mapper y este parece ser el caso (los mismos bordes que en la idea 3, pero un cambio dentro del mosaico). La imagen se debe estirar para poder usar las coordenadas EPSG: 4326.


Segunda idea: EPSG: 3857
Este mosaico utiliza una proyección "web mercator" (alias proyección de mapa de Google), 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 están desplazados hacia el norte. Parece ser la idea correcta.


Tercera idea: EPSG: 3857 a EPSG: 4055
Leí que "web mercator" usa coordenadas WGS84 pero considérelas como si fueran coordenadas esféricas. Debido a la diferencia entre una latitud geodésica y geocéntrica (ver 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 coordenadas esféricas en una esfera basada en WGS84.

Conversión de 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 todavía en el elipsoide (EPSG: 4326) y las convierto 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, ¿existe 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 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 espacialreference.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 utilizo la definición de Spatialreference.org, obtuve las coordenadas correctas en un paso (bueno, todavía no las tengo 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?

Respuesta:

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

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

Para el mosaico que cubre Norteamérica con zoom 1, estaría mirando el siguiente contenido del archivo mundial:

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 negativo, porque en los archivos de imagen aumentar y corresponde a 'abajo' mientras que en el sistema de coordenadas, aumentar 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 píxel para llevarlo al centro.
  • Línea 6: Ídem, solo coordenada Y de la parte superior izquierda.

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

(Nota: no tuve tiempo de probar esto, así que todo está improvisado … ¡pero es de esperar que sea correcto en el primer intento de todos modos!;)

Leave a Comment

Your email address will not be published.

Scroll to Top

istanbul avukat

-

web tasarım