Der schnellste Weg, um eine Punkt-CSV mit einem Polygon-Shapefile räumlich zu verbinden

19

Ich habe eine 1-Milliarde-Punkte-CSV-Datei und ein Shapefile mit etwa 5.000 Polygonen. Was wäre der schnellste Weg, um Punkte und Polygone räumlich zu verbinden? Für jeden Punkt muss ich die enthaltene Polygon-ID abrufen. (Polygone überlappen sich nicht.)

Normalerweise würde ich beide Datensätze in PostGIS laden. Gibt es einen schnelleren Weg, um die Arbeit zu erledigen?

Ich suche eine Open-Source-Lösung.

Underdunkel
quelle

Antworten:

16

Wenn "am schnellsten" die Menge Ihrer aufgewendeten Zeit einschließt , hängt die Lösung davon ab, mit welcher Software Sie vertraut sind und sie schnell verwenden können. Die folgenden Ausführungen konzentrieren sich konsequent auf Ideen zur Erzielung möglichst schneller Rechenzeiten .

Wenn Sie ein vorgefertigtes Programm verwenden, können Sie die Polygone mit ziemlicher Sicherheit vorverarbeiten, um eine Point-in-Polygon-Datenstruktur wie einen KD-Baum oder einen Quadtree einzurichten, deren Leistung in der Regel 0 (log (V ) * (N + V)) wobei V die Gesamtzahl der Eckpunkte in den Polygonen und N die Anzahl der Punkte ist, da die Datenstruktur mindestens O (log (V) * V) benötigt, um erstellt zu werden, und dann müssen für jeden Punkt zu einem Punktpreis O (log (V)) abgetastet werden.

Sie können wesentlich bessere Ergebnisse erzielen, indem Sie zuerst die Polygone gittern und dabei die Annahme ausnutzen, dass keine Überlappungen auftreten. Jede Gitterzelle befindet sich entweder vollständig in einem Polygoninneren (einschließlich des Inneren des "universellen Polygons"), in welchem ​​Fall die Zelle mit der ID des Polygons beschriftet wird, oder sie enthält eine oder mehrere Polygonkanten. Die Kosten für diese Rasterung, die der Anzahl der Gitterzellen entsprechen, auf die beim Rastern aller Kanten verwiesen wird, betragen O (V / c), wobei c die Größe einer Zelle ist, die implizite Konstante in der Big-O-Notation jedoch klein ist.

(Ein Vorteil dieses Ansatzes ist, dass Sie Standard-Grafikroutinen ausnutzen können. Wenn Sie beispielsweise ein System haben, das (a) die Polygone auf einem virtuellen Bildschirm mit (b) einer bestimmten Farbe für jedes Polygon zeichnet und (c) zulässt Wenn Sie die Farbe eines Pixels lesen möchten, das Sie ansprechen möchten, haben Sie es gemacht.)

Wenn dieses Raster aktiviert ist, überprüfen Sie die Punkte vorab, indem Sie die Zelle berechnen, die jeden Punkt enthält (eine O (1) -Operation, die nur wenige Takte erfordert). Sofern die Punkte nicht um die Polygongrenzen gruppiert sind, verbleiben normalerweise nur etwa O (c) -Punkte mit mehrdeutigen Ergebnissen. Die Gesamtkosten für den Aufbau des Netzes und die Vorabschirmung betragen daher O (V / c + 1 / c ^ 2) + O (N). Sie müssen eine andere Methode (wie eine der bisher empfohlenen) verwenden, um die verbleibenden Punkte (dh diejenigen, die sich in der Nähe der Polygongrenzen befinden) zu verarbeiten. Die Kosten hierfür betragen O (log (V) * N * c). .

Wenn c kleiner wird, befinden sich immer weniger Punkte in derselben Gitterzelle mit einer Kante, und daher erfordern immer weniger die nachfolgende O-Verarbeitung (log (V)). Demgegenüber ist es notwendig, O (1 / c ^ 2) Gitterzellen zu speichern und O (V / c + 1 / c ^ 2) Zeit damit zu verbringen, die Polygone zu rastern. Daher wird es eine optimale Gittergröße geben. C. Unter Verwendung dieser Methode betragen die Gesamtkosten für die Berechnung 0 (log (V) * N), aber die implizite Konstante ist in der Regel die Methode betragen aufgrund der O (N) -Geschwindigkeit des Vorscreenings in der kleiner als bei Verwendung der Festprozeduren.

Vor 20 Jahren habe ich diesen Ansatz getestet (unter Verwendung gleichmäßig verteilter Punkte in ganz England und vor der Küste und unter Ausnutzung eines relativ groben Gitters von rund 400K-Zellen, das von den damaligen Videopuffern angeboten wurde) und zwei Größenordnungen schneller als mit dem besten veröffentlichten Algorithmus, den ich konnte, erzielt finden. Selbst wenn die Polygone klein und einfach sind (wie Dreiecke), können Sie sich auf eine Beschleunigung in der Größenordnung verlassen.

Nach meiner Erfahrung war die Berechnung so schnell, dass der gesamte Vorgang durch die Daten-E / A-Geschwindigkeiten und nicht durch die CPU begrenzt wurde. Vorausgesetzt, dass E / A der Engpass sein könnte, erzielen Sie die schnellsten Ergebnisse, wenn Sie die Punkte in einem so komprimierten Format wie möglich speichern, um die Datenlesezeiten zu minimieren. Überlegen Sie auch, wie die Ergebnisse gespeichert werden sollen, damit Sie die Anzahl der Schreibvorgänge auf der Festplatte begrenzen können.

