GIS Daten und WGS Koordinaten

[Peche1985]

Hallo zusammen!

Habe hier im Forum schon einige interessante Diskussionen betreffend Karten- Projektionen / Interpolationen gefunden. Da ein Informatik- Kollege und ich aufgrund einer Diplomarbeit momentan mit diesem Thema beschäftigt sind, wollte ich euch um Rat fragen.
Wir arbeiten an einem Projekt, bei welchem wir ein Programm zur Flugrouten- Planung entwickeln. Hierfür dient eine ICAO- Karte, welche die Schweiz im Format 1:500'000 mit WGS Koordinaten darstellt. Zusätzlich ist es unsere Aufgabe ein kleines GIS zu entwickeln, so dass durch Eingabe eines Ortes ein Wegpunkt an entsprechender Koordinate in der Karte auf dem Bildschirm erzeugt wird.

Nun zu meinen 2 Fragen bezüglich dieser Ausgangslage.

1. GIS:
Einen ersten Prototypen dieses Programms haben wir bereits. Die Schweizer GIS- Daten (Ortsname und Longitude/Latitude- Koordinaten) haben wir von diesem Link http://ecotect.com/downloads/locationdata" target="_blank. Leider sind die Angaben etwas veraltet und zum Teil nicht ganz genau. Hat jemand eine Idee wo man aktuelle GIS- Daten (evt. kostenfrei) beziehen kann? Evt. auch mit Höhenangaben der einzelnen Orte..??

2. Umrechnung Pixel - WGS Koordinaten:
Seit einigen Wochen haben wir uns etwas in Themen wie geografische Projektion und Interpolation eingearbeitet. Die ganze Thematik ist jedoch umfassender als wir erwartet haben, und uns fehlt momentan ein entscheidender Ansatz. Evt. kannst du uns hierbei mit deinen Erfahrungen etwas helfen. Für Berechnungen ideal wäre ja eine zylindrisch projezierte Karte (z.B. Mercator Projektion). Mit der ICAO Karte haben wir aber eine Lamberts winkeltreue Kegelprojektion.
Durch eine Ortsangabe oder durch eine manuelle Eingabe der Ortskoordinaten (z.B. 8.267 und 47.083 für Luzern/CH) soll nun am entsprechenden Ort auf der ICAO Karte ein farbiger Punkt erzeugt werden. Hierfür ist aber eine vorgängige Kalibration nötig, damit das Programm weiss, welchem Pixel welche Koordinate zugeornet werden muss. Wir haben uns gedacht eine Kalibration von 4 vorgegebenen Punkten an den Ecken der Karte zu machen, bei welcher man für die entsprechende Koordinate in die Karte klickt. So dass das Programm weiss, welcher Pixelabstand zu welchem Koordinatenabstand gehört. Aber wir wissen nicht genau, wie wir von dem linearen X-/Y Pixelkoordinatensystem auf die WGS Koordinaten der angezeigten ICAO Karte kommen sollen. Hast du hier einen Lösungsansatz? Geht es evt. mit Einsatz einer zusätzlichen linearen Projektion (z.B. unser CH-1903)?

Ich wäre total froh, wenn ich diese Problemstellung in den Griff bekommen könnte. Für jeden Ratschlag bin ich also sehr dankbar!

Freundliche Grüsse,
Peter aus der Schweiz

[Christian Gründler]

Wir haben uns gedacht eine Kalibration von 4 vorgegebenen Punkten an den Ecken der Karte zu machen, bei welcher man für die entsprechende Koordinate in die Karte klickt. So dass das Programm weiss, welcher Pixelabstand zu welchem Koordinatenabstand gehört. Aber wir wissen nicht genau, wie wir von dem linearen X-/Y Pixelkoordinatensystem auf die WGS Koordinaten der angezeigten ICAO Karte kommen sollen.

Hallo Peter,

verstehe ich es richtig: Ihr habt eine Karte in Schnittkegelprojektion, die digital als Bitmap vorliegt, und wollt von Pixelkoorditaten auf Geodätische Koordinaten bzw. umgekehrt umrechnen? Steffen hat etwas ähnliches für Deutschland gemacht ( viewtopic.php?f=18&t=8450 ) und dort eine Genauigkeit von 5 km erreicht. Die Schweiz ist kleiner, also könnte mit dem gleichen Verfahren eine Abweichung von 1 km realistisch sein. Falls das reicht, geht man folgendermaßen vor:

a) Für 4 gut identifizierbare Punkte in den Ecken der Karte (Orte, Flugplätze, o.Ä.) sowohl die Geodätischen Koordinaten als auch die Pixelkoordinaten ermitteln. (Diese Punkte müssen nicht genau in der Ecke liegen und auch nicht ein Rechteck darstellen, wichtig ist ihre Identifizierbarkeit.)

