Location Based Services sind ohne Frage auf dem Vormarsch - was liegt also näher, als auch auf der eigenen Website Geo-Services anzubieten. Im einfachsten Fall geht es darum, die Entfernung zwischen zwei Orten / Personen zu messen, um Trefferlisten nach Nähe zu sortieren.

Genau dieses Problem hatten wir - in diesem Artikel zeige ich, wie man es löst.

Die Aufgabe

Bei unserer Aufgabenstellung geht es wie oben beschrieben darum, Trefferlisten zu sortieren: Die dem Webuser nächsten Personen sollen zuerst erscheinen. Eine hausnaummern-genaue Bestimmung der Geo-Koordinaten ist für uns nicht relevant; uns reicht die Postleitzahl. Dies hat zudem den Vorteil, das man auch vom Webuser nicht alle Daten erfragen muss - die Angabe der Postleitzahl genügt.

Die Idee besteht nun darin, zu allen Postleitzahlengebieten die geometrischen Mittelpunkte zu bestimmen. Anhand dieser Koordinaten lässt sich dann recht einfach die Distanz zwischen zwei Postleitzahl-Gebieten errechnen.

Bei meiner Recherche nach passenden Daten und Frameworks bin ich zunächst auf OpenGeoDB gestossen. In OpenGeoDB liegen nicht nur Postleitzahlen vor, sondern auch die Zuordnungen der Postleitzahlen zu den Verwaltungsstrukturen (Gemeinden, Bundesländern etc.). Leider sind jedoch die Ortsbezüge(Koordinaten) den Städten und nicht den Postleitzahlen zugeordnet - alle Bewohner Berlins haben die gleiche Koordinate. Das hilft natürlich wenig!

Freie Daten

Bei der weiteren Recherche bin ich dann auf MapBender gestossen. MapBender stellt Geo-daten übersichtlich in Maps zusammen, die zudem navigierbar sind. Für uns relevant: MapBender hat unter SourceForge die Geo-Daten der Postleitzahlen freigegeben (Im Bereich Data/PLZ nachschauen!): zu jeder Postleitzahl findet sich dort der Umriss des Gebietes als Polygonzug. Die Daten liegen als Postgres/PostGIS Daten vor (PostGIS erweitert postgres um GIS-Funktionen. PostGIS ist OpenGIS-konform).

Die Datenbank

Also: Postgres Installieren, PostGIS installieren, Daten einspielen! Mit den debian-postgres Paket kamen wir nicht weit - wir scheiterten immer an dem

undefiend symbol: nth

Als wir aber postgres aus den Quellen kompilierten, funktionierte alles fehlerfrei. Nach dem Einspielen der Postleitzahlen in PostGres liegen die Daten dann in der Tabellle post_code_areas:

postgis=# \d post_code_areas
                                TABLE "public.post_code_areas"
  COLUMN  |       Type        |                           Modifiers                           
----------+-------------------+---------------------------------------------------------------
 gid      | integer           | NOT NULL DEFAULT NEXTVAL('post_code_areas_gid_seq'::regclass)
 plz99    | character varying | 
 plz99_n  | integer           | 
 plzort99 | character varying | 
 the_geom | geometry          | 
Indexes:
    "post_code_areas_pkey" PRIMARY KEY, btree (gid)
CHECK constraints:
    "$1" CHECK (srid(the_geom) = 4326)
    "$2" CHECK (geometrytype(the_geom) = 'MULTIPOLYGON'::text OR the_geom IS NULL)

postgis=#

Aus Fläche mach Punkte

Mit den Polygonzügen der Postleitzahlen wollen wir nicht arbeiten; wir erlauben uns eine weitere Vereinfachnung: Wir reduzieren die Gebiete auf den geometrischen Mittelpunkt. Mit PostGIS ist dies ziemlich einfach; die Funktion centroid macht genau das gewünschte:

postgis=# select (point(centroid(the_geom))) from post_code_areas limit 20;
                point                
