Wie kann ich eine Web-Mercator-Kachel mit gdal korrekt georeferenzieren?

16

Als Beispiel nehme ich die folgende Kachel http://a.tile.openstreetmap.org/3/4/2.png und speichere sie als "4_2.png".

Die WGS84-Koordinaten dieser Kachel können berechnet oder gelesen werden, indem Sie auf die entsprechende Kachel klicken:

0 66.51326044311185 45 40.97989806962013 (West North East South)

So wird die Kachel korrekt georeferenziert (mithilfe von gdal wird ein Geotiff oder ein anderes georeferenziertes Format erstellt), sodass:

  • Die Bitmap muss nicht gestreckt werden (= die Pixel im Geotiff sind genau die gleichen wie in der ursprünglichen Bitmap)
  • Das resultierende Bild wird an der richtigen Stelle in einem GIS-Viewer / Editor geöffnet (wie zum Beispiel in TatukGIS Free Viewer ).

(Bearbeitet am 19. September 2011, um meine Frage klarer zu machen und meine Schlussfolgerungen zu berücksichtigen.)

Meine Schlussfolgerung:

Ich dachte zuerst, dass die dritte Idee (siehe unten) die richtige war. Ich habe den Geotiff in einem GIS-Viewer geöffnet und die angezeigten Koordinaten mit den erwarteten verglichen. Der Geotiff aus der zweiten Idee scheint um 2 Pixel nach Norden verschoben zu sein. Deshalb hielt ich Idee 3 (oder 4) für die richtige.
Wenn Sie jedoch eine Kachel mit einer viel höheren Zoomstufe verwenden, wird der Geotiff aus Idee 3 definitiv nach Süden verschoben. Es war albern, Koordinaten auf einer Kachel der Zoomstufe 3 zu vergleichen. Die Ländergrenzen bei dieser Zoomstufe sind vereinfacht, so dass der Vergleich keine guten Ergebnisse liefert.

Dan S. hatte recht, das Kachelbild ist bereits in EPSG: 3857. Die zweite Idee ist dann die richtige (und liefert auch bei hohen Zoomstufen gute Ergebnisse)

Erste Idee: EPSG: 4326
Der EPSG-Code für WGS84-Koordinaten lautet EPSG: 4326 . Daher benutze ich einfach die WGS84-Koordinaten, um die Kachel mit gdal_translate als Geotif zu georeferenzieren :

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

Die resultierende Karte wird an der richtigen Stelle angezeigt, aber ich befürchte, dass die Projektion nicht korrekt ist und es zu einer Verschiebung in der Mitte der Kachel kommen könnte. Nachdem ich lange versucht hatte, dies zu überprüfen, indem ich die Karte mit gdalwarp neu projizierte, lud ich eine Demoversion von Global Mapper herunter. Dies scheint der Fall zu sein (gleiche Rahmen wie in Idee 3, aber eine Verschiebung innerhalb der Kachel). Das Bild sollte gestreckt sein, um die EPSG: 4326-Koordinaten verwenden zu können.

Zweite Idee: EPSG: 3857
Diese Kachel verwendet eine "Web Mercator" -Projektion (alias Google Map-Projektion), die jetzt einen EPSG-Code hat: EPSG: 3857 (alias EPSG: 900913). Ich konvertiere einfach die Koordinaten mit 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

Meine Koordinaten in Metern sind: 

0 10018754.1713946 5009377.08569731 5009377.08569731 (West North East South)

Jetzt kann ich mit gdal_translate einen Geotiff erzeugen:

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

Mein Eindruck ist, dass dies nicht korrekt ist, da die Grenzen der Karten nach Norden verschoben sind. Es scheint die richtige Idee zu sein.

Dritte Idee: EPSG: 3857 bis EPSG: 4055
Ich habe gelesen, dass "web mercator" WGS84-Koordinaten verwendet, sie aber so betrachtet, als ob sie sphärische Koordinaten wären. Aufgrund des Unterschieds zwischen einem geodätischen und einem geozentrischen Breitengrad (siehe Wikipedia zum Breitengrad ) sind die Breitengradwerte auf einem Ellipsoid oder einer Kugel nicht gleich. Ich fand, dass EPSG: 4055 der Code für Kugelkoordinaten auf einer Kugel ist, die auf WGS84 basiert.

Umrechnung der Koordinaten in 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

Die entsprechenden Kugelkoordinaten sind dann:

0 66.3722684317026 45 40.7894557844857 (West North East South)

Dann mache ich so, als ob diese Koordinaten noch auf dem Ellipsoid wären (EPSG: 4326) und konvertiere sie in 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

Die resultierenden Koordinaten unterscheiden sich von denen von idea2:

0 9979483.26733298 5009377.08569731 4981335.86590183 (West North East South)

Jetzt muss ich nur noch die Koordinaten in die Karte schreiben:

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

