die OSMF hat ja vor einiger Zeit den GPS Planet veröffentlicht. Da 2,8 Mrd Punkte etwas viel sind um die irgendwo einzulesen und die auch in einem komischen Format sind, habe ich dafür ein Programm geschrieben, dass die Daten filtert und wieder in normale Koordinaten umwandelt. Funktioniert auch, habs mit einer BBox getestet. Nun habe ich auch eine Filterung mit Bounding Polygon geschrieben, was aber ewig braucht…mal schauen obs jemals fertig wird. Ein Versuch mit ein paar Testpunkten war erfolgreich. Die “Is Point in Polygon” Routine ist ja speziell bei detailierten Poly-Files recht aufwendig zu berechnen und 2.8 Mrd Punkte sind fast doppelt so viele Nodes wie im OSM Planet. Ich dachte mir jetzt, ich filtere mit einer BBox aus B-Poly(maxX,maxY,minX,minY) vor bevor ich die Punkte in die “Is Point in Polygon” Funktion schiebe. Weiß jemand, wie osmosis das macht…(Umlaufzahl oder Kreuzungszahl) und welche Abkürzungen es da gibt? Ich habe leider nicht den Teil des Source Codes gefunden, der dafür verantwortlich ist
Wie wäre es, die Daten zuerst in eine vernünftige, statische Datenstruktur zu packen und darauf Suchoperationen zu starten? Das Einfachste wäre wohl ein einfacher Region-Quadtree mit effizienten Unterdatenstrukturen ab einem gewissen Level oder Gridfiles.
edit: Nochmal gelesen. Du möchtest also Polygone als Grenzen… Da macht es in der Tat lediglich Sinn, rechteckig vorzusortieren. Eventuell lässt sich das durch Ausnutzen der Datenstruktur optimieren (z.B. diagonale Zellen im Quadtree → wenige Schnittberechnungen).
Eine neue Datenstruktur zu erstellen macht mMn nur Sinn, wenn man die Daten regelmäßig braucht. Ohne Vorfiltern mit einer BBox um die Aussenbereiche des Bounding-Polygons, ist die Rechenzeit erheblich größer. Also beschleunigt das schonmal. In Tests war die Rechenzeit ca. 10 Mal größer, kann aber auch mehr sein. Im wesentlichen hängt das von der Anzahl der Nodes im Bounding Polygon ab (Pi mal Dauem: Bounding Polygon: Rechnschritte = Nodes in den Daten * Nodes im Bounding Polygon vs Rechenschritte=Nodes in Daten *4).
Jetzt ist mir noch eine Idee gekommen. Man könnte ja auch innerhalb des Bounding-Polygon eine Bounding-Box definieren. Alles was in dieser ist, muss dann auch nicht durch die aufwendige Point-In-Polygon Prüfung. Das muss ich aber mal nen Mathematiker fragen.
kleinstes Rechteck berechnen, wo das Polygon gerade so reinpasst. Dann größtes Rechteck suchen, das gerade so ins Polygon reinpasst. Die Differenz musst du dann gegen das Polygon prüfen.
Das hab ich schon. Ich habe den GPS Planet mal um ca. Faktor 1000 verkleinert. Sonst dauern die Tests ewig.
Mit Bounding Box ums Polygon ca. 1 Sekunden
thomas@suncobalt:~/pascal/bbox$ time ./GPXParserV2 -i=gps_parsed.txt -o=eislingen.txt -p=eislingen.poly
Input File : gps_parsed.txt
Output File : eislingen.txt
Outputformat : LATLON
Filter : Bounding Polygon
|-> Filename : eislingen.poly
|-> Vertices in BPoly : 196
BBox prefilter : yes
Top : 48.7277500
Left : 9.6810880
Bottom : 48.6735800
Right : 9.7528120
2,192,988 Points in total, filtered 243
real 0m1.084s
user 0m1.012s
sys 0m0.068s
Ohne Vorfilterung mittels BBox 13 Sekunden.
thomas@suncobalt:~/pascal/bbox$ time ./GPXParserV2 -i=gps_parsed.txt -o=eislingen.txt -p=eislingen.poly -t
Input File : gps_parsed.txt
Output File : eislingen.txt
Outputformat : LATLON
Filter : Bounding Polygon
|-> Filename : eislingen.poly
|-> Vertices in BPoly : 196
BBox prefilter : no
2,192,988 Points in total, filtered 243
real 0m13.101s
user 0m12.997s
sys 0m0.100s
Genau das ist mathematisch nicht trivial. Evtl mache ich es mit Brute Force…grob gesagt, einfach ein paar hundert Rechtecke suchen und das größte nehmen, das gerade so ins Polygon reinpasst-
Als ich 2010 mit osm-Programmierung anfing, war Java für mich die erste Wahl , und ich hab das trotz kontroverser Diskussionen nicht bereut.
Ich sehe auch keinen Sinn darin, einen in PostGIS und auch in Osmosis enthaltenen Algorithmus neu zu erfinden.
nunja, ich schreibe in Pascal. Habe nach 15 Jahren wieder mit dem Programmieren angefangen. Ich wollte halt ein paar Punkte aus dem GPS Planet filtern. Jetzt kann ich mich hinsetzen und warten bis jemand mit Java Kenntnissen ein Programm mit den “perfekten” Routinen schreibt oder ich mache es selber. Und es funktioniert ja recht gut. Es geht momentan nur noch darum, es schneller zu machen…
Danke für die vielen Tipps. Ich werde mal schauen, was sich davon umsetzen lässt. Am Anfang werde ich mal schauen, ob osmosis mit den existierenden Libaries wirklich so viel schneller ist.
Was Du vermutlich meinst, ist “nur” eine Seminararbeit.
Das Paper selbst stammt von Daniels & Milenkovic.
Da von 1998 dürfte - irgendwo - auch sourcecode dazu rumkullern …
Egal wie schnell der Filter auch sein mag: alleine das Einlesen von 2 Milliarden Punkten braucht Zeit.
Also erst gar nicht einlesen, was nicht gebraucht wird.
1a. Ohne Umstrukturierung des GPXPlanets:
Die Punkte sind nach Breite sortiert. Also einmalig einen Index aufbauen mit der Startposition in der Datei für jedes Grad oder jedes zentel Grad Breite.
Beim Filtern dann gleich ein Seek auf MinLat; und Abbruch wenn der erste Punkt mit lat>Maxlat gelesen wird.
1b. Mit Umstrukturierung:
Mit mehr Aufwand kann man die Punkte der Quelldatei in Kacheln (zB. 1024x1024) umsortieren, dazu einen Index anlegen, und dann nur noch die Kacheln lesen, die einen nicht leeren Schnitt mit der Boundingbox haben.
Damit sollte man auch effizient einen Kreis oder ein Land aus dem ganzen Planetfile lesen können.
Schnellfilterung auf außerhalb der Boundingbox wie bereits gemacht.
Boundingbox aufteilen.
Die Boundingbox wird in N×M Rechtecke aufgeteilt, und für jedes Rechteck wird die Information “ganz im Poly”, “teils im poly”, “garnicht im Poly” abgelegt, initialisiert auch “unbekannt”.
Für jeden Punkt wird das Teilrechteck bestimmt.
→ ganz im Poly: Punkt aufnehmen
→ nicht im Poly: Punkt verwerfen.
→ teils im Poly: normalen inPolygon-Test aufrufen.
→ unbekannt:
→→ Prüfen ob eine der Grenzen des Rechtecks sich mit einer der Polygonlinien schneidet (einfacher Algorithmus).
→→ Prüfen ob Punkt im Polygon.
→→ wenn Schnitt → Status=teils im Polygon
→→ wenn kein Schnitt und Punkt im Polygon → Status=ganz im Polygon
→→ wenn kein Schnitt und Punkt nicht im Polygon → Status=garnicht im Polygon.
Man spart in ganz oder garnicht Rechtecken ab dem zweiten Punkt eine “inPoly”-Abfrage, das bezahlt man aber mit der “gibt es Schnitt mit Polygon-Linie”-Routine je nicht leerem Rechteck. Die beste Granularität wird man ausprobieren müssen.
Ich brauche keinen Spoiler. Ich habe ein (fertiges) Programm, das mir aus dem OSM GPS Planet die gewünschten Daten ausschneidet (Bounding Polygon oder Bounding Box) und dieses komische Format des OSM Planet GPS in “normale” Koordinaten (LatLon oder LonLat oder osm-xml) umwandelt… kann man einstellen.
Wenn es jemand braucht, kann ich es auch verteilen. Aber ich will es schneller machen…auch ohne Java Libaries. Es macht Spaß solche Datenmengen behandeln zu können und man lernt viel, wenn man sowas schreibt. Gerne schaue ich in Java Libs rein ob ich Abkürzungen implementieren kann. Momentan ist der Code etwas über 500 Zeilen und wenn man ein Extrembeispiel nimmt (ein DE Bounding Polygon mit 10x mehr Punkten als das der Geofabrik), dann ist osmosis noch doppelt so schnell (ohne Zeit, die es dauert den GPS Planet in osm-xml umzuwandeln). Durch das Vorfiltern eines Bounding Polygons mit einer äusseren Box habe ich das Programm um den Faktor 10 beschleunigt. Wenn ich jetzt mit einer inneren Box noch den Faktor 2 hinbekomme, liege ich mit osmosis gleich auf… Code-Optimierung und evtl Multithread mal ausser acht gelassen.
@all, Netzwolf…auch wenn ich nicht drauf eingehe, ich grübele über jeden Tipp. Wie gesagt, das ist mein erstes richtige Programm seit 15 Jahren. Da geht alles etwas langsamer Ehrlich gesagt bin ich schon froh, dass es überhaupt funktioniert.
Ehrlich gesagt bin ich schon froh, dass es überhaupt funktioniert.