SpatialPolygonsDataFrame 的裁剪答案

【问题标题】：Crop for SpatialPolygonsDataFrameSpatialPolygonsDataFrame 的裁剪
【发布时间】：2012-12-08 14:22:30
【问题描述】：

我有两个SpatialPolygonsDataFrame 文件：dat1、dat2

extent(dat1)
class       : Extent 
xmin        : -180 
xmax        : 180 
ymin        : -90 
ymax        : 90 


extent(dat2)
class       : Extent 
xmin        : -120.0014 
xmax        : -109.9997 
ymin        : 48.99944 
ymax        : 60

我想使用 dat2 的范围裁剪文件 dat1。我不知道该怎么做。我之前只是使用“crop”功能处理光栅文件。

当我对我当前的数据使用这个函数时，会出现以下错误：

> r1 <- crop(BiomassCarbon.shp,alberta.shp)
Error in function (classes, fdef, mtable)  : 

 unable to find an inherited method for function ‘crop’ for signature"SpatialPolygonsDataFrame"’

【问题讨论】：

这是没有希望的，请处理涉及 dat1 和 dat2 或其他问题的问题的详细信息

标签： r crop spatial

【解决方案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")

【讨论】：

感谢您提供代码示例。我只是没有时间发布我的回复。 +1 为您使用 raster 包创建边界多边形的巧妙解决方案。
@JeffreyEvans -- 是的。我发现 raster 在为用户提供用户友好的转换功能/方法方面比 sp 做得更好。我很欣赏 sp 作者可能希望他们的包保持备用和轻量级（因为它旨在支撑许多其他包），但我仍然希望他们最终添加更多实用功能。
整洁：as(extent(130, 180, 40, 70), "SpatialPolygons") 谢谢！
crop 很棒 - 但是光栅中是否有反向裁剪可以做同样的事情，但与第二个 shapefile 相反？
@jebyrnes 查看?raster::erase 中的示例，了解其中的一种方法。

【解决方案2】：

这里是一个如何使用rgeos 以世界地图为例的示例

这来自R-sig-Geo mailing list 上的 Roger Bivand。 Roger 是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")

【讨论】：

很酷的答案。 gIntersection() 比 gIntersects() 更简单。

【解决方案3】：

您不能对 sp 多边形对象使用裁剪。您需要创建一个表示 dat2 的 bbox 坐标的多边形，然后可以使用 rgeos 库中的 gIntersects。

编辑：此评论与 2012 年可用的版本有关，现在不再如此。

【讨论】：

【解决方案4】：

见 ?crop

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

x 光栅* 对象

y Extent 对象，或任何可以从中获取 Extent 对象的对象提取（见详情

您需要使用光栅包中的rasterize 函数对第一个 SpatialPolygon 进行光栅化

我创建了一些数据来展示如何使用光栅化：

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"’

现在使用光栅化

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)

【讨论】：