Internet-Karten in ein GIS importieren

Die Grundlage für viele eigene Karten ist meist eine gute Hintergrundkarte. Beispielsweise erleichtert eine solche die Visualisierung eigener GPS-Tracks erheblich. Eine einfache Möglichkeit, an solche Karten für das eigene GIS zu gelangen sind die zahlreichen Webmapping-Dienste (z.B. Openstreetmap), die es inzwischen gibt. Die Übertragung solcher Kartendaten in ein eigenes GIS (z.b. QGis oder MapWindow) verlangt jedoch ein wenig Knowhow.

Karte mit Mobile Atlas Creator laden

Die Karten verschiedener Online-Kartendienste können mit dem Programm Mobile Atlas Creator heruntergeladen werden. Aus Urheberrechtsgründen wurde die Zahl der Kartenprovider allerdings reduziert. Openstreetmap, Google Maps, Yahoo Maps und Bing Maps wurden entfernt. Bei der Installation ist darauf zu achten, dass Java installiert ist. In Mobile Atlas Creator können dann u.a. die Datenquelle, der Zoomlevel und der Kartenausschnitt eingestellt werden. Für den eigentlichen Zweck des Programmes, Karten für die Mobiltelefon-Navigation mit dem sehr empfehlenswerten Programme TrekBuddy bereitzustellen, gibt es mehrere Atlasformate, in dem die Karten gespeichert werden können. Wenn die Karte an einem Stück gespeichert werden soll, muss jedoch als Atlas-Format "png+Worldfile" gewählt werden. Der Kartenausschnitt kann dann direkt von einem GIS-Programm wie QGis oder MapWindow geladen werden, ohne dass weitere Schritte nötig sind. Zu beachten ist lediglich, dass sich die Karte in der Projektion "Web Mercator" befindet. eingestellt werden. Frühere Versionen boten nur die Option "OziExplorer" an. Dabei wird ebenfalls eine große png-Datei erstellt und zusätzlich ein map-File, in dem die geographischen Koordinaten (Länge/Breite in Grad) der vier Eckpunkte des Bildes gespeichert sind. Diese kann ebenfalls so weiterverarbeitet werden, dass das png-File in ein GIS geladen werden kann, siehe nachfolgdende Absätze. Bei Verwendung von "png+Worldile" entfallen jedoch die folgenden Schritte.

Karte georeferenzieren

Das von Mobile Atlas Creator gespeicherte map-File kann benutzt werden, um die Karte in einem GIS-Programm zu georeferenzieren. Man muss dazu wissen, dass Google und Co. in ihren Webanwendungen zwei verschiedene Koordinatensysteme verwenden. Die von Mobile Atlas Creator verwendeten Koordinaten sind Längen/Breitengrad und beziehen sich auf den WGS84-Ellipsoid. Die Karten-Bilddateien sind hingegen in einer speziellen Mercator-Projektion erstellt. Diese wurde von ESRI "Web Mercator" genannt und unterscheidet sich von der normalen "World Mercator"-Projektion (EPSG:3395) dadurch, dass kein Ellipsoid sondern eine Kugel ("Sphere") als Form der Erde unterstellt wird (EPSG:3857). Um die Karte lagerichtig zu laden, gibt es dann mehrere Möglichkeiten, je nachdem welche Software zur Verfügung steht. Das optimale Ergebnis erhält man, wenn die Karten in dem Koordinatensystem georeferenziert werden, in dem sie projiziert sind, hier also in "Web Mercator". Für nicht ganz so große Gebiete liefert jedoch auch eine Georeferenzierung in "World Mercator" noch gute Ergebnisse. Beispielsweise beträgt die Abweichung einer Deutschlandkarte in der Mitte Deutschlands bis zu 120 m. An der oberen und unteren Kante stimmt die Lage auch bei World Mercator exakt, da dort ja die Referenzpunkte liegen.

Georeferenzieren mit Grad-Koordinaten mit ESRI ArcMap

