punti r nei poligoni

Question

Ho un milione di punti e un file di forma grande -8GB, che è troppo grande per essere caricato in memoria in R sul mio sistema. Il file shape è single-layer quindi un dato x, colpirà al massimo un poligono, purché non si trovi esattamente su un confine! Ogni poligono è etichettato con un severity - ad es. 1, 2, 3. Sto usando R su una macchina Ubuntu 64-bit con 12 GB di RAM.

Qual è il modo più semplice per essere in grado di "tag" il frame di dati al poligono severity modo che ho un data.frame con una colonna in più, vale a dire x, y, severity?

Answer 1

Solo perché tutto quello che hai è un martello, non significa che ogni problema è un chiodo.

Carica i tuoi dati in PostGIS, crea un indice spaziale per i poligoni e fai un singolo overlay spaziale SQL. Esportare i risultati di nuovo a R.

A proposito, dire che lo shapefile è 8 Gb non è una informazione molto utile. Gli shapefile sono costituiti da almeno tre file, il .shp che è la geometria, il .dbf che è il database e il .shx che collega i due. Se il tuo .dbf è 8Gb, puoi facilmente leggere le forme stesse sostituendole con un diverso .dbf. Anche se il file .shp è di 8 GB, potrebbe essere solo tre poligoni, nel qual caso potrebbe essere facile semplificarli. Quanti poligoni hai e quanto è grande la parte .shp dello shapefile?

Answer 2

Penso che dovresti preprocessare i dati e creare una struttura che elenca i poligoni possibili per le aree rettangolari in una griglia. In questo modo, puoi ridurre i poligoni che dovrai controllare rispetto ai punti e la struttura aggiuntiva si adatterà alla memoria poiché ha solo puntatori ai poligoni.

Ecco un'immagine per illustrare:.

Si vuole verificare quale poligono il punto giallo è in voi di solito controlla contro tutti i poligoni, ma con l'ottimizzazione di rete (le linee arancioni, non hai disegnato l'intera griglia, solo uno dei suoi campi) devi solo controllare i poligoni pieni come sono tutti all'interno o parzialmente all'interno del campo della griglia.

Un modo simile sarebbe non memorizzare tutti i dati poligono in memoria, ma solo i poligoni che delimitano caselle che richiedono solo 2 invece di 3 coppie X/Y per ogni poligono (e un puntatore aggiuntivo ai dati del poligono reale), ma ciò non risparmia tanto spazio quanto il primo suggerimento.

Answer 3

Non ho una risposta davvero buona, ma lasciami dare un'idea. Puoi capovolgere il problema e invece di chiedere quale poli ogni punto si adatta, invece come "quali punti sono in ogni poli?" Magari potresti distruggere il tuo shapefile in, per esempio, 2000 contee, poi incrementare progressivamente ogni contea e controllare ogni punto per vedere se è in quella contea. Se un punto si trova in una determinata contea, allora lo tagghi e lo toglierà dalla tua ricerca la prossima volta.

Lungo le stesse linee, è possibile interrompere lo shapefile in 4 regioni. Quindi puoi adattare una singola regione e tutti i tuoi punti in memoria. Quindi basta ripetere 4 volte attraverso i dati.

Un'altra idea sarebbe quella di utilizzare uno strumento GIS per ridurre la risoluzione (numero di nodi e vettori) dello shapefile. Ciò dipende, naturalmente, da quanto sia importante la precisione per il tuo caso d'uso.

Answer 4

Mi interessava vedere questo e mi chiedevo se avessi fatto qualche progresso su questo fronte. Dal momento che hai posto la domanda, immagino che il tuo hardware e il software che puoi utilizzare per eseguire questa operazione relativamente semplice siano migliorati al punto che la soluzione (se ancora necessaria!) Potrebbe essere abbastanza semplice, anche se potrebbe richiedere molto tempo per elaborare un milione di punti. Potresti provare qualcosa del tipo:

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

Speriamo che questo quadro possa darti quello che vuoi. Se puoi farlo o no con il computer/RAM che hai a disposizione ora è un'altra questione!

Answer 5

Mi piacerebbe dare fastshp pacchetto una prova. Nei miei test rapidi, batte in modo significativo lo other methods per gli shapefile reading. Ed è specializzata nella funzione inside che soddisfa le tue esigenze.

codice dovrebbe essere in qualche modo simile a:

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

dove x ed y sono coordinate.

Se ciò non funziona, opterei per la soluzione PostGIS menzionata da @Spacedman.

Answer 6

Per rispondere alla mia domanda ... e in ringraziamento a tutti per il loro aiuto - La soluzione finale adottata è stata l'utilizzo di gdal da python che è stato relativamente facilmente adattato sia ai raster che ai file di forma. Alcuni di questi raster hanno funzionato a circa 40 GB e alcuni file di forma hanno superato gli 8 GB, quindi non c'era modo di adattarsi alla memoria su nessuna delle macchine che avevamo in quel momento (ora ho accesso a una macchina con 128 GB di RAM - ma mi sono trasferito su nuovi pascoli!). Il codice python/gdal era in grado di taggare 27 milioni di punti in da 1 minuto a 40 minuti a seconda delle dimensioni dei poligoni negli shapefile - se ci fossero molti piccoli poligoni era incredibilmente veloce - se ce ne fossero di massicci (250k punti) poligoni negli shapefile era incredibilmente lento! Tuttavia, per confrontarlo, lo abbiamo usato in precedenza su un database spaziale Oracle e ci sarebbero voluti circa 24 ore + per taggare i 27 milioni di punti, o la rasterizzazione e l'etichettatura richiederebbero circa un'ora. Come suggerito da Spacedman, ho provato a usare postgis sulla mia macchina con uno ssd, ma il tempo di turn-around è stato un po 'più lento rispetto all'uso di python/gdal poiché la soluzione finale non richiedeva il caricamento degli shapefile in postgis. Quindi, per riassumere, il modo più rapido per fare questo stava usando Python/gdal:

• copiare i file di forma e la x, y CSV per SSD
• modificare il file di configurazione per lo script python a dire dove i file erano e che livello si intende la modifica contro
• corsa un paio di strati in parallelo - come è stato CPU limitato, piuttosto che i/o limitato