2012-03-22 13 views
6

Voglio sottoporre a serie uno shapefile (il file .shp e i file associati sono here) in un altro delimitato da un insieme di coordinate, diciamo tra long [80,90] e lats [20,30], quindi scrivi questo come un altro shapefile. Se uso il pacchetto maptools:R/GIS: Come suddividere uno shapefile con un bounding box lat-long?

df = readShapeLines("/path/asia_rivers.shp")

e poi guardare la struttura del file con as.data.frame(df), non riesco a trovare alcun modo ovvio di sottoinsiemi per coordinate. I posso utilizzare il pacchetto PBSmapping al sottoinsieme:

df = importShapefile("/path/asia_rivers.shp") 
df_sub = subset(df, X>=80 & X<=90 & Y >=20 & Y <=30) 

ma poi non riesco a essere in grado di forzare questo in un frame di dati che possono essere esportati tramite writeSpatialShape() in maptools. Continuo a ricevere questo errore: Error in PolySet2SpatialLines(df_sub) : unknown coordinate reference system. Sicuramente mi manca qualcosa di molto semplice e ci dovrebbe essere un modo semplice per inserire i geo-dati per geo-coordinate?

risposta

6

si potrebbe provare la seguente:

library(rgeos) 
rivers <- readWKT("MULTILINESTRING((15 5, 1 20, 200 25), (-5 -8,-10 -8,-15 -4), (0 10,100 5,20 230))") 
bbx <- readWKT("POLYGON((0 40, 20 40, 20 0, 0 0, 0 40))") 

rivers.cut <- gIntersection(rivers, bbx) 

plot(rivers, col="grey") 
plot(bbx, add=T, lty=2) 
plot(rivers.cut, add=T, col="blue") 
+0

Il problema non è nel tracciare, ma in subsetting in un formato che può essere esportato come shapefile. – user702432

+0

Hai provato a leggere il tuo shapefile, intersecandolo con un riquadro di delimitazione appropriato utilizzando gIntersection e quindi esportalo, ad esempio con maptools :: writeShapeLines? L'ho solo tracciato per mostrarti che le linee sono state ritagliate con dati fittizi, dal momento che il file di forma che hai fornito non è valido (uno shapefile consiste in più di un semplice file * .shp; http://en.wikipedia.org/ wiki/Shapefile). – johannes

+0

Ciao jmsigner ... Mi dispiace, mio ​​male. Ho cambiato il link alla cartella zippata che contiene i file .dbf, .lyr, .prj, .sbn, .sbx, .shp e .shx originali. Grazie. – user702432

2

So che questo è stato risposto, ma penso che si può fare esattamente quello che si desidera utilizzare PBSmapping. PBSmapping ha una funzione per ritagliare Polysets (per i dati poligonali e le linee) in modo da poter provare:

df <- importShapefile("/path/asia_rivers.shp") 
df_sub <- clipLines(df, xlim = c(80 , 90) , ylim = c(20 , 30), keepExtra = TRUE) 
dfSL <- PolySet2SpatialLines(df_sub) 

L'extra mastio consente di mantenere le colonne non standard quando si esegue il ritaglio (suppongo per i dati degli attributi).

1

Un altro modo:

library(raster) 

s <- shapefile("/path/asia_rivers.shp") 

sub <- crop(s, extent(80, 90, 20, 30)) 

shapefile(sub, 'cropped_rivers.shp') 
Problemi correlati