Auf der Suche nach der schnellsten Lösung für die Point-in-Polygon-Analyse von 200 Millionen Punkten

Ich habe eine CSV mit 200 Millionen Beobachtungen in folgendem Format:

id,x1,y1,x2,y2,day,color
1,"-105.4652334","39.2586939","-105.4321296","39.2236632","Monday","Black"
2,"-105.3224523","39.1323299","-105.4439944","39.3352235","Tuesday","Green"
3,"-104.4233452","39.0234355","-105.4643990","39.1223435","Wednesday","Blue"

Für jeden Koordinatensatz (x1 / y1 und x2 / y2) möchte ich den US-amerikanischen Zensus-Trakt oder den Zensus-Block zuweisen, in den er fällt (das TIGER-Shapefile des Zensus-Trakts habe ich hier heruntergeladen: ftp://ftp2.census.gov/ geo / tiger / TIGER2011 / TRACT / tl_2011_08_tract.zip ). Daher muss ich für jede Beobachtung zweimal eine Punkt-in-Polygon-Operation ausführen. Es ist wichtig, dass die Übereinstimmungen sehr genau sind.

Was ist der schnellste Weg, dies zu tun, einschließlich der Zeit zum Erlernen der Software? Ich habe Zugriff auf einen Computer mit 48 GB Speicher - für den Fall, dass dies eine relevante Einschränkung darstellt.

In mehreren Threads wird die Verwendung von PostGIS oder Spatialite empfohlen (Spatialite sieht einfacher aus - ist es jedoch so effizient wie PostGIS?). Wenn dies die besten Optionen sind, muss unbedingt ein Raumindex (RTree) erstellt werden? Wenn ja, wie macht man das (zB mit dem Census Tract Shapefile)? Ich wäre sehr dankbar für alle Empfehlungen, die Beispielcode (oder einen Zeiger auf Beispielcode) enthalten.

Mein erster Versuch (vor dem Auffinden dieser Site) bestand darin, mit ArcGIS eine räumliche Verknüpfung (nur x1 / y1) von Teilproben der Daten (100.000 Punkte) im US-Volkszählungsblock durchzuführen. Das hat über 5 Stunden gedauert, bis ich den Prozess beendet habe. Ich hoffe auf eine Lösung, die in weniger als 40 Stunden Rechenzeit für den gesamten Datensatz implementiert werden kann.

Entschuldigung, dass Sie eine Frage gestellt haben, die bereits gestellt wurde. Ich habe die Antworten gelesen und habe mich gefragt, wie ich die Empfehlungen umsetzen soll. Ich habe noch nie SQL, Python, C und ArcGIS verwendet - ich bin ein absoluter Anfänger.

ST_DWithin war in meinem Test schneller als ST_Intersects. Das ist überraschend, zumal der vorbereitete Geometrie-Algorithmus in solchen Fällen Einzug halten soll. Ich denke, es besteht die Möglichkeit, dass dies viel schneller ist, als ich hier gezeigt habe.

Ich habe noch einige Tests gemacht und zwei Dinge haben die Geschwindigkeit fast verdoppelt. Zuerst habe ich einen neueren Computer ausprobiert, aber immer noch einen ganz normalen Laptop, vielleicht mit Ausnahme von SATA3-SSD-Disks.

Dann dauerte die Abfrage unten 18 Sekunden anstelle von 62 Sekunden auf dem alten Laptop. Als nächstes stellte ich fest, dass ich mich völlig geirrt hatte, als ich schrieb, dass der Index auf der Punktetabelle nicht notwendig sei. Mit diesem Index verhielt sich ST_Intersects wie erwartet und die Dinge wurden sehr schnell. Ich habe die Punktzahl in der Punktetabelle auf 1 Million Punkte erhöht und die Abfrage:

CREATE TABLE points_ct AS
SELECT imported_ct.gid as ct_gid, t.gid as point_id 
FROM imported_ct , t WHERE ST_Intersects(imported_ct.geom , t.geom);

Läuft in 72 Sekunden. Da es 1249 Polygone gibt, werden 1249000000 Tests in 72 Sekunden durchgeführt. Das macht ungefähr 17000000 Tests pro Sekunde. Oder testen Sie fast 14000 Punkte gegen alle Polygone pro Sekunde.

Ab diesem Test sollten Ihre 400000000 zu testenden Punkte ungefähr 8 Stunden dauern, ohne dass die Last auf mehrere Kerne verteilt werden muss. PostGIS hört nie auf, mich zu beeindrucken :-)

Um das Ergebnis zu visualisieren, können Sie der resultierenden Tabelle zunächst die Punktgeometrie hinzufügen, sie beispielsweise in QGIS öffnen und sie mit eindeutigen Werten im Feld import_ct formatieren.

Zweitens, ja, Sie können die Punkte, die außerhalb eines Polygons liegen, auch erhalten, indem Sie eine rechte (oder linke) Verknüpfung wie folgt verwenden:

CREATE TABLE points_ct AS
SELECT imported_ct.gid as ct_gid, t.gid as point_id 
FROM imported_ct right join t ON ST_Intersects(imported_ct.the_geom , t.geom);

