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

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.

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.

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.

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.

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 ...

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.

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

könnte versuchen, csv2shp csv2shp

Neugierig, in welcher Branche sich der Milliardenpunkt CSV befindet?