Crop for SpatialPolygonsDataFrame Crop for SpatialPolygonsDataFrame r r

Crop for SpatialPolygonsDataFrame


Streamlined method added 2014-10-9:

raster::crop() can be used to crop Spatial* (as well as Raster*) objects.

For example, here's how you might use it to crop a SpatialPolygons* object:

## Load raster package and an example SpatialPolygonsDataFramelibrary(raster) data("wrld_simpl", package="maptools")## Crop to the desired extent, then plotout <- crop(wrld_simpl, extent(130, 180, 40, 70))plot(out, col="khaki", bg="azure2")

Original (and still functional) answer:

The rgeos function gIntersection() makes this pretty straightforward.

Using mnel's nifty example as a jumping off point:

library(maptools)library(raster)   ## To convert an "Extent" object to a "SpatialPolygons" object.library(rgeos)data(wrld_simpl)## Create the clipping polygonCP <- as(extent(130, 180, 40, 70), "SpatialPolygons")proj4string(CP) <- CRS(proj4string(wrld_simpl))## Clip the mapout <- gIntersection(wrld_simpl, CP, byid=TRUE)## Plot the outputplot(out, col="khaki", bg="azure2")

enter image description here


Here is an example of how to do this with rgeos using the world map as an example

This comes from Roger Bivand on R-sig-Geo mailing list. Roger is one of the authors of the sp package.

Using the world map as an example

library(maptools)data(wrld_simpl)# interested in the arealy bounded by the following rectangle# rect(130, 40, 180, 70)library(rgeos)# create  a polygon that defines the boundarybnds <- cbind(x=c(130, 130, 180, 180, 130), y=c(40, 70, 70, 40, 40))# convert to a spatial polygons object with the same CRSSP <- SpatialPolygons(list(Polygons(list(Polygon(bnds)), "1")),proj4string=CRS(proj4string(wrld_simpl)))# find the intersection with the original SPDFgI <- gIntersects(wrld_simpl, SP, byid=TRUE)# create the new spatial polygons object.out <- vector(mode="list", length=length(which(gI)))ii <- 1for (i in seq(along=gI)) if (gI[i]) {  out[[ii]] <- gIntersection(wrld_simpl[i,], SP)  row.names(out[[ii]]) <- row.names(wrld_simpl)[i]; ii <- ii+1}# use rbind.SpatialPolygons method to combine into a new object.out1 <- do.call("rbind", out)# look here is Eastern Russia and a bit of Japan and China.plot(out1, col = "khaki", bg = "azure2")

enter image description here


You cannot use crop on sp polygon objects. You will need to create a polygon representing the bbox coordinates of dat2 and then can use gIntersects in the rgeos library.

Edit: This comment was in relation to the version available in 2012 and this is no longer the case.