Die Daten aus dem map-File können direkt als Referenzpunkte z.B. für ESRI ArcMap verwendet werden. Es müssen lediglich die Grad und Minutenangaben zu Dezimalgrad umgerechnet werden, Längen- und Breitengradangabe vertauscht werden, und die y-Wert-Angaben der Bildkoordinaten mit -1 multipliziert werden (Ozi misst von Nord nach Süd, während ESRI von Süd nach Nord misst). Dies kann durch folgendes Perl-Skript automatisiert werden.

map2georef.pl
Aufruf:
perl map2georef.pl infile.map [mercwgs|mercsphere]

Das dabei entstehende Textfile kann in ArcMap bei der Georeferenzierung direkt geladen werden. Ohne Projektion der Koordinaten bleibt man jedoch bei den geographischen Koordinaten, so dass das Projekt auf "WGS84" eingestellt sein muss. Dies ist aus o.g. Gründen nicht optimal.

Georeferenzierdatei mit projizierten Koordinaten

Durch Anhängen des Parameters "mercwgs" bzw. "mercsphere" projiziert das Perl-Skript map2georef.pl die Koordinaten in das World Mercator bzw. Web-Mercator-System. Vorzugsweise sollte Web-Mercator verwendet werden, sofern das eigene GIS dies untersützt. Damit kann die Karte dann in dem Koordinatensystem georeferenziert werden, in dem sie dargestellt ist. Falls das GIS-Programm wie ältere Versionen von Quantum-GIS keine Projektion von Raster-Layern unterstüzt, kann es sinnvoll sein, den Layer gleich in der Projektion World Mercator zu georeferenzieren. (Bei sehr kleinen Kartenausschnitten könnte wegen der Konformitätseigenschaft der Mercator-Projektion die Karte sogar direkt im UTM- oder Gauss-Krüger-System georeferenziert werden. Dazu müssten die Kartenkoordinaten im map-File entsprechend umgerechnet werden, dies ist im Skript nicht implementiert.)

Worldfile für das Koordinatensystem "WGS84" erstellen

Mit Hilfe eines "Worldfiles" kann die Grafikdatei direkt georeferenziert geladen werden. Worldfiles müssen den gleichen Dateinamen wie die Bilddatei mit einer speziellen Worldfile-Endung besitzen (.wld oder Bilddatei-spezifisch), damit sie von den GIS-Programmen erkannt werden. Ein solches Worldfiles kann aus dem map-File mit Hilfe des folgenden Perl-Skriptes erstellt werden.

map2world.pl
Aufruf:
perl map2world.pl infile.map [ext] [mercwgs|mercsphere|utm [zone]]

Dabei kann mit ext der Dateityp der zugehörigen Bilddatei angegeben werden. Dann erhält das Worldfile gleich die richtige Endung. Wird der Parameter nicht angegeben, wird .wld verwendet. Mit dem dritten Parameter wird gesteuert, ob die Koordinaten umprojiziert werden sollen. Wird kein dritter Parameter angegeben, wird das Worldfile bezogen auf die geographischen Koordinaten in Grad erstellt. Beim Laden im GIS-Programm muss die Projektion des Projektes dann auf WGS84 eingestellt sein. Wie oben erwähnt, ist diese Methode nicht optimal, da die Karten von Google etc. in der Projektion "Web Mercator" erstellt sind.
Das Skript verwendet die ersten drei im map-File gespeicherten Punkte um die Parameter des Worldfiles zu bestimmen. Diese definieren eine Abbildung der Bildkoordinaten auf das geographische Koordinatensystem. Der vierte oder weitere vorhandene Referenzpunkte werden zur Kontrolle der Abbildung verwendet, im Falle einer Abweichung wird eine Warnung ausgegeben.

Worldfile für die Projektion "World Mercator" oder "Web Mercator" erstellen

