r - 如何在 R 中用多边形剪辑 WorldMap？

Question

我从 www.GADM.org 导入了一个世界地图数据集使用 R 包栅格。我想将它剪辑到我创建的多边形以减小 map 的大小。我可以检索数据并且可以毫无问题地创建多边形，但是当我使用“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)

错误信息:

Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, "rgeos_intersection") : Geometry collections may not contain other geometry collections In addition: Warning message: In RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, "rgeos_intersection") : spgeom1 and spgeom2 have different proj4 strings

Answer 1

您不需要使用 PBS(我已经学会了这一点，因为 @flowla 发布的 r-sig-geo 链接是我最初发布的问题!)。这段代码展示了如何在 rgeos 中完成这一切，感谢 various不同postings从罗杰比万德。这将是更规范的子集化方式，无需强制转换为 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，这些文件会保持最新状态，并且会不断检查错误(因为它们是开源的，如果您发现错误，请通过电子邮件发送给他们)。国家边界位于文化按钮下方。您会看到它们也更轻量级，您可以选择适合您需要的分辨率。