-------------------------------------
 (13.724898815155,51.05735206604)
 (13.7359957695007,51.0390377044678)
 (13.7440705299377,51.0672779083252)
 (13.8562636375427,51.0991764068604)
 (13.7630171775818,51.128568649292)
 (13.7324242591858,51.0800037384033)
 (13.7086491584778,51.09889793396)
 (13.6885514259338,51.0847797393799)
 (13.6662459373474,51.0664196014404)
 (13.7031207084656,51.0407161712646)
 (13.6616768836975,51.0388126373291)
 (13.7001214027405,51.0265789031982)
 (13.6995062828064,51.0130062103271)
 (13.7365202903748,51.014440536499)
 (13.7652487754822,51.0151195526123)
 (13.8000407218933,51.0164356231689)
 (13.7834916114807,51.0029048919678)
 (13.8136086463928,50.9896984100342)
 (13.8462166786194,50.9910869598389)
 (13.8025555610657,51.0345325469971)
(20 rows)

postgis=#

Zum bequemeren Arbeiten haben wir uns eine Tabelle gebaut, die nur noch Postleitzahlen und Geo-Koordinaten besitzt; Die Umrechnung der Geokoordinaten in Kugelkoordinaten ist bei der Berechnung der Entfernung hilfreich und wird gleich mit erledigt:

postgis=# CREATE TABLE plz_coord (
postgis(# breite float,
postgis(# laenge float,
postgis(# plz character varying
postgis(# );
CREATE TABLE
postgis=# insert into plz_coord select (point(centroid(the_geom)))[1]*pi()/180.0,(point(centroid(the_geom)))[0]*pi()/180.0, plz99 from post_code_areas;
INSERT 0 8270
postgis=# select * from plz_coord limit 20;                                                                                         
       breite      |      laenge       |  plz  
-------------------+-------------------+-------
 0.239544669741332 | 0.891118885701478 | 01067
 0.239738348128304 | 0.890799239797094 | 01069
 0.239879279281836 | 0.891292124327433 | 01097
 0.241837418374185 | 0.891848858143817 | 01099
 0.240209960659799 | 0.892361849540043 | 01109
 0.239676013515103 | 0.891514231958557 | 01127
 0.239261059734833 | 0.891843997874911 | 01129
  0.23891028813574 | 0.891597588899254 | 01139
 0.238520983925282 | 0.891277144046557 | 01157
 0.239164570081235 | 0.890828534568585 | 01159
 0.238441238896751 | 0.890795311634553 | 01169
 0.239112222321951 | 0.890581792697798 | 01187
 0.239101486453989 | 0.890344904522877 | 01189
 0.239747502744395 | 0.890369938236696 | 01217
 0.240248909390077 | 0.890381789303345 | 01219
  0.24085614339777 | 0.890404759067355 | 01237
 0.240567306937722 | 0.890168603261721 | 01239
 0.241092948348909 | 0.889938106673306 | 01257
 0.241662065864151 | 0.889962341438813 | 01259
 0.240900035620738 | 0.890720609967251 | 01277
(20 rows)

postgis=#

Entfernungen messen

Nun benötigen wir nur noch eine Formel zur Berechnung der Entfernung von zwei Kugelkoordinaten:

  dist = 6378.137 * ARCCOS( SIN(a.breite)*SIN(b.breite) + COS(a.breite)*COS(b.breite)*COS(b.laenge-a.laenge) )

Hat man zwei Postleitzahlen, so kann man dies direkt in SQL ausdrücken

postgis=# select 6378.137 * ACOS( SIN(a.breite)*SIN(b.breite) + COS(a.breite)*COS(b.breite)*COS(b.laenge-a.laenge) ) as dist from plz_coord a, plz_coord b where a.plz=10405 and b.plz=10551 ;
       dist       
------------------
 5.84523133900308
(1 row)

Prima: Das sind die Kilometer, die ich mit dem Rad von Moabit nach Prenzlberg radel!