Ich habe einige Tests durchgeführt, um zu überprüfen, ob es möglich scheint, PostGIS.

Zuerst etwas, was ich nicht verstehe. Sie haben zwei Punkte pro Reihe. Befinden sich immer beide Punkte im selben Polygon? Dann reicht es aus, die Berechnungen an einem der Punkte durchzuführen. Wenn sie sich in zwei verschiedenen Polygonen befinden können, müssen Sie eine Möglichkeit finden, eine Punktreihe mit zwei Polygonen zu verbinden.

Aus den Tests scheint es machbar, aber Sie benötigen möglicherweise eine kreative Lösung, um die Last auf mehr als einen CPU-Kern zu verteilen.

Ich habe auf einem 4 Jahre alten Laptop mit Dual-Core-Centrino-CPU (etwa 2,2 GHz, glaube ich) und 2 GB RAM getestet. Wenn Sie 48 BG RAM haben, haben Sie vermutlich auch viel mehr CPU-Leistung.

Was ich getan habe, war eine zufällige Punktetabelle mit 100000 Punkten wie folgt zu erstellen:

CREATE TABLE t AS
WITH r AS
(SELECT ST_Extent(the_geom)::geometry ext FROM imported_ct)
SELECT ST_Point(x,y) AS geom FROM 
(SELECT GENERATE_SERIES(1,100000)) s,
(SELECT ST_Xmin(ext)+(random()*(ST_Xmax(ext)-ST_Xmin(ext))) x, ST_Ymin(ext)+(random()*(ST_Ymax(ext)-ST_Ymin(ext))) y FROM r
) f;

Dann füge ich ein gid hinzu wie:

ALTER TABLE t ADD COLUMN GID SERIAL;

Dann läuft:

CREATE TABLE points_ct AS
SELECT imported_ct.gid as ct_gid, t.gid as point_id FROM imported_ct , t WHERE ST_Dwithin(imported_ct.the_geom , t.geom,0);

dauert etwa 62 Sekunden (vergleiche mit deinem ArcGIS-Ergebnis mit der gleichen Anzahl von Punkten). Das Ergebnis ist eine Tabelle, die die Punkte in meiner Tabelle t mit der Tabelle mit dem Zensus-Trakt verbindet.

Mit dieser Geschwindigkeit erledigen Sie 200 Mühlpunkte in ungefähr 34 Stunden. Wenn es also ausreicht, einen Punkt zu überprüfen, kann mein alter Laptop dies mit einem Kern tun.

Wenn Sie jedoch beide Punkte überprüfen müssen, ist dies möglicherweise schwieriger.

Dann können Sie die Last manuell auf mehrere Kerne verteilen, indem Sie mehrere Sitzungen mit der Datenbank starten und verschiedene Abfragen ausführen.

In meinem Beispiel mit 50000 Punkten und zwei CPU-Kernen habe ich versucht:

CREATE TABLE t1 as
SELECT imported_ct.gid as ct_gid, t.gid as point_id FROM imported_ct , t WHERE t.gid >50000 and  ST_Dwithin(imported_ct.the_geom , t.geom,0);

in einer db-session zur selben zeit wie running:

CREATE TABLE t2 as
SELECT imported_ct.gid as ct_gid, t.gid as point_id FROM imported_ct , t WHERE t.gid <=50000 and  ST_Dwithin(imported_ct.the_geom , t.geom,0);

auf einer anderen db-session.

Das hat ungefähr 36 Sekunden gedauert, ist also etwas langsamer als das erste Beispiel, was wahrscheinlich vom gleichzeitigen Beschreiben der Disc abhängt. Aber da meine Kerne gleichzeitig arbeiten, dauerte es nicht länger als 36 Sekunden.

Zur Vereinigung der Tabellen t1 und t2 wird versucht:

CREATE TABLE t3 AS 
SELECT * FROM t1
UNION ALL
SELECT * FROM t2;

mit etwa einer halben Sekunde.

Mit frischerer Hardware und einer Verteilung der Last auf viele Kerne sollte dies absolut möglich sein, auch wenn die reale Welt langsamer als der Testfall sein wird.

Erwähnenswert ist, dass das Beispiel von Linux (Ubuntu) stammt. Die Verwendung von Windows wird eine andere Geschichte sein. Aber ich habe alle anderen täglichen Anwendungen ausgeführt, so dass der Laptop von vor ziemlich stark geladen ist. Das könnte also den Windows-Fall recht gut simulieren, ohne etwas anderes als pgadmin zu öffnen.

Der wahrscheinlich einfachste Weg ist mit PostGIS. Im Internet gibt es einige Tutorials zum Importieren von CSV / TXT-Punktdaten in PostGIS. Link1

Ich bin mir nicht sicher, ob die Point-in-Polygon-Suche in PostGIS durchgeführt wird. Es sollte schneller als ArcGIS sein. Der von PostGIS verwendete räumliche GIST-Index ist ziemlich schnell. Link2 Link3

Sie können auch den MongoDB-Geodatenindex testen . Dies erfordert jedoch etwas mehr Zeit, um loszulegen. Ich glaube, dass MongoDB sehr schnell sein kann. Ich habe es nicht mit Point-in-Polygon-Suchen getestet, kann also nicht sicher sein.