Прикрепить пространственный объект к ограничительной рамке в R

13

Для данного Пространственного объекта в R, как я могу обрезать все его элементы, чтобы они лежали внутри ограничительной рамки?

Я хотел бы сделать две вещи (в идеале я бы знал, как сделать обе, но любая из них является приемлемым решением моей текущей проблемы - ограничение шейп-файла полигона для континентальной части США).

  1. Удалите каждый элемент не полностью в пределах ограничительной рамки. Это кажется bbox()<-логичным, но такого метода не существует.

  2. Выполните настоящую операцию обрезки, чтобы на границе были обрезаны элементы, не являющиеся бесконечно малыми (например, многоугольники, линии) . sp::bboxотсутствует метод присваивания, поэтому единственный способ, которым я придумал, - это использовать overили gContains/ gCrossesв сочетании с объектом SpatialPolygons, содержащим прямоугольник с координатами нового ограничивающего прямоугольника. Затем при обрезке объекта многоугольника вы должны выяснить, какие из них содержатся в сравнении с крестом, и изменить координаты этих многоугольников, чтобы они не превышали рамки. Или что-то вроде gIntersection. Но наверняка есть более простой способ?

Хотя я знаю, что существует много проблем с ограничивающими прямоугольниками и что пространственное наложение на многоугольник, определяющий интересующую область, обычно предпочтительнее, во многих ситуациях ограничивающие прямоугольники работают нормально и проще.

Ари Б. Фридман
источник
Просто чтобы прояснить: если ваш пространственный объект расширен (полигоны или линии), вы хотите разрезать его так, чтобы он возвращал только тот кусок, который находится внутри вашего заданного экстента? Я не думаю, что есть более простой способ.
Spacedman
@Spacedman Уточнил, что я интересуюсь обоими, но более простой версии будет достаточно для настоящего вопроса.
Ари Б. Фридман
Вы уже реализовали решение (2), используя rgeos? Похоже, вы хотя бы пытались. Не могли бы вы дать нам это решение и пример, чтобы хотя бы у нас было что сравнить с «простотой»? Потому что, если честно, это кажется довольно простым.
Spacedman
@Spacedman Все просто; просто занимает время .... :-) Я попробовал это gIntersectionи придумал Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, "rgeos_intersection") : TopologyException: no outgoing dirEdge found at 3 2.5 Нет времени для отладки сегодня; написал небрежную версию и исправлю в будущем.
Ари Б. Фридман

Ответы:

11

Я создал небольшую функцию для этой цели, и она была использована другими с хорошими отзывами!

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

Это должно решить вашу проблему. Дальнейшее объяснение здесь: http://robinlovelace.net/r/2014/07/29/clipping-with-r.html

Созданный фиктивный многоугольник b_polyне имеет строки proj4, что приводит к « Предупреждение: spgeom1 и spgeom2 имеют разные строки proj4 », но это безвредно.

RobinLovelace
источник
Я sp, maptools, rgdalи rgeosзагружен. Error in .class1(object) : could not find function "extent"Возможно, у меня проблема с версией R / пакета?
gregmacfarlane
Обратите внимание на строку library(raster)в моем уроке: robinlovelace.net/r/2014/07/29/clipping-with-r.html дайте нам знать, как вы ладите ! Приветствия.
RobinLovelace
Это выдает мне предупреждение: spgeom1 и spgeom2 имеют разные строки proj4. Добавление proj4string (b_poly) <- proj4string (shp) должно решить это?
Матифу
7

Вот неаккуратная граничная версия (достаточная для того, чтобы успеть вовремя к завтрашнему мини-сроку :-)):

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

Если вам нужен ограничивающий прямоугольник для проекта, версия здесь добавляет interpolateаргумент, так что в результате проецируются окно изогнут.

Ари Б. Фридман
источник