Clip Spatial-Objekt an Begrenzungsrahmen in R

13

Wie kann ich bei einem Raumobjekt in R alle seine Elemente so beschneiden, dass sie in einem Begrenzungsrahmen liegen?

Ich möchte zwei Dinge tun (im Idealfall weiß ich, wie man beides macht, aber beides ist eine akzeptable Lösung für mein aktuelles Problem - die Beschränkung eines Polygon-Shapefiles auf die kontinentalen USA).

  1. Legen Sie jedes Element nicht vollständig im Begrenzungsrahmen ab. Dies scheint bbox()<-der logische Weg zu sein, aber es gibt keine solche Methode.

  2. Führen Sie eine echte Clip-Operation durch, sodass nicht infinitesimale Elemente (z. B. Polygone, Linien) an der Grenze abgeschnitten werden . sp::bboxDa mir eine Zuweisungsmethode fehlt, ist die einzige Möglichkeit, die ich mir ausgedacht habe, die Verwendung overoder gContains/ gCrossesund die Verknüpfung mit einem SpatialPolygons-Objekt, das eine Box mit den Koordinaten der neuen Begrenzungsbox enthält. Wenn Sie dann ein Polygonobjekt beschneiden, müssen Sie herausfinden, welche Polygone im Vergleich zu Kreuz enthalten sind, und die Koordinaten dieser Polygone so ändern, dass sie das Feld nicht überschreiten. Oder so ähnlich gIntersection. Aber es gibt doch einen einfacheren Weg?

Obwohl ich weiß, dass es viele Probleme mit Begrenzungsrahmen gibt und dass eine räumliche Überlagerung eines Polygons, das den interessierenden Bereich definiert, in der Regel vorzuziehen ist, funktionieren Begrenzungsrahmen in vielen Situationen einwandfrei und sind einfacher.

Ari B. Friedman
quelle
Um ganz klar zu sein, wenn Ihr räumliches Objekt erweitert ist (Polygone oder Linien), möchten Sie es so schneiden, dass es nur den Teil davon zurückgibt, der sich in Ihrer vorgegebenen Ausdehnung befindet? Ich glaube nicht, dass es einen einfacheren Weg gibt.
Spacedman
@Spacedman Klargestellt, dass ich an beiden interessiert bin, aber die einfachere Version würde für die vorliegende Frage ausreichen.
Ari B. Friedman
Haben Sie die Lösung bereits implementiert, um (2) Rgeos zu verwenden? Es hört sich so an, als hättest du es zumindest versucht. Können Sie uns diese Lösung und ein Beispiel geben, damit wir zumindest der Einfachheit halber etwas Vergleichbares haben? Denn um ehrlich zu sein, scheint das ziemlich einfach zu sein.
Spacedman
@Spacedman Alles ist einfach; braucht einfach Zeit .... :-) Ich habe es mit probiert gIntersectionund hatte Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, "rgeos_intersection") : TopologyException: no outgoing dirEdge found at 3 2.5 heute keine Zeit zum Debuggen; schrieb eine schlampige Version und wird in Zukunft beheben.
Ari B. Friedman

Antworten:

11

Ich habe zu diesem Zweck eine kleine Funktion erstellt, die von anderen mit guten Rezensionen verwendet wurde!

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

Dies sollte Ihr Problem lösen. Weitere Erklärungen finden Sie hier: http://robinlovelace.net/r/2014/07/29/clipping-with-r.html

Das erstellte Dummy-Polygon b_polyhat keine proj4-Zeichenfolge, was zu der Meldung " Warnung: spgeom1 und spgeom2 haben unterschiedliche proj4-Zeichenfolgen " führt. Dies ist jedoch harmlos.

RobinLovelace
quelle
Ich habe sp, maptools, rgdal, und rgeosgeladen. Ich habe Error in .class1(object) : could not find function "extent"vielleicht ein Problem mit der R / Paket-Version.
Gregmacfarlane
Beachten Sie die Zeile library(raster)in meinem Tutorial: robinlovelace.net/r/2014/07/29/clipping-with-r.html Lassen Sie uns wissen, wie Sie vorankommen ! Prost.
Robin Lovelace
Dies erzeugt eine Warnmeldung für mich: spgeom1 und spgeom2 haben unterschiedliche proj4-Zeichenfolgen. Das Hinzufügen von proj4string (b_poly) <- proj4string (shp) sollte es lösen?
Matifou
7

Hier ist eine schlampige Boundary-Version (ausreichend, um meine Bedürfnisse rechtzeitig für die Mini-Deadline von morgen zu erfüllen :-)):

#' 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 Ausschnitt

Wenn Sie den Begrenzungsrahmen zum Projizieren benötigen, fügt die Version hier ein interpolateArgument hinzu, sodass der resultierende projizierte Rahmen gekrümmt ist.

Ari B. Friedman
quelle