Um eine höhere Lagegenauigkeit zu erhalten, sollte die Karte direkt in der Projektion geladen werden, in der sie erstellt wurde, also bei den Karten aus dem Internet "Web Mercator". Dazu benötigt man ein World-File, dessen Werte sich ebenfalls auf dieses Koordinatensystem beziehen. Durch Anhängen des entsprechenden Parameters "mercsphere" konvertiert das Skript die Längen/Breitengradangaben aus dem Map-File in die entsprechend projizierten Koordinaten. Im GIS kann dann die Karte mit übereinstimmender Projektion der Bilddatei und des World-Files geladen werden. Falls das GIS-Programm die Projektion "Web Mercator" nicht untersützt oder die Karte dann nicht in eine andere Projektion transformieren kann, ist es praktischer, die Karte gleich in der "World Mercator"-Projektion zu laden. Natürlich muss die Projektion im World-File und die Projektionseinstellung im GIS übereinstimmen.
Das Perl-Skript map2world.pl kann die Koordinaten von Länge/Breite nach "World Mercator" und "Web Mercator" projizieren. Dazu muss der Parameter "mercwgs" bzw. "mercsphere" angefügt wird. Ältere Versionen von QuantumGIS (< 1.0) verfügen jedoch auch noch nicht über die World Mercator-Projektion. Diese kann jedoch manuell hinzugefügt werden. Die dafür notwendigen Parameter im proj4-Format findet man bei spatialreference.org.

Notlösung: Worldfile für die UTM-Projektion

Falls man nicht über die Möglichkeit verfügt, Rasterkarten umzuprojizieren (z.B. bei älteren Versionen von QGIS), kann man auch ein Map-File erstellen, das das Kartenbild direkt in der Projektion UTM georeferenziert. Dies kann natürlich nicht exakt sein, da sich eine Mercator- und eine UTM-projizierte Karte ja unterscheiden. Lediglich für kleine Kartenausschnitte von wenigen Kilometern ist dieses Verfahren wegen der Konformität der Mercatorprojektion akzeptabel. Der Lagefehler beträgt z.B. bei Karten von 7 km Ausdehnung ca. 10 m.
Zu diesem Zweck verfügt das Skript map2world.pl die Option "utm". Als weiterer Parameter kann die UTM-Zone angegeben werden, fehlt diese Angabe, wird die Zone automatisch bestimmt. Da die Abbildung des World-Files keine korrekte Umprojizierung der Karte wiedergeben kann, erscheint eine Warnung, in der der Fehler des vierten Punnktes abgelesen werden kann.

Eingescannte Karten für das Smartphone (Trekbuddy) aufbereiten

Um Karten aller Art, insbesondere solche auf Papier auch im Mobiltelefon nutzen zu können, sind einige Arbeitsschritte nötig.
1. Karte einscannen (möglicherweise abschnittsweise). Tipp für ESRI/ArcMap: Mit einem Bildbearbeitungsprogramm sicherstellen, dass drei Kanäle (und kein Alphakanal) vorhanden sind. Tipp für QGIS: Es darf ein Alpha-Kanal vorhanden sein, dies hat Vorteile (s.u.).
2. Karte georeferenzieren (in der Projektion, in der die Karte dargestellt ist, also z.B. UTM). Falls auf der Karte keine Referenzpunkte vorhanden sind, die der Projektion entsprechen (wenn beispielsweise geographische Koordinaten eingezeichnet sind, obwohl die Karte in der UTM-Projektion dargestellt ist), so müssen diese zuerst in die Projektion projiziert werden. Eine Möglichkeit dafür ist, die Koordinaten der anvisierten Referenzpunkte in eine Tabelle zu schreiben, diese als x/y-Datei im WGS-Rferenzsystem im GIS zu landen und dann das Shapefile in die gewünschte Projektion (z.B. UTM) zu projizieren. Danach fügt man die UTM-Koordinaten der Attributtabelle zu und kann die Werte wieder auslesen.
3. Evtl. Teilkarten zu einer zusammensetzen ("Mosaic"). Tip für ESRI/ArcMap: Wenn beim Georeferenzieren das ESRI-Grid-Format gewählt wird, wird bei einer schief stehenden Karte kein schwarzer Rand ergänzt. Tipp für QGIS: Das Zusammenfügen im Format GeoTIF funktioniert problemlos. Dabei wird der Alpha-Kanal genutzt, wenn die Option "-n 0" und "-a_nodata 0" angegeben wird. Die Teilbilder werden in der Reihenfolge kombiniert, in der sie im Befehl angegeben werden, also die letzten liegen oben.
4. Karte in die Projektion "Web Mercator" mit einem Shperoid als Datum projizieren (Entsprechend dem System von Google/Yahoo/Bing-Maps: EPSG 3857; Esri: WGS 1984 Web Mercator (Auxiliary Sphere); QGIS: WGS 84 / Pseudo Mercator (EPSG:3857)).
5. Evtl. Karte beschneiden (z.B. schwarze Ränder beschneiden, die beim umprojizieren entstanden sind). Sicherstellen, dass die Farbtiefe nur 8 Bit pro Kanal beträgt, evtl. durch Kopieren (ESRI: "Copy Raster") die Farbtiefe verringern. In QGis bei "Speichern-Unter" die Einstellung "wie dargestellt" wählen.
6. Evtl. Auflösung der Karte ändern. Eine Auflösung von 20m/Pixel entspricht ungefähr der Zoomstufe 13 bei Google-Maps (in mittleren geographischen Breiten). Eine Zoomstufe entspricht einer Änderung der Auflösung um den Faktor 2. Erfahrungsgemäß können Karten im Maßstab 1:50.000, die mit 300dpi eingescannt wurden, einmal um den Faktor 2 (z.B. von 10 auf 20 m/pixel) verkleinert werden, so dass immer noch alles lesbar ist und der Überblick besser wird. Bei Karten im Maßstab 1:25.000 sollte evtl. nur auf 10 m/pixel vergröbert werden, damit die Schrift noch lesbar bleibt, dies entspricht dann Zoomstufe 14.
Tipp: In QGis können Punkte 4-6 in einem Schritt durchgeführt werden.
7. World-File erstellen. ESRI/ArcMap: Tool in der Toolbox; QGIS: Kommandozeilentool: "listgeo -tfw Karte.tif"
8. Map-File erstellen. Dazu gibt es das folgende Perl-Skript:

world2map.pl

Es wandelt das World-File in ein Map-File um.

Aufruf:
perl world2map.pl infile.ext [mapsizeXxmapsizeY] [mercwgs|mercsphere]

Dabei ist Infile.ext das Rasterbild, dessen World-File in ein Map-File konvertiert werden soll. Der Dateiname des World-Files wird aus dem Namen des Bildes abgeleitet (.png -> .pgw; .tif -> .tfw etc.). Wird das World-File nicht unter diesem Namen gefunden, wird nach dem generischen Worldfile "infile.wld" gesucht. Danach kann die Größe des Rasterbildes in Pixels, getrennt durch "x" angegeben werden. Wird diese nicht angegeben, wird mit Hilfe des Programms Imagemagick versucht, die Größe des Bildes automatisch zu bestimmen. Der Parameter mercwgs oder mercshere gibt an, von welcher World-Mercator-Projektion die Koordinaten in geographische Koordinaten zurückprojiziert werden sollen. Hier sollte die Projektion angegeben werden, in die man die Karte projiziert hat (Schritt 4, also in der Regel "mercsphere"). Im Map-File werden also wieder geographische Koordinaten (Längen/Breitengrad) gespeichert. Dies ergibt wieder das "hybride" System, wie bei Trekbuddy und Google Maps üblich.

9. Die Perl-Version von TBCutter (siehe Punkt 10) mag keine Tif-Karten mit GeoTiff-Metadaten. Diese müssen zunächst entfernt werden, z.B. mit dem Befehl
convert InMap.tif OutMap.tif 2> /dev/null
Das Programm convert ist Teil von Imagemagick.
10. Danach kann die Rasterkarte mit TBCutter für Trekbuddy aufbereitet werden. TBCutter schneidet das Rasterbild in viele kleine Kacheln mit z.B. 512x512 Pixeln und erstellt eine passende Index-Datei dazu. Es empfiehlt sich, die Kacheln nicht allzu klein zu machen, um deren Zahl geringer zu halten (z.B. 1024x1024 Pixel).
Die gesamten Ausgabedateien von TBCutter können nun auf das Mobiltelefon kopiert werden und die Karte steht dem Programm Trekbuddy zur Verfügung.
11. Evtl. einen Atlas erstellen: Mehrere Karten können zu einem Atlas zusammengefasst werden, so dass man leichter zwischen den Karten wechseln kann. Am einfachsten kann dazu ein (ungetarter) bestehender Atlas hergenommen werden und die neue Map als Unterverzeichnis hinzufügen. Es gibt auch ein Programm um dies direkt zu bewerkstelligen: TBMapper.
12. Evtl. die Karte oder den atlas "taren", also mit tar packen. Gepackte Karten haben den Vorteil, dass die vielen Kartenkacheln in einer Datei stecken und von manchen Mobiltelefonen besser verwaltet werden können. Um Atlanten und Maps zu packen und auszupacken gibt es das Programm jTBtar.