b) Umrechnung für alle anderen Koordinaten durch bilineare Interpolation zwischen diesen vier Punkten; das geht in beide Richtungen.

Eine genauere Lösung werdet Ihr nur erreichen können, wenn Euch die Parameter der Kartenprojektion bekannt sind, d.h. Lage und Maßstab der beiden Standard-Breitenkreise in Pixelkoordinaten.

M.f.G. Christian

[Steffen81]

Richtig. Die geringe Genauigkeit von 1 bis 2 km (wenige Extremfälle bis ca. 5km) liegt an der mir nicht exakt bekannten Kartenprojektion des Kartenmaterials und dass dieses zum Teil wahrscheinlich nicht topografisch korrekt ist (demografische Karten).

GIS-Daten kann man bei GoogleMaps abfragen. Solange das nicht zu kommerziellen Zwecken genutzt wird, ist das auch kostenfrei mit (z.Z. glaube ich 15.000 Anfragen pro Tag und Nutzer). Dort lässt sich ein JSON-Objekt abrufen, das ne ganze Menge Informationen enthält. Klappt prima: GoogleMaps API
Haken dabei nur, man muss das was man sucht ungefähr kennen oder versucht es mit flächendeckendem Reverse Geocoding.

Edit: Orthographie.

[Peche1985]

Hallo zusammen!

Danke für die rasche und aufschlussreichen Antworten! Bin momentan mit der Google Maps API beschäftigt, um diese GIS- Daten zu bekommen.

@Christian,
du hast mein Problem genau erfasst.... Auch dein Vorgehensvorschlag erscheint mir logisch. Kannst du mir evt. Schritt b) noch etwas genauer erläutern? Ich habe bis jetzt noch nie eine bilineare Interpolation durchgeführt. Hast du hier ein gutes Beispiel/Ansatz wie ich das in unserem Fall angehen könnte? bzw. wo ich mir Informationen hierzu beschaffen könnte?

Freundliche Grüsse,
Peter

[Steffen81]

Das Verfahren der bilinearen Interpolation hatte er in diesem Beitrag erläutert.

Die Referenzpunkte sind dort für Deutschland gewählt. Sind euch die Standardparallelen der ICAO-Karte bekannt? Für eine deutsche ICAO-Karte, die ich hier gefunden habe, liegen diese auf dem 50. und 54. Breitengrad. Für die Schweiz dürfte(n) diese(r) natürlich woanders liegen, jenach dem ob es sich um eine einfache oder konforme Kegelprojektion handelt.

Es ist zwar redundant aber hier nocheinmal meine Methode, die Christians Berechnungsvorschlag umsetzt:

Code: Alles auswählen

public function projectByLongLat($L,$B)
{
	$L1	= $this->reference["nw"]["long"];
	$B1	= $this->reference["nw"]["lat"];
	$x1	= $this->reference["nw"]["x"];
	$y1	= $this->reference["nw"]["y"];
	
	$L2	= $this->reference["ne"]["long"];
	$B2	= $this->reference["ne"]["lat"];
	$x2	= $this->reference["ne"]["x"];
	$y2	= $this->reference["ne"]["y"];

	$L3	= $this->reference["sw"]["long"];
	$B3	= $this->reference["sw"]["lat"];
	$x3	= $this->reference["sw"]["x"];
	$y3	= $this->reference["sw"]["y"];
	
	$L4	= $this->reference["se"]["long"];
	$B4	= $this->reference["se"]["lat"];
	$x4	= $this->reference["se"]["x"];
	$y4	= $this->reference["se"]["y"];
	
	$B11	= $B1 + ($L - $L1) / ($L2 - $L1) * ($B2 - $B1);
	$x11	= $x1 + ($L - $L1) / ($L2 - $L1) * ($x2 - $x1);
	$y11	= $y1 + ($L - $L1) / ($L2 - $L1) * ($y2 - $y1);
	
	$B12	= $B3 + ($L - $L3) / ($L4 - $L3) * ($B4 - $B3);
	$x12	= $x3 + ($L - $L3) / ($L4 - $L3) * ($x4 - $x3);
	$y12	= $y3 + ($L - $L3) / ($L4 - $L3) * ($y4 - $y3);

	$x	= $x11 + ($B - $B11) / ($B12 - $B11) * ($x12 - $x11);
	$y	= $y11 + ($B - $B11) / ($B12 - $B11) * ($y12 - $y11);

	return array($x,$y);
}

[Christian Gründler]

Was mir noch eingefallen ist: möglicherweise ist auf der Karte ja das geodätische Koordinaten-Gitter eingetragen. Dann könnte man sich natürlich bei den Referenzpunkten an die Schnittpunkte der Gitterlinien halten. Aber das ist vermutlich zu schön, um wahr zu sein.

