Обрезка для SpatialPolygonsDataFrame

Ответ 1

Добавлен оптимизированный метод 2014-10-9:

raster::crop() можно использовать для обрезки объектов Spatial* (а также Raster*).

Например, здесь, как вы можете использовать его для обрезки объекта SpatialPolygons*:

## Load raster package and an example SpatialPolygonsDataFrame
library(raster) 
data("wrld_simpl", package="maptools")

## Crop to the desired extent, then plot
out <- crop(wrld_simpl, extent(130, 180, 40, 70))
plot(out, col="khaki", bg="azure2")

Ответ на исходный (и все еще функциональный):

Функция rgeos gIntersection() делает это довольно простым.

Используя отличный пример mnel в качестве точки сбрасывания:

library(maptools)
library(raster)   ## To convert an "Extent" object to a "SpatialPolygons" object.
library(rgeos)
data(wrld_simpl)

## Create the clipping polygon
CP <- as(extent(130, 180, 40, 70), "SpatialPolygons")
proj4string(CP) <- CRS(proj4string(wrld_simpl))

## Clip the map
out <- gIntersection(wrld_simpl, CP, byid=TRUE)

## Plot the output
plot(out, col="khaki", bg="azure2")

Ответ 2

Вот пример того, как это сделать с помощью rgeos с использованием карты мира в качестве примера

Это происходит от Роджера Биванда в список рассылки R-sig-Geo. Роджер является одним из авторов пакета sp.

Использование карты мира в качестве примера

library(maptools)

data(wrld_simpl)

# interested in the arealy bounded by the following rectangle
# rect(130, 40, 180, 70)

library(rgeos)
# create  a polygon that defines the boundary
bnds <- cbind(x=c(130, 130, 180, 180, 130), y=c(40, 70, 70, 40, 40))
# convert to a spatial polygons object with the same CRS
SP <- SpatialPolygons(list(Polygons(list(Polygon(bnds)), "1")),
proj4string=CRS(proj4string(wrld_simpl)))
# find the intersection with the original SPDF
gI <- gIntersects(wrld_simpl, SP, byid=TRUE)
# create the new spatial polygons object.
out <- vector(mode="list", length=length(which(gI)))
ii <- 1
for (i in seq(along=gI)) if (gI[i]) {
  out[[ii]] <- gIntersection(wrld_simpl[i,], SP)
  row.names(out[[ii]]) <- row.names(wrld_simpl)[i]; ii <- ii+1
}
# use rbind.SpatialPolygons method to combine into a new object.
out1 <- do.call("rbind", out)
# look here is Eastern Russia and a bit of Japan and China.
plot(out1, col = "khaki", bg = "azure2")

Ответ 3

Вы не можете использовать обрезку для объектов poly polygon. Вам нужно будет создать многоугольник, представляющий координаты bbox dat2, а затем использовать gIntersects в библиотеке rgeos.

Изменить: этот комментарий был в отношении версии, доступной в 2012 году, и это уже не так.

Ответ 4

увидеть? crop

corp (x, y, filename = "", snap = 'near', datatype = NULL,...)

x Растровый объект

     

y Объект Extent или любой объект, из которого может быть объект Extent   извлечен (см. подробности

Вам нужно растрировать первый SpatialPolygon, используя rasterize функцию из растрового пакета

Я создаю некоторые данные, чтобы показать, как использовать rasterize:

n <- 1000
x <- runif(n) * 360 - 180
y <- runif(n) * 180 - 90
xy <- cbind(x, y)
vals <- 1:n
p1 <- data.frame(xy, name=vals)
p2 <- data.frame(xy, name=vals)
coordinates(p1) <- ~x+y
coordinates(p2) <- ~x+y

если я попробую:

 crop(p1,p2)
 unable to find an inherited method for function ‘crop’ for signature ‘"SpatialPointsDataFrame"’

Теперь, используя rasterize

r <- rasterize(p1, r, 'name', fun=min)
crop(r,p2)

class       : RasterLayer 
dimensions  : 18, 36, 648  (nrow, ncol, ncell)
resolution  : 10, 10  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 
data source : in memory
names       : layer 
values      : 1, 997  (min, max)