Exkurs 1: Die Abbildungsgleichung von Worldfiles

In einem Worldfile stehen sechs Parameter in folgender Reihenfolge untereinander:
A
D
B
E
C
F
Diese definieren eine Abbildung, die durch folgende Gleichung definiert ist:

x ' y ' = A B D E x y + C F left [ matrix{x' ## y'} right ] = left [ matrix{A # B ## D # E} right ] times left [ matrix{x ## y} right ] + left [ matrix{C ## F} right ]

Dabei sind x und y die Koordinaten eines Punktes in der Bilddatei (Pixel) und x' und y' die Koordinaten auf der Karte.
Bei drei gegebenen Bild- und Kartenpunkten ist die Abbildung eindeutig bestimmt. Um aus den Punkten die Abbildungsparameter A bis F zu berechnen, muss folgendes lineare Gleichungssystem gelöst werden:

x 1 ' x 2 ' x 3 ' y 1 ' y 2 ' y 3 ' = x 0 y 0 1 0 0 0 x 1 y 1 1 0 0 0 x 2 y 2 1 0 0 0 0 0 0 x 0 y 0 1 0 0 0 x 1 y 1 1 0 0 0 x 2 y 2 1 A B C D E F left [ matrix{x_{1}' ## x_{2}' ## x_{3}' ## y_{1}' ## y_{2}' ## y_{3}' } right ] = left [ matrix{x_0 # y_0 # 1 # 0 # 0 # 0 ## x_1 # y_1 # 1 # 0 # 0 # 0 ## x_2 # y_2 # 1 # 0 # 0 # 0 ## 0 # 0 # 0 # x_0 # y_0 # 1 ## 0 # 0 # 0 # x_1 # y_1 # 1 ## 0 # 0 # 0 # x_2 # y_2 # 1 } right ] times left [ matrix{A ## B ## C ## D ## E ## F} right ]

Dies erfolgt im obigen Skript durch Invertierung der großen Matrix und Multiplikation mit dem linken Kartenvektor. (Anmerkung: Da die große Matrix reduzibel ist, wäre auch eine Formulierung mit zwei getrennten Gleichungssystem für A, B, C bzw. D, E, F möglich gewesen.) Diese Methode ist so allgemein, dass die Berechnung auch dann korrekt funktioniert, wenn die Punkte nicht auf den Ecken des Bildes liegen. Sogar Drehungen und Scherungen wären möglich, dies ist jedoch mit den map-Files von Mobile Atlas Creator nicht nötig.

Exkurs 2: Die Gleichungen für die Mercator-Projektion

Für die Projektion "World Mercator" mit Bezug auf den WGS84 Ellipsoid sind folgende Größen notwendig:

a = 6378137 (große Halbachse des WGS84-Ellipsoids in Meter)
b = 6356752,314 (kleine Halbachse in Meter)
f = a - b a f = {a - b} over a  (Abflachung)
e = 2 f - f 2 e = sqrt{2 cdot f-f^2}   (Exzentrizität)

Die Abbildung folgt dann folgenden Gleichungen, wobei λ den Längengrad und φ den Breitengrad (in Rad) bedeutet. Das Ergebnis besteht auf dem Ostwert E und dem Nordwert N (in Meter).

E = a λ E = a cdot %lambda
N = a ln ( tan ( π 4 + φ 2 ) [ 1 - e sin ( φ ) 1 + e sin ( φ ) ] e 2 ) N = a cdot ln(tan({%pi} over {4} + {%phi} over {2} )[{1-e cdot sin(%phi) } over {1 + e cdot sin(%phi) } ]^ {e over 2} )

"Web Mercator" unterscheidet sich davon, dass als kleine Halbachse (b) ebenfalls der Wert der großen Halbachse (a) verwendet wird. Die Gleichung vereinfacht sich dadurch erheblich (e=0).

Exkurs 3: Die Gleichungen für die UTM-Projektion

Die Gleichungen für die UTM-Projektion unterscheiden sich deshalb von der normalen Mercator-Projektion, da die Zentralmeridiane nicht mit dem Äquator zusammenfallen und deshalb "schräg" zur Abplattung des Ellipsoids liegen. Die angegebenen Formeln sind deshalb auch nicht exakt durch Winkelfunktionen etc. ausdrückbar, sondern müssen über Reihen approximiert werden. Die Berechnung gilt wie bei der Mercator-Projektion für den WGS84-Ellipsoid und erfolgt ebenfalls über mehrere Zwischenschritte. Die Formeln für a und e sind aus Exkurs 2 zu entnehmen. λ bedeutet wiederum den Längengrad und φ den Breitengrad (in Rad). λ0 ist der Längengrad des Zentralmeridians, als z.B. bei Zone 32: 9°, siehe Wikipedia.

e ' 2 = e 2 1 e 2 e'^2=e^{2} over{1-e^2}
C = e 2 ( cos φ ) 2 1 e 2 C={e^{2}(cos`%phi)^2} over {1-e^2}
ν = a 1 e 2 ( sin φ ) 2 %nu=a over sqrt{1 - e^2 (sin`%phi)^2}
T = ( tan φ ) 2 T=(tan`%phi)^{2}
A = ( λ λ 0 ) cos φ A=(%lambda - %lambda_0) cos`%phi
M = a [ ( 1 e 2 4 3 e 4 64 5 e 6 256 ... ) φ ( 3 e 2 8 + 3 e 4 32 + 45 e 6 1024 + ... ) sin 2 φ + ( 15 e 4 256 + 45 e 6 1024 + ... ) sin 4 φ ( 35 e 6 3072 + ... ) sin 6 φ + ... ] M=a [(1 - e^2 over 4 - {3 e^4} over 64 - {5 e^6} over 256 - ...)%phi - ({3 e^2} over 8 + {3 e^4} over 32 + {45 e^6} over 1024 + ...)sin 2%phi + ({15 e^4} over 256 + {45 e^6} over 1024 + ...)sin 4 %phi - ({35 e^6} over 3072 + ...)sin 6%phi + ...]
FE = 500000 (falsche Ostung in Meter)
FN = 0 (falsche Nordung auf der Nordhalbkugel in Meter; auf der Südhalbkugel ist FN = 10000000)
k0 = 0,9996 (Skalierung)

Damit kann dann der Rechtswert und der Hochwert berechnet werden:

E = FE + k 0 ν [ A + ( 1 T + C ) A 3 6 + ( 5 18T + T 2 + 72C 58 e ' 2 ) A 5 120 ] E=FE+k_0%nu[A+ {(1 - T + C)A^3} over 6 + {(5-18T+T^2+72C-58e'^2)A^5} over 120]
N = FN + k 0 [ M ν tan φ ( A 2 2 + ( 5 5 T + 9C + 4 C 2 ) A 4 24 + ( 61 58T + T 2 + 600C 330 e ' 2 ) A 6 720 ) ] N=FN+k_0[M - %nu cdot tan %phi(A^2 over 2 + {(5 -5 T + 9C + 4 C^2)A^4} over 24 + {(61 - 58T + T^2 + 600C - 330e'^2)A^6} over 720)]

Es gibt für beide Projektionen auch Gleichungen für die Rückwärts-Projektion. Diese sind in den angegebenen Quellen dokumentiert.

Quellen:

International Association of Oil & Gas Producers (Ed.) (2009): Coordinate Conversions and Transformations including Formulas, Surveying and Positioning Guidance Note Number 7, part 2, Internet: http://www.epsg.org/guides/

Snyder, J.P. (1987): Map Projections - A Working Manual, U.S. Geological Survey Professional Paper 1395, Washington DC, USA.

Hauptseite

Dieses Dokument darf unter Nennung des Autors unter gleichen Bedingungen weitergegeben werden (Creative Commons CC-BY-SA) Creative Commons Lizenzvertrag
Joachim Aurbacher (joalbach@gmx.net) – 2011-2021

VG-Wort Zählpixel