Fixer l'objet spatial au cadre de sélection en R

13

Étant donné un objet spatial dans R, comment puis-je découper tous ses éléments pour qu'ils se trouvent dans un cadre de sélection?

Il y a deux choses que j'aimerais faire (idéalement, je saurais faire les deux, mais l'une ou l'autre est une solution acceptable à mon problème actuel - restreindre un fichier de formes polygonales aux États-Unis continentaux).

  1. Déposez chaque élément pas complètement dans le cadre de sélection. Cela semble bbox()<-être la voie logique, mais aucune telle méthode n'existe.

  2. Effectuez une véritable opération de découpage, de sorte que les éléments non infinitésimaux (par exemple les polygones, les lignes) soient coupés à la frontière . sp::bboxne dispose pas d'une méthode d'affectation, donc la seule façon dont je suis venu serait d'utiliser overou gContains/ gCrossesen conjonction avec un objet SpatialPolygons contenant une boîte avec les coordonnées de la nouvelle boîte englobante. Ensuite, lors de la découpe d'un objet polygone, vous devez déterminer ceux qui sont contenus par rapport à la croix et modifier les coordonnées de ces polygones afin qu'ils ne dépassent pas la zone. Ou quelque chose comme ça gIntersection. Mais il y a sûrement un moyen plus simple?

Bien que je sache qu'il existe de nombreux problèmes avec les boîtes englobantes et qu'une superposition spatiale à un polygone qui définit la région d'intérêt est généralement préférable, dans de nombreuses situations, les boîtes englobantes fonctionnent bien et sont plus simples.

Ari B. Friedman
la source
Juste pour être clair, si votre objet spatial est étendu (polygones ou lignes), vous voulez le couper de telle sorte qu'il n'en retourne que la partie qui se trouve à l'intérieur de votre étendue donnée? Je ne pense pas qu'il existe un moyen plus simple.
Spacedman
@Spacedman Clarifié que je m'intéresse aux deux mais la version plus simple suffirait pour la présente question.
Ari B. Friedman
Avez-vous déjà implémenté la solution pour (2) à l'aide de rgeos? Il semble que vous ayez au moins essayé. Pourriez-vous nous donner cette solution et un exemple afin qu'au moins nous ayons quelque chose à comparer pour la «simplicité»? Parce que, pour être honnête, cela semble assez simple.
Spacedman
@Spacedman Tout est simple; prend juste du temps .... :-) J'ai essayé avec gIntersectionet j'ai trouvé Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, "rgeos_intersection") : TopologyException: no outgoing dirEdge found at 3 2.5 Pas de temps pour déboguer aujourd'hui; a écrit une version bâclée et corrigera dans le futur.
Ari B. Friedman

Réponses:

11

J'ai créé une petite fonction à cet effet et elle a été utilisée par d'autres avec de bonnes critiques!

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

Cela devrait résoudre votre problème. De plus amples explications sont disponibles ici: http://robinlovelace.net/r/2014/07/29/clipping-with-r.html

Le polygone factice b_polyqui est créé n'a pas de chaîne proj4, ce qui entraîne « Avertissement: spgeom1 et spgeom2 ont des chaînes proj4 différentes », mais cela est inoffensif.

RobinLovelace
la source
Je sp, maptools, rgdalet rgeoschargé. Je reçois Error in .class1(object) : could not find function "extent"problème version R / paquet peut-être?
gregmacfarlane
Notez la ligne library(raster)dans mon tutoriel: robinlovelace.net/r/2014/07/29/clipping-with-r.html faites-nous savoir comment vous vous débrouillez! À votre santé.
RobinLovelace
Cela produit un message d'avertissement pour moi: spgeom1 et spgeom2 ont des chaînes proj4 différentes. Ajouter proj4string (b_poly) <- proj4string (shp) devrait le résoudre?
Matifou
7

Voici une version limite bâclée (suffisante pour répondre à mes besoins à temps pour le mini-délai de demain :-)):

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

détourage bbox

Si vous avez besoin de la boîte englobante pour projeter, la version ajoute ici un interpolateargument, de sorte que la boîte projetée résultante soit courbée.

Ari B. Friedman
la source