Tell me more ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.
up vote 3 down vote favorite
1

Following on from an answer about intersecting polygons with lines to chop up a polygon into smaller polygon units (in QGIS), I wanted to try the same thing in R. However, I cannot seem to find a method that works!

over() doesn't have a method for polygons intersecting with lines; I found gIntersection() from rgeos but it fails:

require(sp)
require(rgeos)

poly <- readShapePoly("polygon.shp")
lines <- readShapeLines("lines.shp")

chopped <- gIntersection(poly, lines)

Giving:

Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, "rgeos_intersection") : 
  UnsupportedOperationException: GeometryGraph::add(Geometry *): unknown geometry type: N4geos4geom18GeometryCollectionE

Update: Here's a link to the files in question.

Update 2: PaulG notes that it works and after updating rgeos and R I got rid of the error above. Thanks PaulG ...

However, gIntersection results in a SpatialLines object no matter whether I put in (poly, lines) or (lines, poly) - whereas the operation I did in QGIS (or Arc, back in the bad ol' days) will divide a polygon with the lines and result in a polygon object, not line.

So, how do I chop up my polygon with the lines and get polygons out?

share|improve this question
I think it would be easier if you could post a copy of the shape files in question so we could run our own tests. Btw, is maptools loaded, too? – R.K. Apr 21 '12 at 16:12
I have found rgeos::gIntersection() to also be quite fussy about the nature of the Spatial* object you pass it. I have also solved very cryptic error messages from this function by using the byid=TRUE parameter. Not sure if this will work in your case. – digitalmaps Apr 22 '12 at 1:38
@R.K. - yes, maptools also loaded; have updated with link to files. – Simbamangu Apr 22 '12 at 11:28
@PaulG: tried it with byid=TRUE also with no success. – Simbamangu Apr 22 '12 at 11:29
It likely comes down to the format of your Spatial* object. Try producing a SpatialLines and a SpatialPolygons object from first principles and compare what you have there to the object you load using readShapePoly(). It could be something like IDs are not continuous, or hole attributes. – digitalmaps Apr 22 '12 at 14:16
show 4 more comments

1 Answer

up vote 1 down vote

You could try using RSAGA. I'm not too familiar with it myself, but the command would be something like:

rsaga.geoprocessor("libshapes_polygons", "Polygon-Line Intersection", list(POLYGONS="polygonshape.shp",LINES="lineshape.shp",INTERSECT="result.shp"))
share|improve this answer
This looks like it might help, but it's really troublesome to install (lots of dependencies and downloads); will try it out though. – Simbamangu Apr 30 '12 at 9:52
Do you have to do the intersection in R? You could just do the intersection in saga and import again afterwards. Anyway, only 2 downloads should be needed for RSAGA: install.packages("RSAGA") and then saga itself. – johanvdw Apr 30 '12 at 9:57
I have other options for doing the intersection (QGIS) but would like to be able to have this available in R. Installing 'saga itself' is the troublesome bit I mention above - OSX has no easy binary install! – Simbamangu May 15 '12 at 5:27

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.