Przypnij obiekt przestrzenny do ramki granicznej w R.

14

Biorąc pod uwagę obiekt przestrzenny w R, w jaki sposób mogę przyciąć wszystkie jego elementy, aby leżały w obwiedni?

Są dwie rzeczy, które chciałbym zrobić (najlepiej wiedziałbym, jak to zrobić, ale jedno lub drugie jest akceptowalnym rozwiązaniem mojego obecnego problemu - ograniczenie kształtu pliku wielokąta do kontynentalnych Stanów Zjednoczonych).

  1. Upuść każdy element nie do końca w obwiedni. Wydaje się, że bbox()<-byłby to logiczny sposób, ale taka metoda nie istnieje.

  2. Wykonaj prawdziwą operację przycinania, tak aby elementy nieskończenie małe (np. Wielokąty, linie) zostały odcięte na granicy . sp::bboxnie ma metody przypisania, więc jedynym sposobem, jaki wymyśliłem, byłoby użycie overlub gContains/ gCrossesw połączeniu z obiektem SpatialPolygons zawierającym pole ze współrzędnymi nowego obwiedni. Następnie, gdy wycinasz obiekt wielokąta, musisz dowiedzieć się, które są zawarte w stosunku do krzyża, i zmienić współrzędne tych wielokątów, aby nie przekraczały pola. Lub coś w tym rodzaju gIntersection. Ale z pewnością istnieje prostszy sposób?

Chociaż wiem, że istnieje wiele problemów z ramkami ograniczającymi i że przestrzenna nakładka na wielokąt, który określa region będący przedmiotem zainteresowania, jest zwykle lepsza, w wielu sytuacjach ramki ograniczające działają dobrze i są prostsze.

Ari B. Friedman
źródło
Żeby było jasne, jeśli twój obiekt przestrzenny jest rozszerzony (wielokąty lub linie), chcesz go wyciąć tak, aby zwracał tylko jego fragment znajdujący się w danym zasięgu? Nie sądzę, że istnieje prostszy sposób.
Spacedman
@Spacedman wyjaśnił, że jestem zainteresowany obydwoma, ale prostsza wersja wystarczyłaby do tego pytania.
Ari B. Friedman
Czy już zaimplementowałeś rozwiązanie (2) za pomocą rgeos? Wygląda na to, że przynajmniej próbowałeś. Czy możesz podać nam to rozwiązanie i przykład, abyśmy mieli przynajmniej coś do porównania pod względem „prostoty”? Ponieważ, szczerze mówiąc, wydaje się to dość proste.
Spacedman
@Spacedman Wszystko jest proste; po prostu wymaga czasu .... :-) Próbowałem gIntersectioni wymyśliłem Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, "rgeos_intersection") : TopologyException: no outgoing dirEdge found at 3 2.5 Nie ma czasu na debugowanie dzisiaj; napisałem niechlujną wersję i naprawię ją w przyszłości.
Ari B. Friedman

Odpowiedzi:

11

W tym celu stworzyłem małą funkcję, która była używana przez innych z dobrymi recenzjami!

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

To powinno rozwiązać twój problem. Dalsze wyjaśnienia znajdują się tutaj: http://robinlovelace.net/r/2014/07/29/clipping-with-r.html

Tworzony fikcyjny wielokąt b_polynie ma łańcucha proj4, co powoduje „ Ostrzeżenie: spgeom1 i spgeom2 mają różne łańcuchy proj4 ”, ale jest to nieszkodliwe.

RobinLovelace
źródło
Ja sp, maptools, rgdal, i rgeoszaładowany. Uzyskać Error in .class1(object) : could not find function "extent"wersji R / Pakiet problem może?
gregmacfarlane
Zwróć uwagę na wiersz library(raster)w moim tutorialu: robinlovelace.net/r/2014/07/29/clipping-with-r.html daj nam znać, jak się masz! Twoje zdrowie.
RobinLovelace,
To generuje komunikat ostrzegawczy: spgeom1 i spgeom2 mają różne ciągi proj4. Dodanie proj4string (b_poly) <- proj4string (shp) powinno go rozwiązać?
Matifou,
7

Oto niechlujna wersja brzegowa (wystarczająca, aby zaspokoić moje potrzeby na jutrzejszy mini-termin :-)):

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

obcinanie skrzynki

Jeśli potrzebujesz projektu ramki granicznej do projektu, tutaj wersja dodaje interpolateargument, dzięki czemu powstałe pole rzutowania jest zakrzywione.

Ari B. Friedman
źródło