Neuere Einträge
09 Jun 2012
Coding PHP Mathematik 5

Umkreissuche

Vor einem knappen Monat ist etwas mir seit eh und je prophezeites passiert: Ich konnte im Studium gelerntes Mathewissen doch tatsächlich außerhalb der Veranstaltung „Mathematik für Ingenieure 1“ praktisch einsetzen! Aber der Reihe nach …

Zu Beginn ein Sicherheitshinweis: Das hier ist nichts für Leute, die Mathe für ein Arschl*ch halten und in dieser Ansicht auch weiterhin nicht beirrt werden möchten! ;-) Wer es nachbauen möchte, sollte Grundlagen von PHP und Datenbanken mitbringen.

Alles begann damit, dass mir die Aufgabe übertragen wurde, doch eine postleitzahlbasierte Umkreissuche für eine Datenbank zu bauen. Und hier möchte ich nun einen Weg demonstrieren wie man so etwas umsetzen kann. Ob es nun der günstigste ist, weiß ich nicht. Aber es ist mir ausreichend flott.

Geographischer Ort

Die erste Aufgabe ist es, eine Adresse bestehend aus Straße, Postleitzahl und Ort in den geographischen Ort umzuwandeln. Hierzu gibt es mehrere Wege nach Rom:
Zum einen bietet Google eine API an, mit der man vollständige Adressen in ziemlich exakte Orte auf unserer Erde umwandeln kann. Sie ist gut dokumentiert und liefert maschinenlesbare, sehr ausführliche Ergebnisse. Jedoch gilt es hierbei zu bedenken, dass man sich bezüglich der Nutzungseinschränkungen tief im Graugebiet befindet: Nutzung nur in Zusammenhang mit einer Anzeige auf Google Maps Karten; keine Nutzung bei Echtzeiteingaben durch den Benutzer; nur 2.500 Anfragen in 24 Stunden. In jedem Fall raten sowohl Google als auch ich zu einem lokalen Caching der Ergebnisse!
Möglichkeit Nummer zwei ist die OpenGeoDB, die aber nicht straßengenau, sondern nur postleitzahlgenau auflöst. Jedoch kann man die Datenbank selbst hosten und für Umkreissuchen legal nutzen.

Welche Methode man nun nutzt, muss jeder für sich selber entscheiden. Von hier an gehe ich davon aus, dass für alle Adressen und den Mittelpunkt der Umkreissuche Längengrad und Breitengrad vorliegen. Die folgenden Codebeispiele werde ich in PHP und SQL abhalten - dürften aber sehr leicht auf andere Sprachen portierbar sein.

Koordinatensystem

Die geographischen Koordinaten geben jeden Ort über ein sogenanntes Kugelkoordinatensystem an. (Zack! Erstes Buzzword aus Mathe 1.) Beispielsweise Hannover-Mitte liegt 52,4°N 9,8°O. In unserem Kugelkoordinatensystem bedeutet das, dass zwischen Hannover und dem Äquator ein Winkel von 52,4° und zwischen dem Nullmeridian und Hannover ein Winkel von 9,8° liegt. Ebenso ist Hannover 6371km vom Erdmittelpunkt entfernt. Die letzte Information lässt man aber bei den geographischen Koordinaten weg, da man davon ausgehen kann, dass sie sich immer auf Orte auf der Kugeloberfläche beziehen.

Für den weiteren Verlauf müssen wir die Ortsangaben vom Kugelkoordinatensystem in ein kartesisches Koordinatensystem überführen.

Das kartesische Koordinatensystem gibt Orte etwas anders an: Während im Kugelkoordinatensystem zwei Winkel und eine Distanz eine Position beschreiben, sind es im kartesischem System drei Distanzen. Hannover-Mitte kann man nämlich auch so erreichen: Wenn man beim Erdkern losgeht ( Vom praktischen Ausprobieren möchte ich an dieser Stelle ausdrücklich abraten! ;-) ), geht man zunächst 3833 km Richtung Golf von Guinea, dann 657 km Richtung Singapur und anschließend 5045 km Richtung Nordpol. Wenn man sich nicht verlaufen hat, stößt man schließlich in Hannover-Mitte durch die Erdkruste.

Mathematisch geschieht diese Überführung durch drei kurze Multiplikationen. Ein kleines Schaubild soll das verdeutlichen:

Programmtechnisch kommt folgendes dabei heraus:

  1. <?php
  2.         // $lat: Breitengrad
  3.         // $lng: Längengrad
  4.         // $r: Erdradius
  5.  
  6.         $a = array(
  7.                 ‘x’ => $r * cos( $lat ) * cos( $lng ),
  8.                 ‘y’ => $r * cos( $lat ) * sin( $lng ),
  9.                 ‘z’ => $r * sin( $lat )
  10.         );
  11. ?>

