如何在R中剪切带有多边形的WorldMap？

Question

我使用R包栅格从www.GADM.org导入了世界地图数据集.我想将其剪切为我创建的多边形以减小地图的大小.我可以检索数据,我可以创建多边形没有问题,但当我使用'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中完成所有操作,这要归功于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,可以保持最新并不断检查错误(因为如果发现错误,它们是开源的,请发送电子邮件).国家边界属于文化按钮.您会发现它们的重量也更轻,您可以选择适合您需求的分辨率.