r Punkte in Polygonen

Question

Ich habe eine Million Punkte und eine große Form Datei-8GB-die zu groß ist, um in den Speicher in R auf meinem System zu laden. Die Shape-Datei ist Single-Layer, also ein gegebenes x, y wird höchstens ein Polygon treffen - solange es nicht genau an einer Grenze ist! Jedes Polygon ist mit einem severity - z.B. 1, 2, 3. Ich verwende R auf einem 64-Bit-Ubuntu-Rechner mit 12 GB RAM.

Was ist der einfachste Weg, um die Lage zu sein, der Datenrahmen zu dem Polygon zu „Tag“ severity so, dass ich ein data.frame mit einer zusätzlichen Spalte erhalten, das heißt x, y, severity?

Answer 1

Nur weil alles, was Sie haben, ein Hammer ist, bedeutet das nicht, dass jedes Problem ein Nagel ist.

Laden Sie Ihre Daten in PostGIS, erstellen Sie einen räumlichen Index für Ihre Polygone und führen Sie ein einzelnes räumliches SQL-Overlay aus. Exportiere die Ergebnisse zurück nach R.

Übrigens ist die Aussage, dass die Shapefile 8 Gb ist, keine sehr nützliche Information. Shapefiles bestehen aus mindestens drei Dateien, dem .shp, der die Geometrie darstellt, der .dbf, die die Datenbank ist, und dem .shx, der die beiden verbindet. Wenn Ihre DBF-Datei 8 GB groß ist, können Sie die Shapes selbst leicht lesen, indem Sie sie durch eine andere DBF-Datei ersetzen. Auch wenn die .shp 8 GB ist, können es nur drei Polygone sein, in diesem Fall könnte es einfach sein, sie zu vereinfachen. Wie viele Polygone hast du und wie groß ist der SHP-Teil des Shapefiles?

Answer 2

Ich denke, Sie sollten die Daten vorverarbeiten und eine Struktur erstellen, die mögliche Polygone für rechteckige Bereiche in einem Raster auflistet. Auf diese Weise können Sie die Polygone reduzieren, die Sie gegen die Punkte prüfen müssen, und die zusätzliche Struktur wird in den Speicher passen, da sie nur Zeiger auf die Polygone hat.

Hier ein Bild ist zu illustrieren.

Sie wollen prüfen, welche der gelbe Punkt Polygon ist Sie in der Regel gegen alle Polygone überprüfen würde, aber mit der Netzoptimierung (die orangefarbenen Linien, zeichnete nicht das gesamte Gitter, nur eines seiner Felder) Sie müssen nur die gefüllten Polygone überprüfen, da sie alle innerhalb oder teilweise innerhalb des Gitterfeldes sind. Ein ähnlicher Weg wäre, nicht alle Polygondaten im Speicher zu speichern, sondern nur die Polygon-Begrenzungsrahmen, die nur 2 anstelle von 3 X/Y-Paaren für jedes Polygon benötigen (und einen zusätzlichen Zeiger auf die tatsächlichen Polygon-Daten)), aber das spart nicht so viel Platz wie der erste Vorschlag.

Answer 3

Ich habe keine wirklich gute Antwort, aber lass mich eine Idee rauswerfen. Können Sie das Problem umdrehen und anstatt zu fragen, in welchen Polygonpunkt jeder Punkt passt, stattdessen "welche Punkte sind in jedem Poly?" Vielleicht könntest du deine Shapefiles zum Beispiel in 2000 Grafschaften aufstoßen, dann inkrementell jede Grafschaft nehmen und jeden Punkt überprüfen, um zu sehen, ob es in dieser Grafschaft ist. Wenn sich ein Punkt in einer bestimmten Region befindet, markieren Sie ihn und entfernen ihn das nächste Mal aus Ihrer Suche.

Entlang der gleichen Linien können Sie das Shapefile in 4 Regionen aufteilen. Dann können Sie eine einzelne Region und alle Ihre Punkte in den Speicher einfügen. Dann iterieren Sie 4 mal durch die Daten.

Eine andere Idee wäre, ein GIS-Tool zu verwenden, um die Auflösung (Anzahl der Knoten und Vektoren) des Shapefiles zu verringern. Das hängt natürlich davon ab, wie wichtig die Genauigkeit für Ihren Anwendungsfall ist.

Answer 4

