Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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.

share|improve this question
    
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 Jan 13 '13 at 15:57
    
@Spacedman Clarified that I'm interested in both but the simpler version would suffice for the present question. –  Ari B. Friedman Jan 13 '13 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 Jan 13 '13 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. –  Ari B. Friedman Jan 14 '13 at 0:22
add comment

1 Answer

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 )
  }
  bboxSP
}


#' 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 ) {
    require(rgeos)
    keep <- gContains( extent, sp,byid=TRUE ) | gOverlaps( extent, sp,byid=TRUE )
    stopifnot( ncol(keep)==1 )
    sp[drop(keep),]
}

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.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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