您不需要使用 PBS(我已经学会了这一点,因为@flowla 发布的r-sig-geo链接是我最初发布的问题!)。这段代码展示了如何在 rgeos 中完成这一切,这要归功于Roger Bivand 的各种不同帖子。这将是更规范的子集化方式,无需强制转换为 PolySet 对象。
错误消息的原因是您不能对 SpatialPolygons 集合进行 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)
如果可以,我建议您查看naturalearthdata。他们拥有高质量的 shapefile,这些文件会保持最新状态,并且会不断检查错误(因为它们是开源的,如果您发现错误,请通过电子邮件发送给他们)。国家边界位于文化按钮下方。您会看到它们也更轻量级,您可以选择适合您需要的分辨率。