Recorte el objeto espacial al cuadro delimitador en R

13

Dado un objeto espacial en R, ¿cómo puedo recortar todos sus elementos para que se encuentren dentro de un cuadro delimitador?

Hay dos cosas que me gustaría hacer (lo ideal sería saber hacer ambas cosas, pero cualquiera de las dos es una solución aceptable para mi problema actual: restringir un archivo de forma poligonal a los EE. UU. continental).

  1. Elimine cada elemento que no esté completamente dentro del cuadro delimitador. Parece que bbox()<- sería la forma lógica, pero no existe tal método.

  2. Realice una verdadera operación de recorte, de manera que los elementos no infinitesimales (por ejemplo, polígonos, líneas) se eliminen en el límite . sp::bbox carece de un método de asignación, por lo que la única forma de hacerlo sería utilizar over o gContains / gCrosses junto con un objeto SpatialPolygons que contenga un cuadro con las coordenadas del nuevo cuadro delimitador. Luego, al recortar un objeto de polígono, tendrías que averiguar cuáles están contenidos frente a cruz, y alterar las coordenadas de esos polígonos para que no excedan la casilla. O algo como gIntersection . Pero seguramente hay una forma más simple?

Aunque sé que hay muchos problemas con los cuadros delimitadores , y que la superposición espacial a un polígono que define la región de interés suele ser preferible, en muchas situaciones, los cuadros de delimitación funcionan bien y son más simples.

    
pregunta Ari B. Friedman 13.01.2013 - 13:48

2 respuestas

11

He creado una pequeña función para este propósito y ha sido utilizada 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 = TRUE)
}

Esto debería resolver su problema. La explicación adicional está aquí: enlace

El polígono ficticio b_poly que se creó no tiene una cadena proj4, lo que da como resultado " Advertencia: spgeom1 y spgeom2 tienen diferentes cadenas proj4 ", pero esto es inofensivo.

    
respondido por el RobinLovelace 06.08.2014 - 17:12
6

Aquí hay una versión de límite descuidada (suficiente para satisfacer mis necesidades a tiempo para la fecha límite mínima 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),]
}

Sinecesitaelcuadrodelimitadorparaproyectar,laversión aquí agrega un argumento interpolate , de modo que el cuadro proyectado resultante está curvado.

    
respondido por el Ari B. Friedman 14.01.2013 - 01:21

Lea otras preguntas en las etiquetas