16

I would like to do a true GIS Clip in R of soils polygons using a series of single boundary polygons, but I cannot find an R function to properly do it. It should work just like the clip function in ESRI's ArcMap. I've tried the over method in sp package but it doesn't seem to work for polys over polys.

One suggestion was to use the gIntersection in rgeos package as a clip using the following code:

#------------------------------------
library(rgeos)
library(maptools)

#Read layers as SpatialPolygonsDataFrame (both the same Albers projection)
Soils_poly = readShapePoly("Soils_polygons")  #Note - Has 400 polygons
clipper_poly = readShapePoly("clipper_polygon")  #Note - Has 1 polygon

#Try gintersection as clip 
Clipped_polys = gIntersection(Clipper_Tile_poly, Soils_poly)

#-----------------------------------

This takes 5 minutes to run (way too slow) and errors with this:

Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_not_poly, "rgeos_intersection") : TopologyException: no outgoing dirEdge found at -721459.77681285271 2009506.5980877089

I also tried this code to check for overlap:

gIntersects(Clipper_Tile_poly, Soils_poly)

and the result was TRUE. clip function in ESRI ArcMap works fine for this data.

Anyone knows of an R function to properly do a clip on spatial polygons using spatial polygons?

  • Try gIntersection with byid=TRUE (I think it's a problem with dangling topologies for certain clips, and sometimes doing it this way helps), for speed check out gUnarySTRtreeQuery() or gBinarySTRtreeQuery() to identify intersecting bounding boxes of pairs of polygons and only intersect those pairs. There's no easy high-level wrappers for all of this afaik – mdsumner Apr 14 '14 at 22:37
20

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

polygons side-by-side

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

enter image description here

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

  • 1
    Works a dream - unbelievably fast. I'm using this to clip polylines. I wanted to make subsets of the shapefile for British rivers, because it's much faster to work with a smaller subset of data. – CJB Apr 6 '16 at 9:12
  • 1
    @AndreSilva, thought I had, but guess I hadn't! Re: community, I wish the R GIS developer community was a bit larger TBH. I will probably have to fall back on ArcMap for some digitization, and other processes that are quick in the GUI. – Rich Pauloo Sep 9 '17 at 3:05
3

You can also use the raster package raster::intersect(spdf1, spdf2) . It has the advantage of maintaining the attributes in case you have a SpatialPolygonsDataFrame.

library(sp); library(rgeos)

coords1 <- matrix(c(-1.841960, -1.823464, -1.838623, -1.841960, 55.663696,
                55.659178, 55.650841, 55.663696), ncol=2)
coords2 <- matrix(c(-1.822606, -1.816790, -1.832712, -1.822606, 55.657887,
                55.646806, 55.650679, 55.657887), ncol=2)
p1 <- Polygon(coords1)
p2 <- Polygon(coords2)
p1 <- Polygons(list(p1), ID = "p1")
p2 <- Polygons(list(p2), ID = "p2")
myPolys <- SpatialPolygons(list(p1, p2))
spdf1 = SpatialPolygonsDataFrame(myPolys, data.frame(variable1 = c(232,
                                                               242), variable2 = c(235, 464), row.names = c("p1", "p2")))
coords1a <- matrix(c(-1.830219, -1.833753, -1.821154, -1.830219, 55.647353,
                 55.656629, 55.652122, 55.647353), ncol=2)
p1a <- Polygon(coords1a)
p1a <- Polygons(list(p1a), ID = "p1a")
myPolys1 <- SpatialPolygons(list(p1a))
spdf2 = SpatialPolygonsDataFrame(myPolys1, data.frame(variable1 = c(2),
                                                  variable2 = c(3), row.names = c("p1a")))

# works but drop the attributes
#gIntersection(spdf1, spdf2, byid=T)

#better to keep attributes
inter1=raster::intersect(spdf1, spdf2)

plot(spdf1, col="red")
plot(spdf2, col="yellow", add=T)
plot(inter1,col="blue", add=T)

enter image description here

Thanks to this question to pointing that out and for sample code.

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.