Klip objek spasial ke kotak pembatas di R

Diberikan objek Spasial dalam R, bagaimana cara klip semua elemennya berada di dalam kotak pembatas?

Ada dua hal yang ingin saya lakukan (idealnya saya tahu bagaimana melakukan keduanya, tetapi keduanya merupakan solusi yang dapat diterima untuk masalah saya saat ini - membatasi polifon shapefile ke benua AS).

1. Jatuhkan setiap elemen tidak sepenuhnya dalam kotak pembatas. Ini sepertinya bbox()<-cara yang logis, tetapi tidak ada metode seperti itu.

2. Lakukan operasi klip yang benar, sehingga elemen non-infinitesimal (mis. Poligon, garis) terpotong pada batas . sp::bboxtidak memiliki metode penugasan, jadi satu-satunya cara saya datang dengan menggunakan overatau gContains/ gCrossesdalam hubungannya dengan objek SpatialPolygons yang berisi kotak dengan koordinat kotak lompatan baru. Kemudian ketika memotong objek poligon, Anda harus mencari tahu mana yang berisi vs silang, dan mengubah koordinat poligon tersebut sehingga mereka tidak melebihi kotak. Atau semacamnya gIntersection. Tapi pasti ada cara yang lebih sederhana?

Meskipun saya tahu bahwa ada banyak masalah dengan kotak pembatas , dan bahwa pelapisan spasial ke poligon yang mendefinisikan wilayah yang diminati biasanya lebih disukai, dalam banyak situasi, kotak pembatas berfungsi dengan baik dan lebih sederhana.

Saya telah membuat fungsi kecil untuk tujuan ini dan telah digunakan oleh orang lain dengan ulasan yang bagus!

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

Ini harus menyelesaikan masalah Anda. Penjelasan lebih lanjut ada di sini: http://robinlovelace.net/r/2014/07/29/clipping-with-r.html

Polygon dummy b_polyyang dibuat tidak memiliki string proj4, yang menghasilkan " Peringatan: spgeom1 dan spgeom2 memiliki string proj4 yang berbeda ", tetapi ini tidak berbahaya.

Berikut ini adalah versi batas yang ceroboh (cukup untuk memenuhi kebutuhan saya tepat waktu untuk tenggat waktu mini besok :-)):

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

Jika Anda membutuhkan kotak pembatas untuk diproyeksikan, versi di sini menambahkan interpolateargumen, sehingga kotak yang diproyeksikan dihasilkan melengkung.