Diese dritte Idee scheint die besten Ergebnisse zu liefern. Aber ich bin nicht sicher, ob es richtig ist. Wenn Idee 3 korrekt ist, gibt es einen EPSG-Code, um diesen Vorgang in einem Schritt durchzuführen?

Vierte Idee: EPSG: 3857 bis towgs84 = 0,0,0,0,0,0,0

gdal (und anscheinend auch epsg) definiert EPSG: 3857 so:

+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

während spatialreference.org so aussieht:

+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

Wenn ich die Definition von Spatialreference.org verwende, habe ich in einem Schritt die richtigen Koordinaten erhalten.

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

Warum unterscheiden sich die Definitionen von EPSG: 3857 so stark?

web-mapping gdal epsg tiles Name
quelle

Antworten:

11

Das Kachelbild ist bereits in EPSG: 3857. Warum nicht einfach eine Weltdatei erstellen, um darauf zu verweisen?

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

Für die Kachel, die Nordamerika bei Zoom 1 abdeckt, sehen Sie sich den folgenden Inhalt der Weltdatei an:

78271.517
0
0
-78271.517
-19998372.6
19998372.6

Woher diese Zahlen kamen:

  • Zeile 1: Breite eines Bildpixels in Weltkoordinaten = 20037508.342789244 Meter / 256 Pixel.
  • Zeile 2 und 3: Rotation, also n / a.
  • Zeile 4: Höhe eines Bildpixels in Weltkoordinaten. Wie Zeile 1, jedoch negativ, da in Bilddateien ein ansteigendes y 'abwärts' bedeutet, während im Koordinatensystem ein ansteigendes y 'aufwärts' bedeutet.
  • Zeile 5: X-Koordinate in Weltkoordinaten der Mitte des linken oberen Pixels. Dies ist -20037508.342789244, wie von den Kacheln A-la-carte-Link gemeldet, plus 1/2 Pixel, um es in die Mitte zu bringen.
  • Zeile 6: Dito, nur Y-Koordinate von oben links.

GDAL sollte Ihre Worldfile (.pgw für die PNG) abholen; In der Befehlszeile muss noch EPSG: 3857 angegeben werden.

(Hinweis: Ich hatte keine Zeit, dies zu testen, also ist alles an der Manschette ... aber hoffentlich gleich beim ersten Versuch richtig!;)

Dan S.
quelle
Vielen Dank und Entschuldigung für die späte Antwort. Eigentlich würde das aber zu meiner zweiten Idee führen, bei der gdaltransform wirklich nur zur Georeferenzierung des Bildes verwendet wird.
Nennen Sie
Sie hatten Recht mit dem Kachelbild, das bereits in EPSG: 3857 enthalten ist. Die Lösung bestand darin, einfach EPSG: 3857 zu verwenden. Auf diese Weise habe ich die falschen Ergebnisse überprüft.
Name
0

Funktionen, die für die Generierung globaler Kacheln im Web erforderlich sind. Es enthält Klassen, die Koordinatenumrechnungen für Folgendes implementieren:

  • GlobalMercator (basierend auf EPSG: 900913 = EPSG: 3785 ) für Google Maps- , Yahoo Maps- und Microsoft Bing Maps-kompatible Kacheln

  • GlobalGeodetic (basierend auf EPSG: 4326) für OpenLayers Base Map- und Google Earth-kompatible Kacheln

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

Mapperz
quelle
Danke, aber es beantwortet meine Frage nicht. Ich möchte keine Kacheln generieren. Ich möchte wissen, was die richtige Georeferenzierung von Fliesen ist / wäre.
Name
Ich habe "Wenn es richtig ist, gibt es einen EPSG-Code, um diesen Vorgang in einem Schritt durchzuführen?" Antwort EPSG: 900913
Mapperz
Sie haben Folgendes nicht getan: EPSG: 900913 funktioniert wie EPSG: 3857 (oder wie EPSG: 3785), und es ist nicht zulässig, die Operation in einem Schritt auszuführen. Sie müssen EPSG: 4055 noch durchgehen. Außerdem hatte ich in meiner ursprünglichen Frage bereits EPSG: 900913 als Alias ​​für EPSG: 3857 erwähnt.
Name
0

Zu meiner zweiten Frage in den vierten Überlegungen: Warum gibt es einen solchen Unterschied in den Definitionen von EPSG: 3857 zwischen der Definition in gdal und auf spatialreference.org :

Der Hauptunterschied besteht darin, dass gdal "+ nadgrids = @ null" und den Raumbezug "+ towgs84 = 0,0,0,0,0,0,0" verwendet. Entsprechend der Dokumentation oder dem PROJ.4-Format werden beide Parameter für Datumstransformationen verwendet. Ich habe einen interessanten Kommentar von Mikael Rittri auf dem Proj4-Listenserver gefunden :

"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."

Die Verwendung von "+ towgs84 = 0,0,0,0,0,0,0" scheint nicht richtig zu sein oder zumindest nur unter bestimmten Bedingungen.

Name
quelle