[GIS] How to intersect lines and polygons in R

geoprocessingr

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?

Best Answer

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

Related Solutions

[GIS] Intersection between two SpatialPolygonsDataFrame using gIntersection function in R

I cannot as yet offer a complete solution, but the immediate issue is that the parcel object has invalid geometry. If you can fix the geometry, we should be able to get the intersections to work. As an example of problematic geometry, this expression is fine:

gIsValid(hr_kud, reason = TRUE)

But this one will consistently crash R:

gIsValid(parcel, reason = TRUE)

So, one or more of the 20,069 polygons in parcel appears to be broken in some way. I don't know how to validate geometry within R but Quantum GIS does have a tool for this. I assume that you have the original shapefile or other data source from which you created parcel but you can recreate this by using the writeOGR function from the rgdal library like so:

library(rgeos)
library(rgdal)

setwd("your_directory_name_here")

load("Data.RData")
class(parcel)
class(hr_kud)
object.size(hr_kud)
object.size(parcel)
parcel@proj4string
hr_kud@proj4string

gIsValid(hr_kud, reason = TRUE)
#gIsValid(parcel, reason = TRUE) # warning! this will crash R

writeOGR(parcel, ".", driver = "ESRI Shapefile", layer = "poly")

This should produce 4 files:

poly.shp
poly.shx
poly.dbf
poly.prj

In QGIS, first load your shapefile by going to the 'Layer' menu and clicking on 'Add Vector Layer..." then select the poly.shp file or the folder in which the file resides. Then go to the 'Vector' menu, click on the 'Geometry Tools' menu item then on the 'Check geometry validity' item. Then click OK on the next dialog box. With more than 20k polys this took quite some time (about half an hour I think), but the result was this: [EDIT: I upgraded to QGIS 1.8 Lisboa with two results - a massive increase in validation speed and the discovery of other errors, as per the revised screenshot below.]

QGIS geometry check result

So, it looks like you have a number of issues there. Before you close the dialog box, click on the error in the dialog box list and QGIS will show you the offending polygon. I'm not experienced in this kind of GIS but the tips in this tutorial might help. Google QGIS and "self intersecting" for more ideas.

[GIS] How to perform a true GIS clip of polygons layer using a polygon layer in R

The hint provided by @mdsummer of using byid=TRUE works accurately.

See the reproducible example, below:

library(rgeos)
library(sp)

#Create SpatialPlygons objects
polygon1 <- readWKT("POLYGON((-190 -50, -200 -10, -110 20, -190 -50))") #polygon 1
polygon2 <- readWKT("POLYGON((-180 -20, -140 55, 10 0, -140 -60, -180 -20))") #polygon 2

par(mfrow = c(1,2)) #in separate windows
plot(polygon1, main = "Polygon1") #window 1
plot(polygon2, main = "Polygon2") #window 2

polygon_set <- readWKT(paste("POLYGON((-180 -20, -140 55, 10 0, -140 -60, -180 -20),",
                     "(-190 -50, -200 -10, -110 20, -190 -50))"))

par(mfrow = c(1,1)) #now, simultaneously
plot(polygon_set, main = "Polygon1 & Polygon2")

clip <- gIntersection(polygon1, polygon2, byid = TRUE, drop_lower_td = TRUE) #clip polygon 2 with polygon 1
plot(clip, col = "lightblue")

enter image description here

GT <- GridTopology(c(-175, -85), c(10, 10), c(36, 18))
gr <- as(as(SpatialGrid(GT), "SpatialPixels"), "SpatialPolygons")
plot(gr)

enter image description here

clip2 <- gIntersection(clip, gr, byid = TRUE, drop_lower_td = TRUE)
plot(clip2, col = "pink")

enter image description here

Edit:

If one wants to intersect while keeping the attribute data, then one option is the raster::intersect function. See jeb's answer.

Best Answer

Related Solutions

[GIS] Intersection between two SpatialPolygonsDataFrame using gIntersection function in R

[GIS] How to perform a true GIS clip of polygons layer using a polygon layer in R

Related Question