Clip de objetos Espaciales para el cuadro delimitador en R

Question

Dado un objeto Espacial en R, ¿cómo puedo clip de todos sus elementos se encuentran dentro de un cuadro delimitador?

Hay dos cosas que me gustaría hacer (idealmente me gustaría saber cómo hacer, pero tampoco es una solución aceptable a mi actual problema de la restricción de un shapefile de polígonos a los Estados Unidos continentales).

1. Colocar cada elemento no está completamente dentro del cuadro delimitador. Esto parece como bbox()<- sería la forma lógica, pero no existe tal método.

2. Hacer una verdadera clip de operación, de tal manera que no infinitesimal de elementos (por ejemplo, polígonos, líneas) se cortan en el límite. sp::bbox carece de un método de asignación, por lo que la única manera que he encontrado sería el uso de over o gContains/gCrosses en conjunción con un SpatialPolygons objeto que contiene un cuadro con la nueva delimitación del cuadro de coordenadas. Luego, cuando el recorte de un objeto poligonal, tendrías que averiguar cuáles son los contenidos vs cruz, y modificar las coordenadas de los polígonos para que no se exceda de la caja. O algo como gIntersection. Pero seguro que hay una manera más sencilla?

Aunque sé que hay muchos problemas con las cajas de contorno, y que una superposición espacial a un polígono que define la región de interés es generalmente preferible, en muchas situaciones, las cajas de contorno funcionan bien y son más simples.

Answer 1

He creado una pequeña función para este fin y que ha sido utilizado por otros con buenas críticas!

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

Esto debería resolver su problema. Otra explicación es la siguiente: http://robinlovelace.net/r/2014/07/29/clipping-with-r.html

Answer 2

He aquí un descuidado límite de versión (cantidad suficiente para satisfacer mis necesidades en el momento de la mini-plazo de mañana :-) ):

#' 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),]
}

Si usted necesita el cuadro delimitador del proyecto, la versión que aquí se añade un interpolate argumento, de modo que el resultado proyectado cuadro es curvo.