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.
我使用R包栅格从www.GADM.org导入了世界地图数据集。我想将它剪辑为我创建的多边形,以减小地图的大小。我可以检索数据,我可以创建多边形没有问题,但当我使用'gIntersection'命令时,我得到一个模糊的错误消息。
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 Message:
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
2 个解决方案
#1
8
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.
你不需要使用PBS(我已经学到了很多东西,因为@flowla发布的r-sig-geo链接是我最初发布的一个问题!)。此代码显示了如何在rgeos中完成所有操作,感谢Roger Bivand的各种不同帖子。这将是更典型的子集化方法,而不需要对PolySet对象进行强制。
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.
错误消息的原因是您无法对SpatialPolygons集合进行gIntersection,您需要单独执行它们。使用gIntersects找出你想要的。然后我使用gIntersection对每个国家多边形进行子集化。我将SpatialPolygons对象列表传递回SpatialPolygons时出现问题,后者将裁剪的shapefile转换为SpatialPolygons,这是因为并非所有裁剪对象都属于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 )
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.
如果可以,我建议你看看naturalearthdata。它们具有高质量的shapefile,并且不断检查错误(因为如果发现错误,它们是开源的,请发送电子邮件)。国家边界属于文化按钮。您会发现它们的重量也更轻,您可以选择适合您需求的分辨率。
#2
2
How about a little intermediate step? I adopted the following code mainly from R-sig-Geo and I think it should do the trick. You'll need both 'maptools' and 'PBSmapping' packages, so make sure you've got them installed. Here's my code:
一点中间步骤怎么样?我主要从R-sig-Geo采用以下代码,我认为应该这样做。你需要'maptools'和'PBSmapping'包,所以一定要安装它们。这是我的代码:
# Required packages
library(raster)
library(maptools)
library(PBSmapping)
# Download world map
WorldMap <- getData('countries')
# Convert SpatialPolygons to class 'PolySet'
WorldMap.ps <- SpatialPolygons2PolySet(WorldMap)
# Clip 'PolySet' by given extent
WorldMap.ps.clipped <- clipPolys(WorldMap.ps, xlim = c(-20, 40), ylim = c(30, 72))
# Convert clipped 'PolySet' back to SpatialPolygons
EuropeMap <- PolySet2SpatialPolygons(WorldMap.ps.clipped, close_polys=TRUE)
I just tested it and it worked without any problems. It took some computation time to transform SpatialPolygons to PolySet, though.
我只是测试它,它没有任何问题。但是,将SpatialPolygons转换为PolySet需要一些计算时间。
Cheers, Florian
#1
8
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.
你不需要使用PBS(我已经学到了很多东西,因为@flowla发布的r-sig-geo链接是我最初发布的一个问题!)。此代码显示了如何在rgeos中完成所有操作,感谢Roger Bivand的各种不同帖子。这将是更典型的子集化方法,而不需要对PolySet对象进行强制。
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.
错误消息的原因是您无法对SpatialPolygons集合进行gIntersection,您需要单独执行它们。使用gIntersects找出你想要的。然后我使用gIntersection对每个国家多边形进行子集化。我将SpatialPolygons对象列表传递回SpatialPolygons时出现问题,后者将裁剪的shapefile转换为SpatialPolygons,这是因为并非所有裁剪对象都属于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 )
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.
如果可以,我建议你看看naturalearthdata。它们具有高质量的shapefile,并且不断检查错误(因为如果发现错误,它们是开源的,请发送电子邮件)。国家边界属于文化按钮。您会发现它们的重量也更轻,您可以选择适合您需求的分辨率。
#2
2
How about a little intermediate step? I adopted the following code mainly from R-sig-Geo and I think it should do the trick. You'll need both 'maptools' and 'PBSmapping' packages, so make sure you've got them installed. Here's my code:
一点中间步骤怎么样?我主要从R-sig-Geo采用以下代码,我认为应该这样做。你需要'maptools'和'PBSmapping'包,所以一定要安装它们。这是我的代码:
# Required packages
library(raster)
library(maptools)
library(PBSmapping)
# Download world map
WorldMap <- getData('countries')
# Convert SpatialPolygons to class 'PolySet'
WorldMap.ps <- SpatialPolygons2PolySet(WorldMap)
# Clip 'PolySet' by given extent
WorldMap.ps.clipped <- clipPolys(WorldMap.ps, xlim = c(-20, 40), ylim = c(30, 72))
# Convert clipped 'PolySet' back to SpatialPolygons
EuropeMap <- PolySet2SpatialPolygons(WorldMap.ps.clipped, close_polys=TRUE)
I just tested it and it worked without any problems. It took some computation time to transform SpatialPolygons to PolySet, though.
我只是测试它,它没有任何问题。但是,将SpatialPolygons转换为PolySet需要一些计算时间。
Cheers, Florian