8

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?

Community
  • 1
Simbamangu
  • 14,773
  • 6
  • 59
  • 93
  • 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 01: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
  • @Simbamangu I just used your files to run the above script, and I found that the line should be changed to: poly <- readShapePoly("poly.shp"). There are no files polygon.* in the zip you provided. It then worked for me, performing the intersection as expected. Could the error be as simple as loading the wrong shape file? (I am using rgeos package version 0.2-1 under R 2.14.2) – digitalmaps Apr 22 '12 at 20:11
  • @PaulG - no, I just renamed it incorrectly when I archived the files; it's polygon.shp on mine, so that's not it. I'm going to try your suggestion to try building some shapes directly (these were made by QGIS) or use another source, see what I get. – Simbamangu Apr 23 '12 at 02:59
  • @PaulG - after carefully re-reading your comments, discovered I was on rgeos 0.1-5 and updated everything (R to 2.15 / rgeos to 0.2-5). gIntersection works, but not as I'd like (see edit above). – Simbamangu Apr 29 '12 at 17:56
  • @Simbamangu it doesn't seem that gIntersection() will really do what you want. If your lines, forming a grid, were polygonized, it might work, but this is probably not what you want. I suggest asking this either on StackOverflow where more spatial R folks hang out, or on the R-sig-geo list where you'll reach the architects of rgeos and related tools. – digitalmaps Apr 30 '12 at 01:34
  • I usually intersect polygons and lines with qgis and it's very fast. Using R, today i'm trying to intercept a sp object with 24500 lines with another sp object, a gridpolygon. The process started at 19:12 LT and now it's 20:08 and it's still dosent finish. This are the commands from here. library(rgdal) library(raster) library(rgeos) library(dismo) ingrid <- intersect(mylines, gridpolygon) I'm not sure if it will help, i'm looking for a similar command in R, but it takes too much time. You could try [@johanvdw answer](http://gis.stackexchange.com/a/23973/1 – Sergio Sep 18 '15 at 23:19

2 Answers2

4

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.
Sergio
  • 948
  • 9
  • 15
  • 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
4

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"))
GIS-Jonathan
  • 6,775
  • 5
  • 29
  • 58
johanvdw
  • 6,207
  • 27
  • 42
  • 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 09: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 09: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 05:27