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