[Roland Ziegler]

Zur Info:

TransDEM arbeitet an dieser Stelle zweistufig, wenn man über 3+1 Punkte georeferenziert. Zunächst werden über die drei einzugebenden Punkte die Parameter einer gewöhnlichen affinen Transformation bestimmt, z.B. über Gitterschnittlinien, zu deren Pixel-Koordinaten die Koordinaten der Kartenprojektion eingegeben werden. Steht diese Transformation, wird die Pixel-Karte auf die Kartenprojektion rücktransformiert. Die Formeln der Kartenprojektion müssen bekannt sein (Erdöl-Dokument). Dies ermöglicht die inverse Projektion in geografische Koordinaten. Zum Schluss muss das geodätische Datum kompensiert werden. Das mach ich über die Wandlung in geozentrische Koordinaten und die Helmert-Transformation. Auch die Parameter für die Dati findet man bei bei den Erdöl-Leuten.

Nachzulesen hier unter "Georeferenzierung". (Der letzte Schritt in der zweiten Box - UTM-Projektion - fiele für die hier beschriebene Aufgabe weg.)

[Peche1985]

Hi all!!!

Merci für die guten Ratschläge!

Bin einer guten Lösung einen Schritt näher gekommen. Einerseits durch eure Beiträge und noch durch einen Tipp von jemandem, der bei der Swisstopo arbeitet (Geoinformationszentrum der Schweiz).

Momentan löse ich das also folgendermassen:

Benutzer öffnet das Programm und bindet seine ICAO Karte (mit beschrifteten Längen- und Breitengraden) ein. Nun erfolgt die Kalibration. Der Benutzer klickt 2 Koordinaten an, die sich vorzugsweise in der Mitte und in einer Ecke der Karte befinden und auf dem Gitternetz liegen. So dass der Benutzer zu diesen Punkten also in einem Dialog auch noch angeben kann, welche WGS84 Koordinaten es genau sind (Longitude und Latitude). Danach rechne ich mithilfe gegebener Formeln von Swisstopo die beiden Punkte vom WGS84- Format in das CH1903 - Format um (ist eine lineare Schweizerprojektion). Danach weiss ich, bei welcher X- und Y- Pixeldifferenz ich welche CH1903- Differenz habe.
Klickt der Anwender nach dieser Kalibration an einen beliebigen Ort in die Karte, berechene ich als erstes die Pixeldifferenz zum kalibrierten Referenzpunkt in der Mitte der Karte. Mit dem durch die Kalibration ermittelten Faktor "CH1903 Distanz pro Pixel Distanz" ermittle ich die entsprechenden CH1903 Koordinaten. Aus diesen wiederrum lassen sich abschliessend die WGS Koordinaten des Ortes ermitteln, welcher der Benutzer gedrückt hat. Et Voila.

Mit dieser Methode erreiche ich eine ziemlich gute Genauigkeit. Erst wenn ich am äusseren Rand der Karte einen Ort anwähle habe ich eine spürbare Abweichung von max. 2 bis 3 km.

Anhand dem Swisstopo- Mitarbeiter wäre die genaueste Methode jedoch, wenn ich direkt die Formeln der Lambert- Projektion programmieren würde. Dies ist aber wohl ziemlich kompliziert und ergäbe für die von uns verwendete Karte (1:500'000) vermutlich nur ein geringfügiges Verbesserungspotential.

Vermutlich werde ich den momentan erarbeiteten Lösungsansatz für dieses Projekt verwenden. Falls von euch aber noch Anregungen kommen, bin ich natürlich weiterhin offen und dankbar!

Freundliche Grüsse,
Peter

[Roland Ziegler]

Die Schweizer Landesprojektion ist eine schiefachsige Mercator-Projektion. Deine Karte aber ist eine Lambertsche Schnittkegelprojektion. Indem Du Deine Lambert-Karte als Swissgrid-Karte interpretierst, baust Du natürlich einen erheblichen systematischen und nicht-linearen Fehler ein. Den von Christian skizizerten Weg hielte ich da für sauberer. Ist natürlich immer ein Frage der erwünschten Genauigkeit. Wenn es beim Ungefähren bleiben soll, kann man sich vielleicht mit solch groben Näherungen behelfen. Wenn die Ansprüche steigen, wird man wohl um ein mehrstufiges Verfahren, wie ich es für TransDEM oder wie andere es für ähnliche Programme umgesetzt haben, nicht umhin kommen.

[Hortz]

Hallo zusammen! Sehr interessanter thread, leider funktioniert der link im Eingangspost nicht mehr, hat da jemand Ersatz? Danke schonmal
PNG
Gruß Manuel