2

I have two polygons shapefiles and I'd like to clip one by the other. I search on google but I could find only clipping by a bounding box or clipping points by polygons, and that's not what I need. I also find something in other programming languague, except in R (http://rosettacode.org/wiki/Sutherland-Hodgman_polygon_clipping#Python). Could you help me?

Thanks Tiago

Tiago
  • 143
  • 1
  • 2
  • 7
  • 3
    Try `rgeos::gIntersection()` as demo'd in (for example) the accepted answer to [this question](http://stackoverflow.com/questions/13982773/crop-for-spatialpolygonsdataframe/13986029#13986029). Alternatively, `raster::intersect()` should also do the trick. – Josh O'Brien Aug 26 '15 at 18:06
  • @JoshO'Brien Thanks for your answer. It helps me to get almost exactly what I want. I still got only one small problem. When I do the `rgeos::gIntersection()` I lost all the lines stored in the attribute table. Is there a way that I can keep it corresponding to what was intersected? – Tiago Aug 27 '15 at 19:02
  • 4
    Perhaps try `raster::intersect()`, which was expressly designed to keep associated data.frame attributes. (And the only reason I say "perhaps" is that occasionally, it doesn't work as well as `rgeos::gIntersection()`, as for example, with the data in the answer I linked to above. But it usually does work quite nicely.) Hope that helps! – Josh O'Brien Aug 27 '15 at 19:16
  • 'Clip' is tricky--in ArcGIS it's one thing, but in R it's two: one for using polygons to subset points, and one for using polygons to subset polygons. – Mox Jan 19 '20 at 17:14

4 Answers4

2

Although, the question is very old, providing a good answer might help anyone in the future landing on this page.

I think what you're trying to do is straight forward. To illustrate, lets assume i'm interested in eastern coastline of Saudi Arabia (SA) and i have a shape file that has the east and west coast of SA [SA[1] and another shapefile of the Gulf Gulf (prominent waterbody on the east coast of SA). We need the SF package to crop the two shape files

Load the SF package 

library(sf)

Then load both shapefiles 

    ksa <- st_read("cropping_shape_file/ksa/saudi_arabia_coastline.shp")
gulf <- st_read("cropping_shape_file/gulf/ihoPolygon.shp")%>%
  st_transform(st_crs(ksa)) # The st_transform added to this file ensure that both files have same CRS ssstem otherwise it will be impossible to crop.

you can also check if their CRS is same 

st_crs(ksa)==st_crs(gulf)

The sf_read outputs shapefile information into three fields but we're interested in just the geometry

you can 

cropped<-st_crop( ksa$geometry, gulf$geometry)
plot(cropped)

You should have an image like cropped

i'm providng the two shape files here https://mega.nz/file/NnxHXajI#KOK_86gwVQGfjNs7X2HOAtvG2NOkIwcU5EhzQU5X6tg, note that a shape file has four component that needs to be kept together(.shp, .prj, .shx, .dbf )

Hammao
  • 801
  • 1
  • 9
  • 28
  • here you are only keeping the `geometry` column. What if I wanted to keep all the columns? Answer: Just get rid of `...$geometry`. – M-- Jan 16 '23 at 20:44
1

Another way for clipping is to make an selection like "polygontoclip"["templatepolygon", ], found this for points (http://robinlovelace.net/r/2014/07/29/clipping-with-r.html), but also works with polygons.

Debroize
  • 29
  • 2
  • 2
    'Clip' is tricky--in ArcGIS it's one thing, but in R it's two: one for using polygons to subset points, and one for using polygons to subset polygons. – Mox Jan 19 '20 at 17:15
1

This doesn't answer your question, but since you mention clipping by bounding box, so this post comes up in search strings:

From r-bloggers: Clipping by a bounding box

setwd("../")
library(rgdal)  
zones <- readOGR("data", "london_sport")

make a bounding box:

b <- bbox(zones)
b[1, ] <- (b[1, ] - mean(b[1, ])) * 0.5 + mean(b[1, ])
b[2, ] <- (b[2, ] - mean(b[2, ])) * 0.5 + mean(b[2, ])
b <- bbox(t(b))
plot(zones, xlim = b[1, ], ylim = b[2, ])

Using a custom function:

library(raster)
library(rgeos)
## rgeos version: 0.3-5, (SVN revision 447)
##  GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921 
##  Polygon checking: TRUE 
gClip <- function(shp, bb){
  if(class(bb) == "matrix") b_poly <- as(extent(as.vector(t(bb))), "SpatialPolygons")
  else b_poly <- as(extent(bb), "SpatialPolygons")
  gIntersection(shp, b_poly, byid = T)
}
zones_clipped <- gClip(zones, b)
## Warning: spgeom1 and  spgeom2 have different proj4 strings
plot(zones_clipped)

Note that due to the if statements in gClip’s body, it can handle almost any spatial data input, and still work.

westminster <- zones[grep("West", zones$name),]
zones_clipped_w <- gClip(zones, westminster)
## Warning: spgeom1 and spgeom2 have different proj4 strings
plot(zones_clipped_w); plot(westminster, col = "red", add = T)
Mox
  • 511
  • 5
  • 15
0

Clip a shapefile by another shapefile

library(raster)
# library(rgdal) # if require
# library(rgeos) # if require

luxembourg <- shapefile(system.file("external/lux.shp", package="raster")) ### load a shape from raster package
plot(luxembourg) 

another_shape <- as(extent(6, 6.4, 49.75, 50), 'SpatialPolygons') ### draw a simple shape / load your own
plot(another_shape, add=T, border="red", lwd=3)

luxembourg_clip <- crop(luxembourg, another_shape) ### crop (SpatialPolygon) luxembourg with another_shape
plot(luxembourg_clip, add=T, col=4) ### plot on 
plot(luxembourg_clip, add=F, col=3) ### just plot

See some images: Map and Clip

Source: https://rdrr.io/cran/raster/man/crop.html

Yuri Gelsleichter
  • 311
  • 2
  • 5