Ich war interessiert, das zu sehen und fragte mich, ob Sie in dieser Hinsicht Fortschritte gemacht hätten. Da Sie die Frage gestellt haben, stelle ich mir vor, dass Ihre Computerhardware und die Software, die Sie für diese relativ einfache Operation verwenden können, sich etwas verbessert haben, so dass die Lösung (wenn sie noch benötigt wird) recht einfach sein kann verarbeiten Sie eine Million Punkte. Vielleicht möchten Sie etwas wie:

# Load relevant libraries 
library(sp) 
library(maptools) 
library(spatstat) 

# Read shapefile. Hopefully you have a .prj file with your .shp file 
# otherwise you need to set the proj4string argument. Don't inlcude 
# the .shp extension in the filename. I also assume that this will 
# create a SpatialPolygonsDataFrame with the "Severity" attribute 
# attached (from your .dbf file). 
myshapefile <- readShapePoly("myshapefile_without_extension",  proj4string=CRS("+proj=latlong +datum=WGS84")) 


# Read your location data in. Here I assume your data has two columns X and Y giving  locations 
# Remeber that your points need to be in the same projection as your shapefile. If they aren't 
# you should look into using spTransform() on your shapefile first. 
mylocs.df <- read.table(mypoints.csv, sep=",", h=TRUE) 

# Coerce X and Y coordinates to a spatial object and set projection to be the same as 
# your shapefile (WARNING: only do this if you know your points and shapefile are in 
# the same format). 
mylocs.sp <- SpatialPoints(cbind(mylocs.df$X,mylocs.df$Y),  proj4string=CRS(proj4string(myshapefile)) 

# Use over() to return a dataframe equal to nrows of your mylocs.df 
# with each row corresponding to a point with the attributes from the 
# poylgon in which it fell. 
severity.df <- over(mylocs.sp, myshapefile) 

Hoffentlich kann dieses Framework Ihnen geben, was Sie wollen. Ob Sie es mit dem Computer/RAM, den Sie jetzt haben, tun können, ist eine andere Sache!

Answer 5

Ich würde fastshp Paket einen Versuch geben. In meinen kursorischen Tests schlägt es deutlich other methods für reading Shapefiles. Und es hat sich spezialisiert inside Funktion, die meine gut Ihren Bedürfnissen entsprechen.

-Code sollte irgendwie ähnlich sein:

shp <- read.shp(YOUR_SHP_FILE, format="polygon")) 
inside(shp, x, y) 

wo x und y-Koordinaten sind.

Wenn das nicht funktioniert, würde ich für die von @ Spacedman erwähnte PostGIS-Lösung gehen.

Answer 6

Um meine eigene Frage zu beantworten ... und danke an alle für ihre Hilfe - Die endgültige Lösung war die Verwendung von gdal aus Python, das relativ leicht an Raster und Shapedateien angepasst werden konnte. Einige der Raster liefen auf etwa 40 GB, und einige der Shape-Dateien überstiegen 8 GB - also gab es keine Möglichkeit, sie in den Speicher der Maschinen einzupassen, die wir damals hatten (Jetzt habe ich Zugriff auf eine Maschine mit 128 GB RAM - aber ich bin auf neue Weiden gezogen!). Der Python/Gdal-Code war in der Lage, 27 Millionen Punkte in 1 Minute bis 40 Minuten zu markieren, abhängig von der Größe der Polygone in den Shapefiles - wenn es viele kleine Polygone gab, war es erstaunlich schnell - wenn es massive (250k Punkte) gab) Polygone in den Shapefiles war es erstaunlich langsam! Zum Vergleich: Früher haben wir es auf einer Oracle-Spatial-Datenbank ausgeführt, und es dauerte etwa 24 Stunden, um die 27 Millionen Punkte zu markieren, oder das Rastern und Markieren würde ungefähr eine Stunde dauern. Wie Spacedman vorgeschlagen habe, habe ich versucht, postgis auf meiner Maschine mit einer ssd zu verwenden, aber die Bearbeitungszeit war ziemlich viel langsamer als die Verwendung von Python/gdal, da die Shapefiles nicht in die Postgis geladen werden mussten. Um es zusammenzufassen, den schnellsten Weg, dies zu tun, wurde Python/gdal mit:

• kopieren Sie die Shape-Dateien und die x, y csv
• ändern Sie die Konfigurationsdatei für den Python-Skript ssd die Dateien zu sagen, wo waren und welche Schicht gegen
• Lauf ein paar Schichten parallel zu markieren - wie es cpu eher begrenzt war als i/o