Given a Spatial object in R, how do I clip all of its elements to lie within a bounding box?

There are two things I'd like to do (ideally I'd know how to do both, but either is an acceptable solution to my current problem--restricting a polygon shapefile to the continental U.S.).

  1. Drop each element not fully within the bounding box. This seems like bbox()<- would be the logical way, but no such method exists.

  2. Do a true clip operation, such that non-infinitesimal elements (e.g. polygons, lines) are cut off at the boundary. sp::bbox lacks an assignment method, so the only way I've come up with would be to use over or gContains/gCrosses in conjunction with a SpatialPolygons object containing a box with the new bounding box's coordinates. Then when clipping a polygon object, you'd have to figure out which are contained vs. cross, and alter the coordinates of those polygons so that they don't exceed the box. Or something like gIntersection. But surely there's a simpler way?

While I know that there are many problems with bounding boxes, and that a spatial overlay to a polygon that defines the region of interest is typically preferable, in many situations, bounding boxes work fine and are simpler.

  • Just to be clear, if your Spatial object is extended (Polygons or Lines) you want to cut it such that it returns only the chunk of it that is inside your given extent? I don't think there is a simpler way.
    – Spacedman
    Commented Jan 13, 2013 at 15:57
  • @Spacedman Clarified that I'm interested in both but the simpler version would suffice for the present question. Commented Jan 13, 2013 at 20:59
  • Have you already implemented the solution to (2) using rgeos? It sounds like you have at least tried. Could you give us that solution and an example so at least we have something to compare against for 'simplicity'? Because, to be honest, that seems pretty simple.
    – Spacedman
    Commented Jan 13, 2013 at 21:18
  • @Spacedman Everything's simple; just takes time.... :-) I tried it with gIntersection and came up with Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, "rgeos_intersection") : TopologyException: no outgoing dirEdge found at 3 2.5 No time to debug today; wrote up a sloppy version and will fix in the future. Commented Jan 14, 2013 at 0:22

2 Answers 2


I've created a small function for this very purpose and it has been used by others with good reviews!

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 = TRUE)

This should solve your problem. Further explanation is here: http://robinlovelace.net/r/2014/07/29/clipping-with-r.html

The dummy polygon b_poly that is created has no proj4 string, which results in "Warning: spgeom1 and spgeom2 have different proj4 strings", but this is harmless.

  • I have sp, maptools, rgdal, and rgeos loaded. I get Error in .class1(object) : could not find function "extent" R/package version issue perhaps? Commented Dec 10, 2014 at 21:05
  • Note the line library(raster) in my tutorial:robinlovelace.net/r/2014/07/29/clipping-with-r.html let us know how you get on! Cheers. Commented Dec 11, 2014 at 1:50
  • This produces a warning message for me: spgeom1 and spgeom2 have different proj4 strings. Adding proj4string(b_poly) <- proj4string(shp) should solve it?
    – Matifou
    Commented Mar 1, 2017 at 2:44

Here's a sloppy boundary version (sufficient to meet my needs in time for the mini-deadline tomorrow :-) ):

#' Convert a bounding box to a SpatialPolygons object
#' Bounding box is first created (in lat/lon) then projected if specified
#' @param bbox Bounding box: a 2x2 numerical matrix of lat/lon coordinates
#' @param proj4stringFrom Projection string for the current bbox coordinates (defaults to lat/lon, WGS84)
#' @param proj4stringTo Projection string, or NULL to not project
#' @seealso \code{\link{clipToExtent}} which uses the output of this to clip to a bounding box
#' @return A SpatialPolygons object of the bounding box
#' @example 
#' bb <- matrix(c(3,2,5,4),nrow=2)
#' rownames(bb) <- c("lon","lat")
#' colnames(bb) <- c('min','max')
as.SpatialPolygons.bbox <- function( bbox, proj4stringFrom=CRS("+proj=longlat +datum=WGS84"), proj4stringTo=NULL ) {
  # Create unprojected bbox as spatial object
  bboxMat <- rbind( c(bbox['lon','min'],bbox['lat','min']), c(bbox['lon','min'],bbox['lat','max']), c(bbox['lon','max'],bbox['lat','max']), c(bbox['lon','max'],bbox['lat','min']), c(bbox['lon','min'],bbox['lat','min']) ) # clockwise, 5 points to close it
  bboxSP <- SpatialPolygons( list(Polygons(list(Polygon(bboxMat)),"bbox")), proj4string=proj4stringFrom  )
  if(!is.null(proj4stringTo)) {
    bboxSP <- spTransform( bboxSP, proj4stringTo )

#' Restrict to extent of a polygon
#' Currently does the sloppy thing and returns any object that has any area inside the extent polygon
#' @param sp Spatial object
#' @param extent a SpatialPolygons object - any part of sp not within a polygon will be discarded
#' @seealso \code{\link{as.SpatialPolygons.bbox}} to create a SP from a bbox
#' @return A spatial object of the same type
#' @example
#' set.seed(1)
#' P4S.latlon <- CRS("+proj=longlat +datum=WGS84")
#' ply <- SpatialPolygons(list(Polygons(list(Polygon(cbind(c(2,4,4,1,2),c(2,3,5,4,2)))), "s1"),Polygons(list(Polygon(cbind(c(5,4,2,5),c(2,3,2,2)))), "s2")), proj4string=P4S.latlon)
#' pnt <- SpatialPoints( matrix(rnorm(100),ncol=2)+2, proj4string=P4S.latlon )
#' # Make bounding box as Spatial Polygon
#' bb <- matrix(c(3,2,5,4),nrow=2)
#' rownames(bb) <- c("lon","lat")
#' colnames(bb) <- c('min','max')
#' bbSP <- as.SpatialPolygons.bbox(bb, proj4stringTo=P4S.latlon )
#' # Clip to extent
#' plyClip <- clipToExtent( ply, bbSP )
#' pntClip <- clipToExtent( pnt, bbSP )
#' # Plot
#' plot( ply )
#' plot( pnt, add=TRUE )
#' plot( bbSP, add=TRUE, border="blue" )
#' plot( plyClip, add=TRUE, border="red")
#' plot( pntClip, add=TRUE, col="red", pch="o")
clipToExtent <- function( sp, extent ) {
    keep <- gContains( extent, sp,byid=TRUE ) | gOverlaps( extent, sp,byid=TRUE )
    stopifnot( ncol(keep)==1 )

bbox clipping

If you need the bounding box to project, the version here adds an interpolate argument, so that the resulting projected box is curved.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.