r puntos en polígonos

Question

Tengo un millón de puntos y un archivo de forma grande-8GB-que es demasiado grande para cargar en la memoria en R en mi sistema. El archivo de forma es de una sola capa, por lo que dado x, y se golpeará como máximo un polígono, ¡siempre que no se encuentre exactamente en un límite! Cada polígono está etiquetado con un severity - p. 1, 2, 3. Estoy usando R en una máquina ubuntu de 64 bits con 12GB ram.

¿Cuál es la forma más sencilla de ser capaz de "etiqueta" de la trama de datos al polígono severity de modo que consiga un data.frame con una columna adicional, es decir x, y, severity?

Answer 1

El hecho de que todo lo que tiene es un martillo, no significa que cada problema es un clavo.

Cargue sus datos en PostGIS, cree un índice espacial para sus polígonos y realice una sola superposición espacial de SQL. Exportar los resultados a R.

Por cierto, decir que el shapefile es de 8 Gb no es una información muy útil. Los archivos Shapefiles están formados por al menos tres archivos, el .shp que es la geometría, el .dbf que es la base de datos y el .shx que conecta los dos. Si su .dbf es de 8 Gb, puede leer fácilmente las formas en sí reemplazándolo con un .dbf diferente. Incluso si el .shp es de 8 Gb, podría ser de solo tres polígonos, en cuyo caso podría ser fácil simplificarlos. ¿Cuántos polígonos tienes y qué tan grande es la parte .shp del shapefile?

Answer 2

Creo que debe preprocesar los datos y crear una estructura que enumere posibles polígonos para áreas rectangulares en una cuadrícula. De esta forma, puede reducir los polígonos que tendrá que verificar contra los puntos y la estructura adicional se ajustará a la memoria ya que solo tiene punteros a los polígonos.

Aquí está una imagen para ilustrar:.

Usted desea comprobar qué polígono el punto amarillo está en que se suele comprobar contra todos los polígonos, pero con la optimización de la rejilla (las líneas de color naranja, no dibujó toda la cuadrícula, solo uno de sus campos), solo debe verificar los polígonos rellenos, ya que todos están dentro o parcialmente dentro del campo de la cuadrícula.

Una manera similar sería no almacenar todos los datos de polígono en la memoria, solo los polígonos que delimitan cuadros que solo requerirían 2 en lugar de 3 pares X/Y para cada polígono (y un puntero adicional a los datos de polígono reales), pero esto no ahorra tanto espacio como la primera sugerencia.

Answer 3

No tengo una muy buena respuesta, pero déjame echar una idea. ¿Puedes darle la vuelta al problema y en lugar de preguntar en qué punto de cada punto encaja, en lugar de "qué puntos están en cada poli?" Tal vez podría reventar su shapefile en, por ejemplo, 2000 condados, luego tomar incrementalmente cada condado y verificar cada punto para ver si se encuentra en ese condado. Si un punto está en un condado dado, entonces lo etiqueta y lo saca de su búsqueda la próxima vez.

En la misma línea, puede dividir el shapefile en 4 regiones. Entonces puede ajustar una sola región más todos sus puntos en la memoria. Luego solo itera 4 veces a través de los datos.

Otra idea sería utilizar una herramienta GIS para reducir la resolución (número de nodos y vectores) del shapefile. Eso depende, por supuesto, de la importancia de la precisión para su caso de uso.

Answer 4

Me interesó ver esto y me pregunté si habían hecho algún progreso en este frente. Desde que hizo la pregunta, imagino que el hardware de su computadora y el software que puede usar para hacer esta operación relativamente simple han mejorado hasta cierto punto, hasta el punto en que la solución (si aún es necesaria) puede ser bastante simple, aunque puede tomar mucho tiempo. procesar un millón de puntos. Quizás desee probar algo como:

# 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) 

Espero que ese marco pueda darle lo que desea. Si usted puede o no puede hacerlo con la computadora/RAM que tiene disponible ahora es otro asunto.

Answer 5

Daría fastshp paquete una prueba. En mis pruebas superficiales, supera significativamente a other methods para reading shapefiles. Y se ha especializado la función inside que mi bien se adapta a sus necesidades.

Código debe ser de alguna manera similar a:

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

donde x e y son las coordenadas.

Si eso no funciona, buscaría la solución PostGIS mencionada por @Spacedman.

Answer 6

Para responder a mi propia pregunta ... y en agradecimiento a todos por su ayuda - La solución final adoptada fue usar gdal de python, que se adaptó con relativa facilidad tanto a los rásteres como a los archivos de formas. Algunos rásteres corrieron a unos 40 gb y algunos de los archivos de forma superaron los 8 gb, por lo que no había manera de que cupieran en la memoria de ninguna de las máquinas que teníamos en ese momento (ahora tengo acceso a una máquina con 128 gb de ram - ¡pero me he mudado a nuevos pastos!). El código python/gdal logró etiquetar 27 millones de puntos en un lapso de 1 minuto a 40 minutos, dependiendo de los tamaños de los polígonos en los shapefiles (si había muchos polígonos pequeños, era increíblemente rápido) si había algunos puntos masivos (250k puntos).) polígonos en los shapefiles fue increíblemente lento! Sin embargo, para comparar, solíamos ejecutarlo previamente en una base de datos espacial de Oracle y tomaría aproximadamente 24 horas + marcar los 27 millones de puntos, o rastrillar y etiquetar demoraría aproximadamente una hora. Como Spacedman sugirió, intenté usar postgis en mi máquina con una ssd, pero el tiempo de respuesta fue un poco más lento que con python/gdal, ya que la solución final no requería que los shapefiles se cargaran en Postgis. Entonces, para resumir, la forma más rápida de hacer esto fue el uso de Python/gdal:

• copiar los archivos de forma y la x, y CSV a SSD
• modificar el archivo de configuración para el script en Python que decir donde están los archivos eran y qué capa de etiqueta contra el
• plazo de un par de capas en paralelo - ya que se cpu limitada en lugar de i/o limitado