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.


Ż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.


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.

Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.