如何在R中用多边形裁剪WorldMap? [英] How to clip WorldMap with polygon in R?

问题描述

我已经使用R包栅格从www.GADM.org导入了世界地图数据集.我想将其裁剪到我创建的多边形上以减小地图的大小.我可以检索数据，也可以创建多边形，但是当我使用"gIntersection"命令时，会收到模糊的错误消息.

I have imported a world map dataset from www.GADM.org using the R package raster. I would like to clip it to a polygon I create to reduce the size of the map. I can retrieve the data and I can create the polygon no problem, but when I use the 'gIntersection' command I get an obscure error message.

关于如何裁剪我的世界地图数据集的任何建议?

Any suggestions on how to clip my World Map dataset?

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

推荐答案

您不需要使用PBS(由于

You don't need to use PBS (I have learnt this the hard way, as the r-sig-geo link posted by @flowla was a question originally posted by me!). This code shows how to do it all in rgeos, with thanks to various different postings from Roger Bivand. This would be the more canonical way of subsetting, without resorting to coercion to PolySet objects.

出现此错误消息的原因是，您无法对SpatialPolygon的集合进行gIntersection，您需要单独进行处理.使用gIntersects找出想要的.然后，我使用gIntersection子集每个国家/地区多边形.我在将SpatialPolygons对象列表返回给SpatialPolygons时遇到问题，这将裁剪的shapefile转换为SpatialPolygons，这是因为并非所有裁剪的对象都属于class SpatialPolygons.一旦我们排除了这些，一切都会很好.

The reason for the error message is that you can't gIntersection on a collection of SpatialPolygons, you need to do them individually. Find out which ones you want using gIntersects. I then subset each country polygon using gIntersection. I had a problem passing a list of SpatialPolygons objects back to SpatialPolygons, which turns the croppped shapefiles into SpatialPolygons, which was because not all cropped objects were of class SpatialPolygons. Once we excluded these everything works fine.

# 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，这些文件始终保持最新状态，并不断检查是否存在错误(因为如果发现错误，它们是开源的，请通过电子邮件发送给他们).国家/地区边界位于 Cultural 按钮下.您会发现它们也更轻巧，您可以选择适合您需要的分辨率.

If you can, I recommend you look at naturalearthdata. They have high quality shapefiles that are kept up-to-date and are constantly checked for errors (since they are open source if you find an error do email them). Country boundaries are under the Cultural buttons. You will see they are also a bit more lightweight and you can choose a resolution appropriate for your needs.

这篇关于如何在R中用多边形裁剪WorldMap?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！