我导入了世界地图数据集www.GADM.org使用 R 包栅格。我想将其剪切到我创建的多边形中以减小地图的大小。我可以检索数据并且可以毫无问题地创建多边形，但是当我使用“gIntersection”命令时，我收到一条晦涩的错误消息。
关于如何剪辑我的世界地图数据集有什么建议吗？
library(raster)
library(rgeos)
## Download Map of the World ##
WorldMap <- getData('countries')
## Create the clipping polygon
clip.extent <- as(extent(-20, 40, 30, 72), "SpatialPolygons")
proj4string(clip.extent) <- CRS(proj4string(WorldMap))
## Clip the map
EuropeMap <- gIntersection(WorldMap, clip.extent, byid = TRUE)
错误信息：
RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, "rgeos_intersection") 中的错误： 几何集合不能包含其他几何集合 另外：警告消息： 在 RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, "rgeos_intersection") 中： spgeom1 和 spgeom2 具有不同的 proj4 字符串
你不需要使用 PBS（我是通过惨痛的教训才知道这一点的，因为r-sig-geo@flowla 发布的链接是我最初发布的问题！）。这段代码展示了如何在 rgeos 中完成这一切，感谢various不同的postings来自罗杰·比万德。这将是更规范的子集化方式，无需强制 PolySet 对象。
出现错误消息的原因是您无法对 SpatialPolygon 集合进行 gIntersection，您需要单独执行它们。找出您想要使用哪些gIntersects。然后我使用每个国家多边形的子集gIntersection。我在将 SpatialPolygons 对象列表传递回 SpatialPolygons 时遇到问题，这会将裁剪后的 shapefile 转换为 SpatialPolygons，这是因为并非所有裁剪后的对象都是class SpatialPolygons。一旦我们排除了这些，一切就正常了。
# This very quick method keeps whole countries
gI <- gIntersects(WorldMap, clip.extent, byid = TRUE )
Europe <- WorldMap[which(gI), ]
plot(Europe)
#If you want to crop the country boundaries, it's slightly more involved:
# This crops countries to your bounding box
gI <- gIntersects(WorldMap, clip.extent, byid = TRUE)
out <- lapply(which(gI), function(x){
gIntersection(WorldMap[x,], clip.extent)
})
# But let's look at what is returned
table(sapply(out, class))
# SpatialCollections SpatialPolygons
# 2 63
# We want to keep only objects of class SpatialPolygons
keep <- sapply(out, class)
out <- out[keep == "SpatialPolygons"]
# Coerce list back to SpatialPolygons object
Europe <- SpatialPolygons(lapply(1:length(out), function(i) {
Pol <- slot(out[[i]], "polygons")[[1]]
slot(Pol, "ID") <- as.character(i)
Pol
}))
plot(Europe)
如果可以的话，我建议你看看自然地球数据。他们拥有保持最新的高质量形状文件，并不断检查错误（因为它们是开源的，如果您发现错误，请通过电子邮件发送给他们）。国家边界位于Cultural纽扣。您会发现它们也更轻量级，您可以选择适合您需求的分辨率。