Fürchterlich komplizieter Umkreissuche

[Rolf B]

Hallo Raketenwissenschaftler,

einerseits hast Du Recht - die Voreinschränkung mit einem „Quadrat“ hilft schonmal gewaltig. D.h. man findet zuerst mal die Intervalle, in denen Längen- und Breitengrad liegen müssen, um überhaupt im Radius enthalten sein zu können.

Für die Breite ist es relativ einfach, der Abstand zwischen zwei Breitengraden beträgt näherungsweise 40000km / 360 = 111.111 km. Bei den Längengraden ist es komplizierter, der Abstand zwischen zwei Längengraden hängt von der Breite ab. Der Radius der "Erdscheibe" in einem bestimmten Breitengrad $$\phi$$ beträgt $$r_\phi = r_E \cos \phi$$, und da der Umfang zum Radius proportional ist, gilt das auch für den Umfang der Erde an diesem Breitengrad. Demnach ist der Abstand zweier Längengrade auf Breite $$\phi$$: $$l_\phi = 111,111 \cos \phi$$. Bleibt man in Deutschland, liegt der südlichste Punkt an der Oberstdorfer Trifthütte, $$\phi=47{,}27^\circ$$, und einem Längengradabstand von 75,4km. Der nördlichste Punkt liegt im Rickelsbüller Koog auf 54,9 Grad mit einem Längengradabstand von 63,9km.

Andererseits kann man sich durch passende Vorausberechnungen auch viel Rechnerei sparen. Ich schrub das ja schon anderswo.

Vorausberechnung Stufe 1: sin und cos von Länge und Breite jedes gespeicherten Ortes bestimmen und in der Ortstabelle abspeichern. Das kann man ab MySQL 5.7 mit generated columns lösen (aber STORED, nicht VIRTUAL, sonst ist's witzlos). Ebenso sin und cos von Länge und Breite des Referenzortes bestimmen.

Die Abstandsformel für Orthodrome steht in der Wikipedia. Sie basiert auf dem Zentriwinkel, d.h. dem Winkel zwischen den beiden Geraden vom Erdmittelpunkt zu den beiden Punkten, deren Abstand zu finden ist. Der Cosinus des Zentriwinkels berechnet sich so:

$$\cos \zeta = \sin(\phi_A) \cdot \sin(\phi_B) + \cos(\phi_A) \cdot \cos(\phi_B) \cdot \cos(\lambda_B - \lambda_A)$$

Das unschöne an dieser Formel ist, dass man den Cosinus der Differenz der Längengrade ausrechnen muss, was sich nicht vorberechnen lässt. Aber Mathematik ist ja die Kunst, das Rechnen zu vermeiden 😉, darum wenden wir zuerst einmal das Additionstheorem für $$\cos(\beta - \alpha)$$ an:

$$\cos \zeta = \sin(\phi_A) \cdot \sin(\phi_B) + \cos(\phi_A) \cdot \cos(\phi_B) \cdot (\cos\lambda_B \cos\lambda_A + \sin\lambda_B \sin\lambda_A) $$

Die Distanz d der beiden Punkte A und B ist $$d = r\cdot\zeta$$, mit r=Erdradius und $$\zeta = \arccos(...)$$ (im Bogenmaß) entsprechend der obigen Formel.

Nun muss man alle Orte finden, für die d einen bestimmten Abstand A unterschreitet. Das führt wieder zu einer Lästigkeit, denn das Berechnen des arccos in der SQL Abfrage ist ebenfalls zeitraubend. Also nutzen wir den bekannten Umstand, dass die arccos-Funktion zwischen x=-1 und x=1 streng monoton fallend ist, dementsprechend gilt

$$ \displaystyle \quad\quad r \cdot \zeta \lt A$$
$$ \displaystyle \Longleftrightarrow \zeta \lt \frac{A}{r}$$
$$ \displaystyle \Longleftrightarrow \cos \zeta \gt \cos\frac{A}{r}$$

Diesen Cosinus kann man wieder vorausberechnen. In der Query berechnet man dann nur noch den Wert innerhalb des arccos-Aufrufs, was 5 Multiplikationen und 2 Additionen sind, und vergleicht ihn mit dem vorausberechneten Referenzwert. Achtung, es ist ein GRÖSSER-Vergleich, weil arccos eine fallende Funktion ist.

Die Hardcore-Lösung wäre übrigens, eine Tabelle mit allen möglichen PLZ Paaren zu bilden und darin den Abstand dieser PLZ zu speichern. Das dauert eine Weile, aber das kann man nach jedem Update der PLZ DB einmal machen und die Abfragen gehen dann Ratz-Superfatz direkt auf den Index, wenn man nach "Von-PLZ" und "Abstand" Indexiert. Das sind 67 Millionen bis 841 Millionen Sätze, je nachdem, ob man nur die 8200 Hauszustellungs-PLZ will oder auch die 16500 Postfach- oder 3100-Großempfänger-PLZ berücksichtigen will. Ob sich das lohnt, hängt vom Mengengerüst der erwarteten Abfragen ab.

Rolf