7

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?

|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
2

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"))
|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
2

I've been facing the same problem several times. I'will focus on two methods that works for me

1) The intersect function from raster package. It returns a SpatialLinesDataframe including the id of each polygon.

library(raster)
newLines <- raster::intersect(SpatialLinesDataframe, SpatialPolygonDataframe)

2) Calling QGIS from R with the package RQGIS, as showed here I think that this method is faster than several other, but you have to take care with the qgis version. Here is the github of the creator. Also, with this package you can use any of the geoalgoritm from qgis.

library(RQGIS)    
qgis_env <- set_env("PATH")
params <- get_args_man(alg = "qgis:intersection", 
                           qgis_env = qgis_env)
params$INPUT <- line_1
params$INPUT2 <- poly_1
out <- run_qgis(alg = "qgis:intersection",
                params = params,
                load_output = params$OUTPUT,
                qgis_env = qgis_env)

EDIT in 2018: I was thinking and you polygon with another polygon, to get smaller polygons: Take a look of this link for the Manual of the new sf R package (Pebesma, 2018).

  • Edzer Pebesma (2018). sf: Simple Features for R. R package version 0.6-3. https://CRAN.R-project.org/package=sf
  • Pebesma, E., 2018. Simple Features for R: Standardized Support for Spatial1995 Vector Data. The R Journal.
|improve this answer
  • Interesting ... however, (1) results in lines, which isn't what I'm after (want to chop up a polygon into smaller polygons); (2) fails with "Error in sys.excepthook", will keep trying though. – Simbamangu Oct 18 '16 at 17:21
  • Maybe an intersection between your original polygon and another polygon formed with your lines would chop your original polygon – Sergio Oct 18 '16 at 17:40
  • Simbamangu, (1) does what you want. The info is stored in newLines$id (where id is the identifier of your polygon features). So if you subset yourPolygon[newLines$id, ] you get exactly the features of the polygon that are touched by the spatialLine. – Fitzroy Hogsflesh Jun 13 '18 at 22:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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