7

I want to subset out a shapefile (the .shp and associated files are here) into another one bounded by a set of coordinates, say between longs [80,90] and lats [20,30], and then write this out as another shapefile. If I use the maptools package:

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

and then look at the structure of the file with as.data.frame(df), I can't find any obvious way of subsetting by coordinates. I can use the PBSmapping package to subset:

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

but then I can't seem to be able to coerce this into a SpatialLines data frame which can be exported via writeSpatialShape() in maptools. I keep getting this error: Error in PolySet2SpatialLines(df_sub) : unknown coordinate reference system. Surely I am missing something very basic and there should be an easy way of subsetting geo-data by geo-coordinates?

user702432
  • 11,898
  • 21
  • 55
  • 70

3 Answers3

6

You could try the following:

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")
johannes
  • 14,043
  • 5
  • 40
  • 51
  • The problem is not in plotting, but in subsetting in a format which can be exported as a shapefile. – user702432 Mar 22 '12 at 17:42
  • Did you try to read your shapefile, intersect it with a an appropriate bounding box using gIntersection and then export the it, e.g. with maptools::writeShapeLines? I only plotted it to show you that the lines have been trimmed with dummy data, since the shape file you provided is not valid (a shapefile consists of more than just a *.shp file; http://en.wikipedia.org/wiki/Shapefile). – johannes Mar 22 '12 at 18:03
  • Hi jmsigner... Sorry, my bad. I changed the link to the zipped folder which contains the original .dbf, .lyr, .prj, .sbn, .sbx, .shp, and .shx files. Thanks. – user702432 Mar 23 '12 at 01:48
  • 1
    Hi, given the size of the shapefile, I would preprocess the data in a GIS or with ogr2ogr from gdal. The following just took a few seconds on my computer: ogr2ogr -spat 80 20 90 30 asisa_rivers_cut.shp asia_rivers.shp – johannes Mar 23 '12 at 08:37
  • Oh boy! Looks like this is going to be a long walk :-) Thanks for gIntersection tip, jmsigner. One thing... Could you please tell me why when I do gIntersection(SpatialPolygon, SpatialLinesDataFrame), the result is simply of class SpatialLines (which cannot be exported with writeSpatialShape())? – user702432 Mar 23 '12 at 10:10
2

I know this has been answered, but I think you can do exactly what you want using PBSmapping. PBSmapping has a function to clip Polysets (for polygon and lines data) so you can try:

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

The keep extra allows you to keep non-standard columns when you do the clipping (I assume for attribute data).

Hakim
  • 3,225
  • 5
  • 37
  • 75
Simon O'Hanlon
  • 58,647
  • 14
  • 142
  • 184
2

Another way:

library(raster)

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

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

shapefile(sub, 'cropped_rivers.shp')
Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63