Nun hat man alle Adressen und den Mittelpunkt der Umkreissuche im kartesischem Format vorliegen. Für den weiteren Vorlauf nenne ich die Koordinaten Vektoren. Sie sind Zeiger, die vom Mittelpunkt der Erde auf den jeweiligen Ort zeigen.

Abstand zweier Orte

Letztlich ist nun die kürzeste Verbindung zweier „Vektorspitzen“ auf einer Kugel gesucht. Der geneigte Mathematiker spricht von Orthodromen. Dazu bestimmen wir zunächst den Winkel, der zwischen den Vektoren liegt. Ein kurzer Blick in mein Mathebuch verrät, dass uns das Skalarprodukt dabei gerne weiterhilft. Der Winkel sagt uns direkt, zu welchem Anteil wir die Erde umkreisen müssen. Beispielsweise würde ein Winkel von 90° bedeuten, dass wir ein Viertel des Erdumkreises ablaufen müssen, um von A nach B zu kommen.

Die Mathematik dahinter sieht dann so aus:

Der dazugehörige Code:

  1. <?php
  2.         // $a: Array mit Koordinaten von Vektor A
  3.         // $b: Array mit Koordinaten von Vektor B
  4.  
  5.         function dproduct( $a, $b ) {
  6.                 return $a['x'] * $b['x'] + $a['y'] * $b['y'] + $a['z'] * $b['z'];
  7.         }
  8.  
  9.         $angle = acos(
  10.                        dproduct( $a, $b ) /
  11.                        ( sqrt( dproduct( $a, $a ) ) * sqrt( dproduct( $b, $b ) ) )
  12.                      );
  13.  
  14.         // Da der Winkel im Bogenmaß vorliegt, reicht eine einfache Multiplikation
  15.         $dist = $angle * 6371;
  16. ?>

Wie gut das funktioniert, zeigt ein Proof of Concept mit der Abstandsbestimmung von Hannover und Gwynneville in Australien: Mit der oben vorgestellten Methode habe ich eine Abweichung von 5,3km gegenüber von Google Earth. Hierfür habe ich übrigens die Google Maps API angezapft.

Umkreissuche

Jetzt endlich kommen wir zur eigentlichen Umkreissuche. Ich fasse noch einmal zusammen, was wir gegeben haben:

  • Eine Datenbank mit allen Adressen und den vorher schon ausgerechneten kartesischen Koordinaten x, y und z.
  • Die kartesischen Koordinaten des Mittelpunkts der Umkreissuche
  • Den Radius der Umkreissuche

Ziel ist es nun eine passende Anfrage an den Datenbankserver zu generieren, damit dieser uns alle Datensätze zurück liefert, die im Umkreis enthalten sind. Hier der Code, der die Anfrage generiert:

  1. <?php
  2.         $query = sprintf( "SELECT * FROM adressen "
  3.                          ."WHERE %s / 6371 >= ACOS("
  4.                          ." ( %s * x + %s * y + %s * z ) /"
  5.                          ." ( %s * SQRT( x*x + y*y + z*z )"
  6.                          .")",
  7.                          $radius,
  8.                          $a['x'],
  9.                          $a['y'],
  10.                          $a['z'],
  11.                          sqrt( $a['x'] * $a['x'] +
  12.                                $a['y'] * $a['y'] +
  13.                                $a['z'] * $a['z'] ) );
  14. ?>

Hier wird auch deutlich, warum ich die Vorgehensweise, die unter dem Artikel Orthodrome gezeigt wird einfach so ignoriert habe: Bei dieser hier vorgestellten Methode muss pro Datensatz nur einen Arkuskosinus (ACOS) und eine Hand voll Multiplikation und Additionen berechnet werden. Bei der Wikipedia-Methode hingegen kommen noch einige Sinus- und Kosinusberechnungen hinzu, was etwas aufwändiger sein dürfte.

Ich hoffe, hiermit jemanden weitergeholfen zu haben. Ebenso hoffe ich, dass ich nun mit dem Text keine Mathematiker wegen der wohl recht unscharfen Ausdrucksweise zur Verzweiflung / Weißglut gebracht habe. Wenn doch, bitte ich um Pardon. ;-) Wenn jemand noch Verbesserungsvorschläge hat oder Fehler findet, darf dieser sehr gerne einen Kommentar schreiben!

Veröffentlich von Juergen Fitschen am 09.06.2012
Permanentlink: http://bequemlichkeitsgruen.de/post/1339252991

5 Kommentare wurden hinterlassen

#110.06.2012 15:46

Tobby:

Also ich habs verstanden! :-D