whuber
quelle
6
Sehr guter Zeitpunkt für die Realisierung der Lösung im Vergleich zur Rechenzeit. Eine lange Zeit in Anspruch zu nehmen, um zu einer optimalen Lösung zu gelangen, ist nur dann von Vorteil, wenn Sie diese Einsparungen durch die Optimierung realisieren (insbesondere aus Sicht des Arbeitgebers).
Sasa Ivetic
5

Ich für meinen Teil würde wahrscheinlich CSV-Daten in eine shp- Datei laden und dann ein Python-Skript mit shapefile und shapely schreiben , um die enthaltene Polygon-ID zu erhalten zu aktualisieren.

Ich weiß nicht, ob Geotools und JTS schneller als Shapefile / Shapely sind ... Hab keine Zeit, es zu testen!

Bearbeiten : Übrigens ist die CSV-Konvertierung in das Shapefile-Format wahrscheinlich nicht erforderlich, da Werte leicht formatiert werden können, um mit räumlichen Objekten aus Ihrem Polygon-Shapefile getestet zu werden.

simo
quelle
4
Ich lade die Daten direkt mit einem CSV-Reader und fülle einen räumlichen Rtree- Index auf. Die Kombination von Rtree und Shapely hat eine beeindruckende Leistung (viel besser als PostGIS; ich kann es nicht mit JTS vergleichen, da ich kein Java kann).
Mike T
2
Gute Idee, vorausgesetzt, Sie müssen nicht alle 1b-Punkte auf einmal speichern. Bei einem Minimum von 16 Bytes pro Punkt (X / Y) sehen Sie Daten im Wert von 16 GB. Wenn Rtree den Index auf lokalem Speicher erstellt, wird die Leistung definitiv verbessert. Das Importieren von 1b-Punkten in ein einzelnes Shapefile funktioniert ebenfalls nicht. Die Status-Shapefiles der OGR-Spezifikation sind auf 8 GB beschränkt (4 GB werden empfohlen). Eine einzelne Punktform verwendet 20 Bytes.
Sasa Ivetic
4

Am Ende habe ich die Polygone in ein Raster konvertiert und an den Punktpositionen abgetastet. Da sich meine Polygone nicht überlappten und keine hohe Genauigkeit erforderlich war (Polygone stellten Landnutzungsklassen dar und ihre Grenzen galten ohnehin als ungewiss), war dies die zeiteffizienteste Lösung, die ich finden konnte.

Underdunkel
quelle
3

Ich würde schnell ein kleines Java - Programm auf der Basis schreiben Shape - Datei Leser von GeoTools und der Betrieb enthält der JTS . Ich weiß nicht, wie schnell es sein kann ...

julien
quelle
1
Wenn Sie die Daten in PostGIS haben, können GeoTools Hauptindizes usw. verwenden.
Ian Turton
3

Verwenden Sie Spatialite .

Laden Sie die GUI herunter. Sie können sowohl Shapefile als auch CSV als virtuelle Tabellen öffnen. Dies bedeutet, dass Sie sie nicht tatsächlich in die Datenbank importieren, sondern sie als Tabellen anzeigen und dass Sie sie schnell verbinden und abfragen können, wie Sie möchten.

Sean
quelle
3

Sie können es ziemlich schnell mit OGR in C / C ++ / Python machen (Python sollte das langsamste der 3 sein). Durchlaufen Sie alle Polygone und setzen Sie einen Filter für die Punkte. Durchlaufen Sie die gefilterten Punkte und Sie wissen, dass jeder der Punkte, die Sie durchlaufen, zum aktuellen Polygon gehört. Hier ist ein Beispielcode in Python mit OGR, der die Polygone durchläuft und die Punkte entsprechend filtert. C / C ++ - Code sieht dem sehr ähnlich, und ich würde mir vorstellen, dass Sie eine deutliche Geschwindigkeitssteigerung gegenüber Python erzielen. Sie müssen einige Codezeilen hinzufügen, um die CSV zu aktualisieren:

from osgeo import ogr
from osgeo.gdalconst import *

inPolyDS = ogr.Open("winnipeg.shp", GA_ReadOnly)
inPolyLayer = inPolyDS.GetLayer(0)
inPointDS = ogr.Open("busstops.vrt", GA_ReadOnly)   
inPointLayer = inPointDS.GetLayerByName("busstops")

inPolyFeat = inPolyLayer.GetNextFeature()
while inPolyFeat is not None:
  inPtFeat = inPointLayer.GetNextFeature()
  while inPtFeat is not None:
    ptGeom = inPtFeat.GetGeometryRef()
    # Do work here...

    inPtFeat = inPointLayer.GetNextFeature()

  inPolyFeat = inPolyLayer.GetNextFeature()

VRT-Datei (busstops.vrt):

<OGRVRTDataSource>
  <OGRVRTLayer name="busstops">
    <SrcDataSource>busstops.csv</SrcDataSource>
    <GeometryType>wkbPoint</GeometryType>
    <LayerSRS>WGS84</LayerSRS>
    <GeometryField encoding="PointFromColumns" x="X" y="Y" reportSrcColumn="FALSE" />
  </OGRVRTLayer>
</OGRVRTDataSource>

CSV-Datei (busstops.csv):

FID,X,Y,stop_name
1,-97.1394781371062,49.8712241633646,Southbound Osborne at Mulvey

CSVT-Datei (busstops.csvt, OGR benötigt sie zum Identifizieren von Spaltentypen, andernfalls wird der räumliche Filter nicht ausgeführt):

Integer,Real,Real,String
Sasa Ivetic
quelle
2
Geht das nicht 5000 Mal durch 1 Milliarde Punkte (einmal für jedes Polygon)?
Underdunkel
Ein räumlicher Index ist ein absoluter Index Muss . Ich habe Rtree schon einmal erwähnt, und ich werde es noch einmal erwähnen!
Mike T