Ich frage mich allerdings, wie man auf die Idee kommt, „Gwynneville“ in Australien als Ziel für die Berechnung zu benutzen. Wer will denn schon nach Gwynneville??? ;)

Schöne Grüße aus Gwynneville!

#210.06.2012 16:12

Jue:

Alles natürlich reiner Zufall. Genau so ist Hannover natürlich völlig aus der Luft gegriffen. Timbuktu wäre ja sogar naheliegender.

#313.06.2012 15:31

Oli:

Erstmal muss ich sagen klasse und verständlicher Artikel ;)

Eine Anmerkung habe ich noch. Da du ja schon extra Vektoren benutzt kannst du durch das Punktprodukt und der Norm (Länge) mit ziemlich geringen Kosten[1] - trotzdem winkelbasierend - rechnen. In der Grafikprogrammierung wird das ebenfalls gerne verwendet. Hier ein Beispiel (C++):

  1. float Class::cosOfAngle(const Vector3f& a, const Vector3f& b) const
  2. {
  3.     /* it returns scalarproduct (a*b) diveded by the lengths.
  4.        such that cos_of_angle(a,b) = cos <| (a, b)
  5.        -> <- = -1
  6.        -> -> = 1
  7.        |^ -> = 0
  8.        <- |^ = 0 */
  9.     return a.dot(b) / (a.norm() * b.norm());
  10. }

D.h. man rechnet einfach mit dem cos(winkel) zwischen 2 Vektoren weiter und setzt nicht nochmal den arccos drauf. Du musst das hier aber anscheinend machen, weil du bei der Distanz nicht einfach sagen kannst:

  1. {
  2. return (a-b).norm();
  3. }

Sondern die Distanz bei einem Weg über die Kugeloberfläche wissen möchtest. Die Frage also ob man diese nicht einfach anders berechnet und sich somit den arccos spart (der ist ja sowieso nur dazu da um den cos(WinkelZwischenAundB) wieder in den WinkelZwischenAundB zu basteln).

Grüße,

Oli

[1] wenige Aufrufe von + und * kostet für den Rechner nix.Trigonometrische Funktionen hingegen schon (wobei außerdem hohe Ungenauigkeit reinkommt. Es kann z.B. passieren, dass du Probleme bekommst bei fast 0° bzw. fast 360° Winkel zwischen Ort A und Ort B)

#413.06.2012 16:50

Jue:

Danke für das Kompliment und deine Anmerkungen, Oli! :)

Ich habe mir mal erlaubt, Syntaxhighlighting um deinen Code zu bauen - sonst geht das schöne Einrücken und vor allem die ASCII-Illustration kaputt.

Ja, wenn ich direkt (a-b).norm() berechne, mag das für kurze Distanzen sicherlich genau genug sein. Aber man muss bedenken, dass man dann durch die Erde hindurch geht und sich nicht auf deren Oberfläche bewegt.

Ich hatte überlegt, ob ich nicht einfach die ACOS auf die andere Seite der Ungleichung ziehe, da ja nur der rechte Teil für jeden Datensatz neu berechnet werden muss. Der linke Teil wird nur einmal für einen Query berechnet. (Zumindest gehe ich davon aus, dass mySQL da so intelligent ist.) Aber ich bin mir damals nicht ganz sicher gewesen, wie sie die Operation COS() auf eine Ungleichung auswirkt und hatte damals dann in Anbetracht der immer näher rückenden Deadline gesagt: „Ach komm, da ist ne potente Xeon-CPU hinter und bei 1000 Datensätzen zeigt sich das ganze noch recht flott.“

Weißt du da mehr? Ich habe mir seit dem nicht mehr so richtig darüber Gedanken gemacht. ;-)

#513.06.2012 17:11

Jue:

... ich habe eben noch mal kurz drüber nachgedacht:

Momentan schaue ich ja, wo der Winkel zwischen Umkreismittelpunktsvektor und Kandidantenvektor eine gewisse Größe nicht überschreitet. Sag wir, ich suche alle Orte, bei denen der Winkel zwischen 0° und 90° liegt. Also ein Umkreisradius von 10.000km.
Mathematisch ausgedrückt: Maximalwinkel 90° >= Kandidatenwinkel

Wenn ich nun aber den COS um beides herumlege und mir den Verlauf ansehe, merke ich, dass jener bei 1 ( cos 0° ) losgeht und letztendlich bei -1 ( cos 180° ) ankommt.
Mathematisch: cos( Maximalwinkel 90° ) <= cos( Kandidatenwinkel )

Das größer/gleich ist somit zum kleiner/gleich geworden.

Müsste so hinhauen, oder?

Schreibe den 6. Kommentar

Ältere Einträge
CC BY-NC-SA 3.0